Uniform tensor clustering by jointly exploring sample affinities of various orders

Throughout the paper, lowercase letters represent scalars, while bold lowercase letters represent vectors, such as $\bm{v}\in\mathbb{R}^{m}$. Bold uppercase letters denote matrices, such as $\bm{M}\in\mathbb{R}^{m\times n}$, while bold calligraphy letters represent tensors, such as, $\mathcal{T}\in\mathbb{R}^{m\times n\times l}$. Let $\mathcal{T}(:,:,i)$, $\mathcal{T}(:,i,:)$, $\mathcal{T}(i,:,:)$ represent the $i$-th frontal, lateral and horizontal slices of an order-3 tensor $\mathcal{T}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$, respectively. $\mathcal{T}(:,:,i)$ can be abbreviated as $\mathcal{T}^{(i)}$. A tensor can be reformatted into a matrix via series operations, called unfolding. For example, the three-order and four-order tensor unfold operations are as follows:

Three matrix multiplications are considered, namely, Hadamard product, Kronecker product and Khatri-Rao product [27].

A popular product between a tensor and a matrix is called a $k$-mode product.

For convenience, we summarize the frequently used notations and definitions in Table I.

Spectral clustering [3] starts by computing a sample-wise affinity matrix $\bm{S}_{2}\in\mathbb{R}^{m\times m}$. Every entry $\bm{S}_{2}(i,j)$ measures the proximity between two samples $x_{i}$ and $x_{j}$. Then, a normalized affinity matrix $\hat{\bm{L}}_{2}$ is estimated with $\hat{\bm{L}}_{2}=\bm{D}_{2}^{-\frac{1}{2}}\bm{S}_{2}\bm{D}_{2}^{-\frac{1}{2}}$, where $\bm{D}_{2}$ denotes the degree matrix of $\bm{S}_{2}$ with $\bm{D}_{2_{ii}}=\sum_{i}\bm{S}_{2{ij}}$. The Laplacian matrix $\bm{L}$ can be generated through $\bm{L}=\bm{I}-\hat{\bm{L}}_{2}$. Spectral graph clustering seeks a low-dimensional embedding of samples by maximizing the cross-entropy while preserving the local proximity of paired samples:

It is a standard eigenvalue problem and thus can be efficiently solved by calculating the eigenvectors $\bm{V}$ of the Laplacian matrix $\bm{L}$, using methods such as singular value decomposition. Typically, the normalized cut of the graph is minimized in this model, but as stated in [28], this model can also be expressed as a maximization problem with the normalized pairwise affinity $\hat{\bm{L}}_{2}$. By using the tensor $k$-mode product, the spectral clustering in Eq. (7) can be rewritten as:

where $\bm{v}_{j}$ is the $j^{th}$ column of the partition matrix $\bm{V}$. This equation establishes the relationship between the membership $\bm{v}_{j}$ and the pairwise affinity. Then, a clustering task is performed on the obtained embedding.

The pairwise affinity is easily broken by noise contamination or concentration effects in HDLSS data. To address this issue, our recent work IPS2 [26] attempted to use high-order affinity among more than two samples. IPS2 used the fourth-order tensor affinity to measure the proximity of two sample pairs. The decomposable fourth-order tensor affinity can be defined as the product of two pairwise similarities:

Let $\hat{\bm{T}}_{4}$ denote the matrix unfolded from a fourth-order decomposable tensor affinity $\mathcal{T}_{4}$; it has been shown in  [26] that

This result establishes a direct connection between a fourth-order decomposable tensor affinity and the usual second-order matrix affinity $\bm{S}_{2}$. Now, consider a normalized affinity matrix $\hat{\bm{L}}_{4}$, computed by $\hat{\bm{L}}_{4}=\hat{\bm{D}}_{4}^{-\frac{1}{2}}\hat{\bm{T}}_{4}\hat{\bm{D}}_{4}^{-\frac{1}{2}}$. Here, the diagonal matrix $\hat{\bm{D}}_{4}$ is the degree matrix of $\hat{\bm{T}}_{4}$ with the diagonal elements being $\hat{\bm{D}}_{4_{ii}}=\sum_{j}\hat{\bm{T}}_{4_{ij}}$. It has been shown that the fourth-order and second-order normalized affinity matrix $\hat{\bm{L}}_{4},\hat{\bm{L}}_{2}$ shares a similar relation as in Eq. (10):

One can reshape the eigenvector $\hat{\bm{v}}$ into a matrix $\bm{V}\in\mathbb{R}^{m\times m}$, and the matrix $\bm{V}$ is an eigenmatrix of the fourth-order tensor affinity $\hat{\bm{T}}_{4}$ in Eq. (10). The established equivalence between the decomposable tensor affinity $\mathcal{T}_{4}$ and any similarity matrix $\bm{S}$ opens a new path to understand high-order affinity. However, the fourth-order decomposable affinity inherits the drawback of the low-order similarity [29] in that it is vulnerable to noise corruption and concentration effects when applied to high-dimensional samples. Alternatively, one can use undecomposable tensor affinity to provide complementary sample affinities that the pairwise matrix affinity could not provide. The IPS2 uses a Fisher-ratio-like tensor affinity:

for $i,j,k,l\in n$, where $d_{ij}$ denotes the distance between samples $x_{i}$ and $x_{j}$ and the parameter $\sigma$ is an scaling constant.

One can calculate the normalized tensor affinity $\mathcal{L}_{4}$ from the fourth-order tensor affinity $\mathcal{\bm{T}}_{4}$. Then, we can learn a low-dimensional embedding $\bm{V}$ from the fourth-order tensor affinity by solving the following maximization problem:

We now formulate a uniform framework to learn a low-dimensional embedding from affinities of various orders. Suppose we have $n$ samples $\bm{X}=[\bm{x}_{1},\bm{x}_{2},...,\bm{x}_{m}]\in\mathbb{R}^{n\times m}$. The aim of clustering is to assign the samples into $k$ disjoint clusters, $C_{1},C_{2},\dots,C_{k}$. We introduce an entropy, called total normalized similarity, to quantify the sample assignment:

Let an indicator matrix $\bm{H}=[\bm{h}_{1},\bm{h}_{2},\dots,\bm{h}_{k}]\in\mathbb{R}^{n\times k}$ be the sample assignment such that $\bm{\bm{H}}_{ij}=|\bm{C}_{j}|^{-1}$ if $x_{i}\in\bm{C}_{j}$ and zero otherwise.

One can maximize the normalized associativity in Eq. (22) to seek an optimal sample assignment ${\bm{C}_{i},\dots,\bm{C}_{k}}$. Via algebraic transformation, this problem is equivalent to solving:

where $\bm{h}_{j}$ denotes the $j^{th}$ column of $\bm{H}$.

One can verify that solving the maximum normalized associativity is NP-hard. Alternatively, one can relax the binary assignment matrix to the orthonormal matrix $\bm{V}\in\mathbb{R}^{n\times k}:\bm{V}^{T}\bm{V}=\bm{I}$ to lower the expectation of binary assignment. The simplified problem then becomes solving

where $\bm{v}_{j}$ represents the $j^{th}$ column of $\bm{V}$.

We now fuse the aforementioned undecomposable third-order and fourth-order tensor affinities with the pairwise affinity matrix to seek a uniform low-dimensional embedding. The high-order tensor affinity can be used to effectively address the challenge of the vulnerability to high-dimensional datasets. The proposed model, named uniform tensor clustering (UTC), effectively fuses the similarity measurements among samples to yield a uniform low-dimensional embedding to approximate the partition matrix. It achieves the goal by maximizing the cluster’s associativity of various orders. The proposed model is formulated as:

where $\bm{v}_{j}$ is the $j^{th}$ column of the matrix $\bm{V}$. This unified model aims to fuse various affinity matrices to yield a uniform low-dimensional embedding that is robust to noise contamination and high-dimensional concentration effects. The proposed model enhances the clustering performance by exploiting spatial information from multiple (more than two) samples. Upon solving this problem, one obtains the embedding $\bm{V}$. The clustering task is accomplished by fusing low- to high-order information and using K-means clustering based on the embedding matrix $\bm{V}$ to obtain the final group. Therefore, once we have the indecomposable tensor affinity $\mathcal{S}$ from different orders, computing the fused low-dimensional embedding yields the different orders’ spatial information to enhance the clustering performance.

The optimization problem in Eq. (25) is convex, and the objective is differentiable. However, the gradient of the objective function involves solving a polynomial function of order three. It has no explicit solution and thus is computationally intractable. To circumvent this issue, we introduce a slack variable $\bm{v}_{2}$ to approximate the term $\bm{v}\otimes\bm{v}$, thereby avoiding having to solve a high-order polynomial function. The model can be rewritten as

This problem can be solved by alternating direction minimization. Its augmented Lagrangian formulation is as follows:

where $\bm{Y}_{1}$ and $\bm{Y}_{2}$ are Lagrange multipliers. The constant $\mu>0$ is a penalty parameter. We then decompose the problem into two subproblems with respect to the variables $\bm{V}_{1}$ and $\bm{V}_{2}$ alternatively. Each variable is solved while fixing other variables. The process is updated iteratively until convergence.

Step 1): Solving the subproblem with respect to the variable $\bm{V}_{1}$

By fixing the variable $\bm{V}_{2}$, the problem can be simplified as:

The gradient of the quadratic function is:

Thus, the problem in Eq. (28) can be solved via gradient descent or by setting the gradient function Eq. (29) to be zero to obtain the solution directly.

Step 2): Solving the subproblem with respect to the variable $\bm{V}_{2}$

By fixing the variable $\bm{V}_{1}$, the augmented Lagrange function can be simplified as:

The gradient of the objective function is:

By setting the gradient to zero, one can obtain the implicit solution as:

Step 3): Updating the multipliers $\bm{Y}_{1}$ and $\bm{Y}_{2}$:

The three steps are iteratively updated until convergence or until a stopping criterion is met: $max(\|\bm{v}_{1}^{k+1}-\bm{v}_{1}^{k}\|_{\infty},\|\bm{v}_{2}^{k+1}-\bm{v}_{2}^{k}\|_{\infty},\|\bm{v}_{1}^{k+1}\otimes\bm{v}_{1}^{k+1}-\bm{v}_{2}^{k+1}\|_{\infty})<\epsilon$. Since the main problem is convex, the numerical scheme is guaranteed to be convergent. The alternate strategy to solve the problem is summarized in Algorithm 1.

We employ five popular clustering methods to benchmark the proposed UTC. A brief introduction to these methods is given below:

• •

  Spectral clustering (SC) [3]: The classic spectral clustering gives a baseline on behalf of pairwise similarity.

• •

  Clustering with Adaptive Neighbors (CAN)  [17]: Dynamically learning the affinity matrix by assigning the adaptive neighbors for each data point based on local distance. The final similarity matrix’s connected components are made precisely equal to the cluster number by imposing a new rank constraint.

• •

  Clique averaging+ncut (CAVE)  [21]: The hypergraph, on behalf of the high-order affinity, is approximated using clique averaging, and the resulting graph is partitioned using the normalized cuts algorithm.

• •

  Pair-to-pair clustering (PPC) [26]: Tetradic undecomposable affinities are used to obtain the low-dimension embedding as the final similarity matrix of the spectral method.

• •

  Integrating tensor similarity and pairwise similarity (IPS2) [26]: Applying the spectral clustering method on the fused similarity that combines the pair-to-pair affinities and pairwise similarity.

• •

  Uniform Tensor Clustering (UTC): The proposed method.

In terms of the type of the affinity, the five methods can be divided into three categories, i.e., pairwise affinity, high-order affinity, and fused-order affinity.

The performance of each method is evaluated with respect to five popular metrics: accuracy (ACC), adjusted Rand index (ARI), F-score, normalized mutual information (NMI), and purity (PURITY). The larger the value is, the better the performance is.

Syndata-1 consists of 40 samples from two different subgroups. Each group is perfectly distributed along a straight line, while the two groups are distributed perpendicularly to each other. We use circles and crosses to denote the samples from the two subgroups, and their distribution is shown in Fig.  3(a).

We applied the spectral clustering, which works on pairwise affinity, on Syndata-1. Since the samples from the two subgroups are perpendicular, the resulting affinity matrix is non-differentiable. Therefore, the SC fails to yield the correct assignment. It mistakenly classifies the samples into two parts, with most in the upper and the rest in the bottom. In comparison, UTC can identify the sample assignments perfectly with $100\%$ accuracy. The results indicate that the proposed undecomposable high-order affinity can explore complementary information to improve the clustering performance.

Syndata-2 consists of three subgroups, each derived from independent and identically distributed normal distributions with an equal standard deviation of 0.5 and mean of 2. The data from different groups are orthogonal, and each group contains 30-40 samples.

The pairwise affinity tends to be ambiguous with expanding feature dimensionality, limiting the clustering performance. We applied the proposed UTC on Syndata-2 to validate its stability in handling high-dimensional data by comparison with popular methods, including SC, CAN and CAVE, PPS, and IPS2. The results are shown in Fig. 4. When the feature dimension varies from 10 to 10000, UTC can maintain robustness to perform accurate clustering. Furthermore, the baseline methods gradually lose the ability to identify the data structure.

Moreover, when the number of data dimensions reaches 10000, Tab. II shows that our proposed UTC methods, including fusing pairwise with triadic affinity, tetradic affinity, and both, substantially outperform all other state of the art methods. The ACC clustering results of fusing triadic and fusing tetradic are higher than those of the second-highest baseline method by 14% and 28%, demonstrating the effectiveness of fusing high-order affinity with pairwise relation under HDLSS data. Additionally, the ACC achieves 100% when utilizing both higher-order affinities, proving that fusing different-order affinities is more helpful in completely explaining the data structure than is merging only a single higher-order affinity.

We also employ the t-SNE to visualize the differences in the embedded features after different methods. Most of the samples mix together when visualizing their embedding by t-SNE on the raw data, as shown in Fig. 5(a). On the other hand, the pairwise affinity, shown in Fig. 5(b), has no clear-cut structure, except the diagonal elements having large values. Therefore, it is challenging to discriminate the three subgroups when using the pairwise affinity matrix directly. Now, we fuse the three-order tensor affinity with pairwise affinity and apply UTC to obtain a local embedding. t-SNE is applied to the resulting embedding to visualize the sample assignment distribution. The visualized result is shown in Fig. 5(c). The samples in one subgroup (Subgroup A) are clearly separated from the others, while the remaining two subgroups (Subgroups B and C) are not well-distinguished. The upper left of the affinity matrix shows an obvious blocky structure, but the remaining part lacks a clear structure, as shown in Fig. 5(d), which confirms the results obtained by t-SNE. Alternatively, if we fuse the fourth-order tensor affinity with the pairwise matrix affinity and apply UTC to obtain another local embedding, the visualization using t-SNE is shown in Fig. 5(e). The samples in Subgroup C clearly are separated from the other two subgroups, while the other two subgroups are merged. The affinity heatmap also validates the observation that Subgroup A in the lower-right corner possesses distinct value compared to the others, as shown in Fig. 5(f). Finally, we fuse the pairwise affinity with the third-order and fourth-order. UTC is employed to learn a uniform embedding that is then visualized via t-SNE, as shown in Fig. 5(g). The three subgroups are perfectly separated in this case. The diagonal elements of the corresponding affinity matrix also show three obvious block-like structures and an affinity difference between the three subgroups. Therefore, an accurate sample assignment can easily be obtained from the embedding by the proposed UTC. In summary, the above experimental results validate that the high-order tensor affinities make up for the insufficiency of the pairwise affinity matrix when corrupted by high-dimensional features.

We then validated the power of UTC by applying it to six public benchmark datasets. The six datasets are chosen to be representative of various sources, including facial image, object image, and bioinformatics data. The statistics for the six datasets are summarized in Tab. III. Additionally, a detailed description is provided.

• •

  Yale dataset has 165 grayscale images covering 15 different individuals. Each individual has 11 different facial and configuration images, such as with glasses, without glasses, center-light, left-light, happy, sad, and so on. We extracted 4096-dimensional raw pixel values of every image in 5 classes for our experiment.

• •

  Coil20 dataset contains 1440 images of 20 categories of objects. Each category has 72 images from different views. We sample 20 samples from 5 classes, and all images have a size of $128\times 128$.

• •

  Lymphoma dataset, one of the most common subtypes of non-Hodgkin’s lymphoma has two molecularly different forms of diffuse large B-cell lymphoma (DLBCL), which have gene expression patterns that indicate different stages of B cell differentiation. The dataset contains a total of 62 samples, 4702 based on expression fragments.

• •

  DriveFace dataset contains images sequences of subjects while driving in real scenarios. It is composed of 606 samples of $6400$ features each, acquired over different days from 3 drivers with various facial features, such as glasses and beards. For each type, we pick 26 samples for the experiment.

• •

  Lung dataset contains a total of 203 samples that can be divided into 5 classes, with 149, 21, 20, 6, and 17 samples, respectively. Each sample has 12,600 genes. We selected 40 samples of the first category and all samples of the remaining four categories for experiments.

• •

  GLI-85 dataset consists of the transcriptional data of 85 diffuse infiltrating gliomas from 74 patients. Those gliomas can be divided into two kinds of tumor subclasses. Each instance has 22283 features.

We applied the proposed UTC as well as five comparative methods on the six datasets, and the results are summarized in Tab. IV. The bold value represents the best result. The clustering method that uses higher-order affinity uniformly outperforms the standard method using pairwise affinity. Therefore, the higher-order affinity constructed from the relationships of multiple samples can better describe the spatial structure of the data, which is advantageous when dealing with HDLSS data. Furthermore, our proposed UTC method significantly outperforms the competitive methods with respect to all five evaluation metrics for all types of data, including facial image, object image, and bioinformatics data.

On the Yale dataset, the ACC clustering result of our method is over 33% higher than that of SC and over 8% higher than the second-highest IPS2. On the GLI-85 dataset, our results are more than 16% and 6% higher than those of the baseline and second-highest IPS2 methods. There are also 1%, 7%, 10%, 2%, and 3% improvements compared with the second-best performance on the Coil, DriveFace, Colon, Lymphoma, and Lung datasets, respectively. IPS2 and UTC achieve the top two results on almost all datasets with respect to nearly all the metrics.

We also visualized the sample spatial distribution after t-SNE[31] and the heatmap of the affinity matrix on the Yale, Coil, and DriveFace datasets in Fig. 6 and the Lymphoma, Lung, and GLI-85 datasets in Fig. 7, respectively. The spatial distributions after t-SNE on the raw data and the resulting embedding are shown on the left side of the two columns, while the heatmaps of the raw data and the affinity matrix fusing with various orders are on the other side. In most figures, the low-dimensional embedding after UTC displays better separation and less overlap, thus making it easy for clustering. In comparison, the heatmap of similarity produced by SC and UTC further demonstrates the effectiveness of our method. Deep yellow represents a large proximity value. The affinity heatmap on the raw samples possesses highly blurred boundaries, and there is no apparent block structure. However, the affinity matrix calculated from the low-dimensional embedding after UTC has clear-cut boundaries. The large value in the matrix generated by UTC is concentrated on the diagonal blocks, implying that the fused affinity can more accurately reflect the real relationship among the samples, thus yielding superior clustering performance in Tab. IV. Taking the Yale dataset as an example, the ACC of the clustering results obtained by UTC is 33% higher than that of SC. The visualization results confirm our method’s superiority intuitively, as shown in Fig. 6(a). Most of the samples are mixed together in the visualization based on pairwise affinity, resulting in an inability to distinguish between different subgroups accurately, as shown on the far left of Fig. 6(a). The heatmap of pairwise affinity confirms that only one cluster is significantly different from the other samples. Correspondingly, in the t-SNE visualization of the resulting embedding, most samples are clustered with their same cluster’s partners. Five obvious blocks are displayed on the diagonal in the heatmap. The right side of Fig. 6(a) confirms that UTC can better reflect the real structure of the data compared with single-order methods, which is critical for clustering.

In summary, the reason for the superior performance of UTC is that it utilizes the complementarity of high- and low-order affinities to comprehensively capture the spatial structure of data, thereby leading to a decisive effect on the clustering performance.

The UTC consists of two major computational components: constructing the high-order tensor affinity and solving the minimization problem in Eq. (27). In the first component, we calculate both the triadic and tetradic affinities by means of their $k$-nearest neighbors (KNN). Given $m$ samples, one can use the KNN strategy to reduce the time complexity from $\mathcal{O}(m^{4})$ to $\mathcal{O}(mk^{4})$. Thus, the two affinity tensors are sparse, with the number of nonzero elements being $s$. For the latter, we solve th model by means of an iterative strategy. The subproblem of $\bm{V}_{1}$ based on the gradient descent strategy has a computational complexity of $\mathcal{O}(ncsm^{2})$, where the $c$ denotes the number of clusters and the maximum iteration number is $n$. Solving the other sub-problem $\bm{V}_{2}$ takes $\mathcal{O}(csm^{2})$ using the Germs algorithm. Thus, each iteration has a total computational cost of $\mathcal{O}(ncsm^{2})+\mathcal{O}(csm^{2})=\mathcal{O}((n+1)ncsm^{2})$. In total, the complexity to learn the uniform embedding in Algorithm 1 is $\mathcal{O}(mk^{4}+K((n+1)ncsm^{2}))$, with $K$ being the number of iterations.