Hierarchical Clustering using Reversible Binary Cellular Automata for High-Dimensional Data

Usually, the real-life datasets to be clustered are in the form of a set of data elements (rows) where each data element has several features (columns) with each feature’s value as a real number. So, to apply CA, we first need to convert each of the data elements containing the feature (attribute) values into a binary string. For our work, we consider only the datasets with quantifiable feature properties such that they can be converted into binary strings, namely the datasets having only continuous and categorical attributes. However, the scheme will work on any other dataset provided the attribute values are intelligently encoded into binary numbers in the preprocessing stage.

This work considers frequency based encoding [6] for the preprocessing. The essence of this encoding scheme is it converts the real numbers into binary maintaining a minimum hamming distance between consecutive intervals. Here we consider that the range of values for the continuous attributes can be divided into a maximum of $4$ disjoint intervals, so two-bit encoding is required for such an attribute. Let $x$ be a feature value and ${\tt M}$ be the encoding function. Then

where [$x_{1},x_{s}$], [$x_{s+1},x_{w}$], [$x_{w+1},x_{q}$] and [$x_{q+1},x_{t}$] are the maximum four disjoint intervals for the feature. Similarly, for a categorical attribute with $k$ distinct values, $k$ bits are needed to maintain minimum hamming distance where $i^{th}$ unique value will have one $1$ in the $i^{th}$ position with all other $k-1$ bits as $0$s. Once the features are encoded, the feature values of every attribute of an object can be concatenated to get a binary string per object. These binary strings can now be used in the CA as configurations of the reversible CA. We name these configurations as target configurations.

As mentioned, we are considering $1$-dimensional $2$-state $2$-radius finite CAs for clustering. The reason for choosing radius (r) as $2$ is – as per our encoding, a minimum of two bits are needed to represent a feature. In the CA, each cell updates its state based on a local rule, also called a rule $\mathcal{R}:\{0,1\}^{5}\rightarrow\{0,1\}$ which takes the neighborhood cells’ values as arguments. Each of these combinations of arguments is also called a Rule Min Term or an RMT. It is often represented by the string generated by concatenating the arguments or its decimal equivalent. For example, the argument to $\mathcal{R}$ where all the neighbors have state $1$ is represented by RMT $11111$ or RMT $31$ which is the decimal equivalent of $11111$. Two RMTs can be grouped if, among the neighborhood combinations, all the neighbors except one have the same state. We can define it formally as follows:

As we are taking $2$-state $2$-radius CA, that is, number of neighbors is $5$, we can get five such different $\mathcal{E}^{k}$, $0\leq k\leq 4$ where each $\mathcal{E}^{k}$ contains 16 sets $\mathcal{E}^{k}_{i},~{}0\leq i\leq 15$. The groups are shown in Table 1.

A rule can be represented by the string generated by concatenating next state values of all RMTs starting from $31$ to $0$, that is, $\mathcal{R}(1,1,1,1,1)$ to $\mathcal{R}(0,0,0,0,0)$, or its decimal equivalent. For instance, rule $267422991$ represents the string $00001111111100001000110100001111$ which means, for this rule, $\mathcal{R}(0,0,0,0,0)=\mathcal{R}(0,0,0,0,1)=\mathcal{R}(0,0,0,1,0)=\mathcal{R}(0,0,0,1,1)=1$, $\mathcal{R}(0,0,1,0,0)=0$ and so on. A snapshot of the states of all cells at any time instant is called the configuration. The global transition function $G_{n}:\mathcal{C}_{n}\rightarrow\mathcal{C}_{n}$ works on the set of configurations $\mathcal{C}_{n}$: for any ${\tt x},{\tt y}\in\mathcal{C}_{n}$, if ${\tt y}=({y_{i}})_{\forall i\in n}$ is the next configuration of ${\tt x}=(x_{i})_{\forall i\in n}$, then ${\tt y}=G_{n}({\tt x})=G_{n}(x_{0}x_{1}\cdots x_{n-1})=(\mathcal{R}(x_{i-2},x_{i-1},x_{i},x_{i+1},x_{i+2}))_{0\leq i\leq n-1}$

For example, Figure 1 represents the cycles of a $5$-cell reversible CA $267422991$. Here, the configurations are represented by their corresponding equivalent decimal numbers. We can see that, configuration $8(01000)$ is reachable from the configuration $13(01101)$ but not from $7(00111)$. Now, let us consider a hypothetical dataset (shown in Table 2) with one continuous and one categorical attribute.

The categorical attribute has three distinct values, so we need three bits to encode it. Whereas, based on the values of the continuous attribute, we can create three disjoint sub-intervals $[5,7]$, $[10,15]$ and $[50,55]$ to be encoded by $00$, $01$, and $11$ respectively. In this way, each data element can be converted into a binary string which can be considered as the configurations of a $5$-cell CA. Observe that, the encoding function is surjective mapping more than one object into the same bits (see Object No. 5 (row 7) and 6 (row 8) of Table 2): that is, the encoding itself provides a basic clustering of data elements. Now, let us cluster these objects using the CA $267422991$ of Figure 1. As per the cycles of this CA, the configurations $1$, $9$, and $12$ are reachable from each other and form one cluster, whereas, $4,25$ & $26$ form another cluster and $2$ belong to a third cluster. So, by CA $267422991$, the objects of this hypothetical dataset are clustered as $(Obj_{1},Obj_{2},Obj_{5},Obj_{6})$, $(Obj_{3})$, $(Obj_{4},Obj_{7},Obj_{8})$.

Let us first briefly discuss some established algorithms that are to be used as a reference point of comparison in this paper. These algorithms are well-documented and work for both larger and smaller datasets. The most classical of these algorithms is $K$-means [9, 12]. Introduced in 1973 and inspiring researchers to come up with improvements since then, this seminal algorithm figures out $K$ number of center points in a data set and then groups each data point into the nearest centroid. The next popular algorithm is DBSCAN [8]. It works on the principle of density-based spatial clustering of applications with noise where the neighborhood within a given radius has to have at least a certain number of points for each point of a cluster. A contemporary algorithm to DBSCAN is BIRCH [19] which clusters huge datasets by first producing a compact summary of a large dataset while maintaining as much information available after which the clustering takes place. A relatively new algorithm is Meanshift [5]. This algorithm uses a pattern recognition procedure, mean shift to assign data points to clusters in an iterative manner while shifting the data points toward the modes of density. The most effective one is agglomerative hierarchical clustering, here abbreviated as hierarchical [3], which iteratively clusters the data points giving the best set of clusters.

Table 3 lists the standard benchmark validation indices that we are using for checking the performance of the clustering algorithms. Here, the Silhouette score is a real number within (-1, 1) where the higher the score, the better the clusters. The same is also true for the Calinski-Harabasz indexing but without any fixed range. However, in the case of the Davies-Bouldin index, a lower score indicates better clustering with the ideal score as zero.

For $2$-state $5$-neighborhood CAs under null boundary condition, there are a total of $226$ rules which are reversible for some $n\in\mathbb{N}$ [15]. Among these CAs, our first criterion is to filter out the CAs having a strict locality property. This indicates that, for these CAs, at the time of evolution, there is not much change in the states of the cells while hopping from one configuration to the next. This can happen when there is less amount of information propagation from one configuration to the next. Also, the rule needs to have a high rate of self-replication, that is, for the rule, there are many RMTs ($x_{i-2}x_{i-1}x_{i}x_{i+1}x_{i+2}$), for which the next state is the state of the cell itself, that is, $\mathcal{R}(x_{i-2},x_{i-1},x_{i},x_{i+1},x_{i+2})=x_{i}$. Here, $x_{i}$ denotes the state of the $i^{th}$ cell. When the configuration contains a cell with such a neighborhood combination, the rule keeps the state of the cell unchanged. Because of these, the configurations that are similar belong to the same cycle.

Now, this information propagation in the CA is calculated with respect to every neighbor. For that, we take each pair of arguments to the rule that differs by only one neighborhood position and see if the corresponding next states are different. According to Definition 2.1, with respect to the $k^{th}$ neighbor, $\mathcal{E}^{k}_{i}$, for each $i$, is such a set. If the next state values of the RMTs of $\mathcal{E}^{k}_{i}$ are different, then that neighborhood combination plays no role in updating the state of the cell under consideration. For example, if $\mathcal{R}(x_{i-2},0,x_{i},x_{i+1},x_{i+2})=\mathcal{R}(x_{i-2},1,x_{i},x_{i+1},x_{i+2})$ for all values of $x_{i-2},x_{i},x_{i+1},x_{i+2}$, that means, the next state of $i^{th}$ cell is independent of any changes in its left neighbor. The cumulative sum of all such pairs per neighbor can be used to derive the rate of information propagation for that neighbor. For instance, in the above example, the rate of information propagation from this neighbor is zero. The detailed derivation of the formula used for calculating it can be found in Ref. [11, 16]. In short, for each neighbor $k$, $0\leq k\leq 4$, information propagation to the $k^{th}$ neighbor due to change in the current cell is

where

and

The information propagation rate considering the cell itself as its neighbor, that is, for $k=2$, is the self-replication rate which is the information propagation for the cell itself. The maximum information propagation rate for any neighbor is 100% and the minimum is 0%. To choose the rules having strict locality property means, the rules need to have high rate of self-replication and a low rate of information propagation. Hence, our first filtering criterion is:
Criterion 1: Select the CA rules such that the rate of information propagation for each neighbor is $\leq 75\%$ but not all of them are equal to $0$ or $75\%$ and the rate of self-replication is $\geq 75\%$.

Table LABEL:tab:revRules lists the $162$ possible candidate rules out of the total 226 rules, which satisfy our criteria. Here, the first column enlists the rule in decimal value and the next five columns note down the rate of information propagation (in percentage) for that rule concerning each of the five neighbors. All of these rules are reversible for all $n\in\mathbb{N}$. But $162$ is still a huge number. So, we need to apply stricter criteria to these rules by analyzing the inherent cycle structure of the CAs to further reduce the candidate rule set size.

A good clustering implies similar data points are placed together. In the case of CA, as cycles are instrumental in deciding which data points will be in the same cluster, we need to ensure that, in the chosen CA, the configurations inside each cycle maintain a minimum distance. As we have chosen frequency-based encoding, the objects have been encoded maintaining minimum hamming distance between consecutive ranges. So, as per this encoding, the target configurations (encoded objects) which have minimum hamming distance between them, need to be placed in the same cluster to make the clustering effective. This inevitably indicates, that in our CA, the hamming distance between the configurations of each cycle needs to be as minimum as possible.

To maintain this, we calculate the average hamming distance ($\mathcal{H}_{avg}$) between the configurations of each cycle. As the configurations are in binary, this average hamming distance is just bit-wise XOR of the configurations. So, for an $n$-cell CA with number of cycles $\mathcal{K}$, where cycle number $i$ ($cycle_{i}$) contains $k_{i}$ elements $c^{i}_{1},c^{i}_{2},\cdots,c^{i}_{k_{i}}$, the average hamming distance ($\mathcal{H}^{i}_{avg}$) is calculated as:

We need our rules to have this $\mathcal{H}^{i}_{avg}$ to be as small as possible for all $i\in\mathcal{K}$. The minimum value is $0$. However, if we restrict the rules to have all $\mathcal{H}^{i}_{avg}=0$ for every $i\in\mathcal{K}$, then we are left with no rule in our selected $162$ rules. If we take the CAs to have all $\mathcal{H}^{i}_{avg}\leq 9$ for every $i\in\mathcal{K}$, then for $n=13$, we have only one rule $4042321935$ where out of the total $56$ cycles, average hamming distance is $0$ for each of the $48$ cycles and $1$ for each of the remaining $8$ cycles. So, restricting all $\mathcal{H}^{i}_{avg}\leq 9$ is not practical. We need to have a good number of cycles with this low average hamming distance such that there are enough rules in the list. This leads us to our second filtering criterion:
Criterion 2: For any $n$, choose those CA rules out of Table LABEL:tab:revRules for which there exists at least ${l}_{1}$ percent of cycles in the CA, such that in each of those cycles ($i$), the average hamming distance $\mathcal{H}^{i}_{avg}\leq 9$.

This $l_{1}$ can be set according to the user’s requirement. For example, considering $l_{1}=40\%$ for $n=13$, we get a list of $45$ CA rules. This list is shown in Table 5. Similarly, for other $n$, we can derive the list. For this work, the value of $l_{1}$ chosen for different $n$ along with the number of rules satisfying that criterion is shown in Table 6 Now, although these rules have a good number of cycles maintaining minimum hamming distance, many of them have a very large number of cycles. To achieve the desired number of clusters in less time, we must have some CAs with a limited number of cycles. So, to include such CAs, we relax our constraint on the hamming distance a little – we allow CAs with a limited number of cycles such that at least 50% of the cycles have an average hamming distance of only two digits. Therefore, our third and final selection criteria is:
Criterion 3: For any $n$, choose those CA rules out of Table LABEL:tab:revRules for which number of cycles in the CA is $\leq l_{2}$ and at least 50% of these cycles have the average hamming distance $\mathcal{H}_{avg}\leq 99$.

Here also, the maximum allowed number of cycles ($l_{2}$) can be chosen based on user requirements. Here, for $n=13$, we fix $l_{2}=40$ and get $4$ such CA rules following Criteria 3. These rules are also shown in Table 2. Table 6 records the value of $l_{2}$ chosen in our work for each $n$ and the number of rules we can get for that $l_{2}$. In this table, for $n=6,7$, a number of rules are marked with * because, as $n$ is very small, there are only a few configurations. For example, for $n=6$, the configurations are $000000(0)$ to $111111(63)$. So, the hamming distance between configurations will always be two digits. Hence, we take only CAs with a number of cycles = 2 such that all cycles maintain average hamming distance $\leq 9$. We get $6$ such rules for $n=6$: $252702735$, $1263225675$, $3789677025$ and $4042321935$ with an average hamming distance of the cycles as $0$ and $260960271$ and $756019215$ where an average hamming distance of both the cycles is $9$. Similarly, for $n=7$, restricting the CA to have only four cycles we get only $9$ rules $252695055,252702735,1263225675$, $3035673735$, $3785744805$, $3789677025$, $4041289185$, $4042310415$ and $4042321935$ which maintain average hamming distance $\leq 9$ for $50\%$ of the cycles.

Our final list of CA rules is the union of the set of rules selected following criteria 2 and 3. These CAs are good candidates for effective clustering having a good mix of rules with very low hamming distance and a limited number of cycles. The number of such rules is shown in the last row of Table 6. For example, for $n=13$, there are $48$ such unique rules (the rule $252691440$ is repeated in both criteria). The detailed list of selected rules for each $n=\{6,\cdots,12\}$ is shown in Table 7. With these rules, we can proceed to our clustering algorithm where we use only these rules for clustering with binary CA.

The dataset consists of $t$ number of objects, where each object is a $p$ bit binary string. Let $F_{D}=\{obj_{1},obj_{2},\cdots,obj_{j},\cdots,obj_{t}\}$ be the set of $t$ objects where $obj_{j}=b_{1}^{j}b_{2}^{j}b_{3}^{j}\cdots b_{i}^{j}\cdots b_{p}^{j}$. Here $b_{i}^{j}$ represents $i^{th}$ bit of $j^{th}$ object. In a low-dimensional dataset, the length $p$ of each object will be small. For such cases, we can take this $p$-bit string directly as the configuration of our CA. But in the case of the high-dimensional dataset, the length $p$ will be large. To reduce computational costs, we must split each object into many splits and then take each split into separate configurations. As shown in Figure 2, each object is split into $s$ configurations such that $obj_{j}=\langle split_{1}|split_{2}|\cdots|split_{s}\rangle$ where the string ‘$b_{1}^{j}b_{2}^{j}b_{3}^{j}\cdots b_{i}^{j}$’ $\in$ $split_{1}$, ‘$b_{i+1}^{j}b_{i+2}^{j}b_{i+3}^{j}\cdots b_{2i}^{j}$’ $\in split_{2}$, $\cdots$,‘$b_{si+1}^{j}b_{si+2}^{j}\cdots b_{p}^{j}$’ $\in$ $split_{s}$. The size of each split is equal to the cell length of the CA we choose, and the number of splits depends on the split size. That means $\lfloor s=\dfrac{p}{n_{1}}\rfloor$, where $s$ is the number of splits, $p$ is the length of the object in binary string format, and $n_{1}$ is the split size. Let $x\restriction_{split_{i}}$ represent the set of strings in the $i_{th}$ split of all objects. If there are $t$ objects, then $\lvert x\restriction_{split_{i}}\rvert$ is $t$.