Fused Orthogonal Alternating Least Squares for Tensor Clustering
Authors: Jiacheng Wang
Dan L Nicolae
Venue: NeurIPS 2022 Accept (2022)
Keywords: High-order tensors
Clustering
Tensor decomposition

Abstract: We introduce a multi-modes tensor clustering method that implements a fused version of the alternating least squares algorithm (Fused-Orth-ALS) for simultaneous tensor factorization and clustering.  The statistical convergence rates of recovery and clustering are established when the data are a noise contaminated tensor with a latent low rank CP decomposition structure. Furthermore, we show that a modified alternating least squares algorithm can provably recover the true latent low rank factorization structure when the data form an asymmetric tensor with perturbation. Clustering consistency is also established. Finally, we illustrate the accuracy and computational efficient implementation of the Fused-Orth-ALS algorithm by using both simulations and real datasets.

Fused Orthogonal Alternating Least Squares for Tensor Clustering
Jiacheng Wang
Dan Nicolae
Department of Statistics University of Chicago jiachengwang @ uchicago.edu
Department of Statistics University of Chicago
nicolae @ statistics www.stics.uciccoo.ddu
Abstract
We introduce a multi-mode tensor clustering method that implements a fused ver- sion of the alternating least squares algorithm Prthe--Ott--ALS) for simultaneous tensor factorization and clustering. The statistical convergence rates of recovery and clustering are established when the data are noise contaminated tensor with a latent low rank CP decomposition structure. Specifically, we show that a modified alternating least squares algorithm can provably recover the true latent low rank factorization structure when the data form an asymmetric tensor with perturbation. Clustering consistency also established. Finally, we illustrate the accuracy and computational efficient implementation of the Fused-Orth-ALs algorithm by using both simulations and real datasets.
1 Introduction
Tensors, or multidimensional arrays, have played an essential role in a wide range of scientific and business applications, including neuroscience (9 [66]. social networks [18] 4 and recommendation systems 3 28). As higher order tensors have been extensively used as a framework for storing and organizing massive data, the associated need to develop methods for tensor multi-mode clustering has sharply increased. Several methods have been proposed for matrix biclustering for columns and rows simultancously 7 22 However, those methods cannot be directly applied to tensors due the complex higher order generalizations of the matrix singular value decomposition (SVD) and principal component analysis (PCA). Motivated by the methodology gap between low and high dimensional arrays, we investigate here the multi mode clustering problem for tensors.
We propose the Fused Orthogonal Alternating Least Squares (Fused-Orthe -rtted--Ot--L)) algorithm for uncov- ering co-clustering structures in high order tensors. Clustering along specific mode is equivalent to clustering the slice matrices for that mode as shown in Figure The core idea is to model true cluster means as having rank-l matrix basis. The method is designed to systematically investigate and achieve the following. First, under mild assumptions on the tensor decomposition structure, our approach encourages smoothness in the decomposed components, leading to clustering performance with higher accuracy. Second, orthogonality is pursued for the decomposed components, avoiding the problem of local-minima convergence of the classical ALS tensor decomposition algorithm, preventing multiple recoveries of the same factors, and achieving faster convergence rates compared to popular tensor decomposition methods. Third, theoretical guarantees are provided for the recovery and clustering consistency with even single tensor sample, which is difficult achieve in vector matrix clustering. The results in synthetic and real data demonstrate the robustness of our approach under mild model misspecification and assumption violations, and the effectiveness under high dimensional settings for extracting underlying clustering structures. Even though we focus on order three tensors, is straightforward to derive analogous results for higher order tensors.
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

Related work. Multiple papers are closely related to, but also clearly distinct from, our proposed methodology for tensor muulti-mode clustering. Tensor clustering utilizing tensor block model appears in Chi et al. g and is further developed in Wang and Zeng (26] A popular alternative to the tensor block model the CANDECOMPPAARAAAG (CP) decomposition that captures structure using the sum of rank-1 tensors. Our paper adopts the CP decomposition because it handles heterogeneity in each mode, learns the clustering patterns across different modes of data in a more independent way, and provides flexibility for clustering a certain mode of the tensor without being affected by correlation with other modes. Our method similar to those in a recent series of papers [27] 21 that use the CP decomposition structure. Note that their estimation algorithms use the framework of tensor power method [. In contrast, our algorithm utilizes Alternating Least Square (ALS), a widely employed decomposition algorithm [4. The ALS algorithm estimates the CP decomposition by specifying factor matrices simultaneouslyy, instead of the column-by-column recovery for each factor matrices implemented in tensor power methods. The ALS algorithm has been proved to be robust and computational efficient (1O, 12 [5]. In particular, compared to Wang et al. [27] which performs clustering over estimated factor matrices, our method encourages smoothness in the factor matrices by imposing generalized LASSO penalty, resulting in apparent clustering structures in factor matrices since data from same cluster will possess the same value by choosing appropriate tuning parameters. In addition, one critical assumption in Sun and Li [21 is that samples from the same cluster are ordered consecutively, an assumption that is rarely valid in applications. Our approach employs a generalized LASSO regularization that is based on graphical distance intuition learned from initial observations. We demonstrate that the generalized LASSO penalty achieves better clustering performance when clustering assignments are randomly distributed. Even though Fused- Orth-AL.S is motivated by the Orth-ALS algorithm [9), we make several advances. First, Orth-ALS only considers tensor decomposition performance under the noiseless situation, and we successfully provide statistical guarantees for Fused- Orth-ALS under any error model, with a mild constraint on error tensor spectral norm. Second, theoretical properties for Orth-ALS are established only for symmetric tensors, while Fused-Orth-All ALS can achieve both recovery and clustering consistency for any asymmetric tensors.
Notations. Calligraphic letters such as ) Rd xd used to denote order three (d1.d2, da)- dimensional tensors. By convention, we use a colon to indicate all elements of a mode. Thus, mode-1, 1. mode-2, mode-3 slices of y can be represented as V.:., V.:: V.kk Bold uppercase and lowercase letters are used to represent matrices and vectors respectively, such as A E Rdixds xd2 and u E Rd = [u,,..., Ud,]T. The ith row and jth column of A are denoted as A; and A.y. For notational convenience, we define the multilinear combination of tensor with three vectors U1, U2 us as Y(u, U2, Ug)= j.k 21, U3 In particular, the spectral norm of an order three tensor y is defined by | = max [b=1|ualle -|=-||==l Σ i.j.A U1 U3 & Vitk. The Frobenious norm of
J' is the natural generalization of the Frobenius norm of a matrix, || F V..j.. Diyk- Lastly, [K] = (1, 2.., K} stands for the whole index set and we o, O represent outer product and Khatri-Rao product. For two positive sequences }. }. On S bn means on I Cb, for some constant C¹ 0 independent of
O.. - 8 Mode-1 siices
Modee
u r6 Mode-3
CoG. Mode-3 slices
Mode-2 s
Mode-2 slices
Figure 1: Fused-Orth-ALS algorithm performs multi-modes clustering through detecting cluster structures of slices along each mode.

---

2 Clustering model via tensor decomposition with regularization
Given observed tensor y x Xda that noisy version of a tensor of interest ) i.c. y ) + E, our goal is to uncover the clustering structures for the three modes of 3 We assume that y" tensor with a rank K CP decomposition structure y ΣlE(K] W; o B C.,.
X B R C E RRdx* are factor matrices with columnwise unit norm, i.e. | = |B,, =I -I,, = 1, E [K], and w = [..,..., wk] TRA captures the weights of the factor matrices. The construction of ) leads to dependence of clustering structures on the factor matrices along each mode. Specifically, suppose there are 8 clusters for the first mode, and each row in A can be expressed as A.: Σド= Mi., Ligay, where Mi.j (I.1.)) "(1.:.....(............. ((I.)). k]] is the mean value for rows in jth cluster and 2 contains all the rows indices that belong to jth cluster over first mode. Consequently, the true cluster means A1, As, , for slices along first mode (Di::..., Vd.:: as indicated in Figure [ can be written as
where A E R
K
K A. Σ (L.) B4 C.
A
Σ
W5)(L,1), B, Ca, 1=1
This reveals the core idea of clustering via CP decomposition: all the true cluster means share the same K rank-1 1 matrix basis B; C.,, Vi E [K]. Thus, for every i, we would expect there exists s E [s1] such that E 2 and as a result, the clustering label for ith element along the first mode will be As shown in Sun and Li (21), the tensor Gaussian mixture model can be viewed
special case of this clustering structure. However, finding the partitions 2 Vi E [81] is a combinatorial hard problem and could be computational intractable. Inspired by regularization term proposed in Chi et al. ( we impose full pairwise difference operator over rows of A to yield a weighted penalty on local differences,
|AAAIi
Σ (E,E
Mh. |A:
A
where S : ( i2)l in in [di), E [d2]] is containing pairs of different row indices and parameters 7d = (7d, (11, 12) E S) are non-negative weights that are critical for controlling the penalty imposed on pairwise row differences for A; we will explain detail how to choose them later. Analogously, μ2,j Bj. AA, y² and μ3.j. Cj. 'A, 7 can be defined to characterize the clustering structure and regularization along the second and third modes.
In summary, we propose the following penalized constrained optimization as an approach to reveal the latent clustering structure,
min || Σ A
wi A B. C + X11 AAAI + AAB | +| AAC|]]]
A.B.C.w
s.t. | = B. =
[K]
(1)
In particular, increases, A; will shrink towards each other, indicating pairwise differences of rows in A will become increasingly sparse. Sparsity in pairwise row differences naturally leads partitioning assignment 2 Vj E [s1]. For simplicity, we choose be the same over three modes. We would like to draw the attention to the difference between the regularization term in (I) and that in Chi al. g. We use é, norm penalty which will yield pairwise difference equal to due to feature selection for & regularization, while Chi et al. [ adopts the Frobenious norm which gradually. As mentioned, dynamic tensor clustering (21) used a fused LASSO penalty. However, this simple fusion structure only works if samples from same cluster have consecutive indices. Obviously, our weighted pairwise difference operator 'A is a generalization, and later we will show how our method outperforms dynamic tensor clustering when the cluster labels are randomly assigned.
only shrinks the pairwise difference
Implementation: fused orthogonal alternating least squares algorithm
3.1 Algorithm
Starting from the alternating least square (ALS) algorithm 14] we propose the fused orthogonal alternating least squares (rused-OOth-ALS)) method described in Algorithm by adding two additional
3

---

steps that are designed for each iteration of ALS, which allows for the avoidance of local optima and more rapid convergence to the true factors. Moreover, a 'Fuse" operator with generalized LASSO regularization A is employed on each column of AL. estimates which in general is defined as
clustering. First, the orthogonalization step is performed before
1
Fuse u1 A, A) arg min vER 2
Σ
(0y
uz)²
A|
j=1
Algorithm 1: Fused-Orth-ALS Algorithm
Input: Tensor tensor CP rank
Output: D = Σx= ú; B', C
Initialize A BB, B
Co with columns randomly from the unit sphere and iteration index t = 0:
while termination condition not satisfied t +++1 1; 1;
Find QR decomposition of 1-1 set At-1 Q. Orthogonalize - and Čr-1 similarly;
X + (C)) ④ B 1---1 ) X : =[ with columns X <- Fuse X A A);
Y (/22((C-11 1+ Y = Y with columns Y f- Fuse Y:: ?A, A);
Z' - (aaa((B @ A' Z' = 12 with columns Z <- Fuse Z A A);
Normalize X Y Zt with columnwise norm unit and store the results Á' B Ö
end
Estimate weights w; = VA 187 B CD), V E [K]
We use the 'Fuse' operator since the regularization term in (I) can be viewed as the sum of 'Fuse" regularizations on each column of factor matrices, e.g., AAllI Σ AA and it be solved efficiently via existing methods [2] 29) After obtaining estimates A, B, C, clus-
tering algorithm such as k-means or hierarchical clustering are used to obtain cluster assignment
2, B,, Ce Vi E [81]..] E [s2], k E [+s]-
From computational complexity perspective, orthogonalization, ALS updates, and Fuse operations in Algorithm I take O (d + d2 d3)). O(Kd1d2d3) and O (dk + + +d)) number of operations respectively. Thus, in a general case, the total computational complexity of Fused-Orth-ALS algo- rithm on (d1, d2, ... d dimeensionll tensor O(K (max (K Σ dy, DI[; d, Σド= d)))). When dj are of the same order d and tensor rank K = o(Da the total complexity can be simplified to O(K DdP) which is the same order as classical ALS algorithm, meaning that extra orthogonalization and "Fuse" steps do not increase the computational complexity. The proposed clustering model is clearly distinct from dynamic tensor clustering (DTC) [21] in terms of implemen- tation of generalized LASSO regularization, ALS updates and added orthogonalization. Our method has the same computational complexity as DTC but orthogonality leads to significant speedups in the iterations required for recovery convergence, which will be discussed in detail in section 5 later.
3.2 Regularization weights and tuning parameters
In practice, the choice of appropriate regularization weights y ,²,72,7" is critical since it affects clustering accuracy and computational efficiency. For instance, Thy.i characterize the shrinkage of row difference between A; and Aj:. From a graph perspective, each row of A can be treated as a node an undirected graph and nodes i are connected by an edge with weight 7h. All the edges compose the edge set S. Clearly, large value of 7hyd2 implies high similarity between A;; and Ag:: resulting in clustering nodes i1, 12 into the same group. In particular, our adopted A reduces to fused LASSO penalty utilized in DTC [21] in the case of chain graph where S = [i (i, i 1)] i 1,,, di 1}. A well-developed strategy [20 5 choosing regularization weights: we perform rank-K CP decomposition approximation to V, utilize the estimated factor matrix A to quantify the 'distance' among different nodes, and then calculate the regularization
espoused for
regularizattion weights through
13.r2 i exp (
(xxp-~11 [-v][A Aix.113)
(2)

---

Here, c,,.ig is an indicator function specifying whether Aig: is among k-nearest neighbors of A that controls density of graph. The smallest that ensures the graph connected is the default choice. The remaining component in (2) is Gaussian kernel with as measure of scale chosen as median Fuclidean distance between the i1th and i2th rows that are n-eerestt neighbors of each other. Lastly, we normalize (2) Viyii d2 + +g to ensure the penalty over three modes are on the same scale. Similar formulas are used for 72 and y².
We choose through the use of the extended Bayesian Information Criterion
21
26
by minimizing
log (1V 1 11
11
d
log d II j=1
d x
df,
J=1
where dfx is the number non-zero elements in estimated factor matrices, characterizing degrees of freedom. Moreover, rank K and clusters 81, 82, 83 are usually unknown and need to be estimated from data. We use the "elbow point' method to choose rank K by plotting the recovery error, and determine 81, by computing gap statistics 23
Cluster recovery and convergence
This section provides convergence analysis for the Fused-Orth-AIS algorithm. For simplicity, set di = d2 = d3 = d and choose equal regularization weights over all modes. Without loss of generality, assume Wmax w1 2 W2.... ... wk Wmin 0.
The technical assumptions needed for our analysis are as follows:
KCC, C. I a/ a/vd for some A1: Define incoherence = m (A A I (B.,
polylog(d) and K = o33 /polylog(d))/3) The spectral norm of y" satisfies ||>*|| S Wmax
A2: Define the initialization error as EQ = max; A A. 221 [ B [12 which satisfies
<
min
Winin
(k -5--1, 1);
Winin
2p(K 1) (K 1)( {
1- (K 1)(e 1)(++
12
72
1 (K 1)E(1 +
(K (K-11) 1)(1
with is constant upper bounded by 10K 10Kumaxx (wmin vd).
A3: Denote the spectral norm of & by vb such that 1b S min (Wmin/ 6, Winax K/d
A4: Denote M maxx
A |l2 and define & max; || C.lll that satisfies
Wmax a + 2 - 1) + 2 - 1)) + wt (2 (2Muniil( Wmin -
Assumption Al relaxes the requirements the orthogonality of columns in factor matrices and imposes weaker condition on the rank than Orth-ALS algorithm [9)] Assumption A2 specifies initialization criteria that are needed only for computational issues. The bounded perturbation is specified in A3. In particular, it can be proved with high probability if elements
are i.i.d from a sub-Gaussiar distribution stated in Tomioka and Suzuki (44. The bounded fusion assumption, A4, restricts clustering complexity on factor matrices. To illustrate its meaning, consider a simple scenario when all clusters along the third mode have the same size d/s3. For C. belonging the same cluster Cj. they should take the same value μ3.j under our model assumption. Consequently, A4 can be rephrased as 9 51 1/sa]); ); in the special case when C has only one cluster, e.g. $3 = 1, A4 reduces to AAC||1 = O. Thus, bounding is equivalent to bounding cluster size $3 when sample size over third mode is kept fixed.
We can now derive recovery error bound for the Fused-Oris Orth-ALS algorithm.
Theorem Assuming assumptions A1-A4 hold, factor matrix estimate C.,, Vi E [K] of Algorithm 1 satisfies the following error bound with high probability
||C C. "maxxp p²
+ up/ W'in
Wmin
when choosing appropriate X. Similar error bounds hold for the other two factor matrices estimates
A. B.
5

---

The error bound derived in Theorem I reveals how the weight w;, incoherence parameter p. rank K and perturbation level e interact with each other on affecting the convergence behavior of Fused- Orth-ALS algorithm. Lower perturbation level 4/, lower incoherence parameter rank K and lower signal ratio of Wmax/ Wmiin all result lower error bounds. Clearly, the error bound for each column in factor matrices has two parts: one is related to perturbation level qb, the other is dependent on the underlying CP decomposition structure, signal ratio Wmax/Wmin Wmin rank K and dimension along each mode. Thus, under high dimensional settings when d + ce, i.e., p² S a² /d - 0, the second part will dominate the error bound while if tensor data is almost noiseless, e.g. & Rs O, the first part will be the main source of error. The following corollary analyzes the relationship between these two parts under a special case when the enror tensor follows a uub-Gaussian distribution.
Corollary 1 Assume assumptions AI-A4 hold. If we further assume each element in error tensor,
Eijk. is independent, with zero mean and satisfies Ele S c
e Wmax/ Winin S C where C is
positive constant and the minimal weight satisfies
21 d (P/(K- - 1)2
Wmin
oo 3d
+ log 8
log 3/2
then C.,, Vi E from Algorithm 1 satisfies the following error bound,
C. C.ill2 2 1K11
with probability at least 1 6. Same error bounds hold for the other two two factor estimates A, B.
Thus, by imposing the sub-Gaussian distribution assumption on error tensor elements, recovery consistency for factor matrices estimated from Fused-Orth-ALS algorithm can be established.
Next, we demonstrate cluster consistency for the clustering algorithm performed on factor matrices recovered by Fused-Orth-Ar Saeed-Orth-ALSS algorithm. Since clustering algorithm is performed on rows of A, B, C,
clustering error is quantified through the true mean value and its estimate || f llz M2.ja
ju2.j 12, |A3.a P3,jall 2|2.
Theorem 2 Assume assumptions in Corollary 7 hold, we have
max
MI 11 S K1.5 /d
)1
max
[2.j 12.jal12 S K 1.5 /d
82
max |aa.s Ma.islla E K /d
hold with probability least - 8. Furthermore, if minj c(
m/m mym' ||
Mi lla 2
K1.5/ /d,
1, || - M2 ll2 2 K k1.5/d, min der Il
2 l2 Z K 1.5 /d, we have 2m = 21",, Bn B E = ej hold with probability at least
1 0.
This theorem shows that clustering consistency holds as long as K k155/d + 0, which allows rank K increase with the dimension d. At first look, the conclusion in Theorem 2 seems to indicate that the cluster mean error bound does not depend on number of clusters over each mode, $1, $2, 83. However, stated before, assumption A4 employs "invisible" bound on the number of clusters which is closely related with the convergence rate derived in Theorem We provide a detailed explanation in supplementary material stating cluster mean error bound k1. /d increases with the number of clusters 8;. Note that Theorem 2 assumes that true rank is known; we leave the effect of an estimated rank K on clustering consistency for further study.
Numerical experiments
5.1 Synthetic datasets
In this section, we investigate the performance of the Fused-Orth-ALS algorithm on small samples and compare the recovery and clustering error with alternative tensor-based clustering methods.

---

Recovery enror is evaluated through || -1 Clustering error measures the frequency of mismatches between the estimated clustering assignment Di and the true M over samples 21 ....... Zn. which is defined as
(G..)
+
< J2 E [m))/// nn(n - 1)/2
and can be easily computed after obtaining the cluster assignment 2m B C, Vm E [si], E [sg]. 1 [sy]. Note that the clustering error is similar to two other commonly used metrics: adjusted Rand index [I]] and variation of information [7]] For the simulations below, elements in the error tensor & are generated independently from a Gaussian distribution with mean and variance o², We compare Fused-Orth-it
Algorithm with three methods: dynamic tensor clustering (DTC) [21]] CP-Kmeans algorithm which performs a rank- K CP decomposition on the tensor observations first and then independently applies k-means algorithm clustering to the rows of derived factor matrices, and multiway clustering for tensor block models (26.)
5.1.1 Single-mode clustering
We first generate an order three tensor with CP-rank 2 for which we would like to do single-mode clustering along the third mode. The dimension of the first and second modes are set the same (d1 d2 = d) and their unnormalized columns are
A, B (M, - 0.5p, 0.55,, O.....0) 0) A. B. 00,0,0,0, (00.00000000000.:.4.0....... 0.5,, -0.5p, O.... 0
8
The third factor matrix with unnormalized & unshuffled columns is generated using
C2 (n....ll "i... T, C:2 --y... -I My..., M, U4.... [ds/2]
[da/2]
[ds/4]
[da/2]
[ds/4]
Then, the rows of C are shuffled randomly. There are four clusters over the third mode with cluster means (p), (J, p.) (-u, p)), (4,,-u) respectively. After normalizing the columns of A, B, C, we calculate the weights w;. Even though this example is inspired from Sun and Li (2I), shuffling the rows of C leads to additional clustering challenges since more sophisticated regularization is required to capture the clustering features. reflects cluster difficulty as it determines the distance between different clusters. A smaller leads more challenging clustering problem. For simplicity, the factor matrices used in this example, A, B, C, are orthogonal, satisfying assumption AL.
We set d = 20, d3 48, 1, vary (0, 0.25, 0.5, 0.75, }} and compare the performance of the four algorithms by reporting recovery and clustering errors. As summarized Figure 2 (a) and outperforms the other three methods under different values of perturbation,
(b), Fused-Orth-Ae achieving the best recovery error and clustering error. Not surprisingly. TBM performs worst due to model misspecification: the formulation of simulated factor matrices leads to underlying cluster structure for mode three instead of a multiway block structure for all three modes. Figure 2 and (d) show that average rumning times are comparable for Fused- fused-Orth-a Orth-ALS and DTC where the orthogonalization step added in Sused-OOtt-AAL. still decrease number of iterations required for convergence. Furthermore, to illustrate the advantage of generalized LASSO regularization imposed on Missed-Orth-A.. over fused LASSO penalty adopted in DTTC, we vary the value for E {0.1, 0.7) to mimic different levels of clustering complexity and set d 8, 0.001, ds 20. It is obvious from Table T that our method outperforms DTC especially as the 'signal' y becomes weak.
(c)
Table 1: Comparison of Fused-Orth-ALS and DTC under different signal levels. Average emoor and standard deviations (in parenthesis) are reported based on 50 replications.
/ = 0.1
1 0.7
Fused-Orth-AI.S
DTC
Fused-Orth-ALS
DTC
Recovery error Clustering error
(..034))
1. 1.0538(0.0224)
0.0011(0.0001) ((0)
0.0069(0.0006)
0.175)
0.3847(0.0394)
0(0)
7

---

(a)
(b)
(c)
(d)
1.00
150-
Running time (sec)
Recovery Emor
Clusterin Emro
0.75
0.3
Iterations
0.50
0.2
0.25
50-
0.00
0.00
0.25
0.75 sigma
0.50
1.00
0.00
0.25
0.50
0.75
1.b0
0.00
0.25
0.75 sigma
0.50
1.00
0.75 sigma
0.50
1.50
sigma
CP-Kne DTC
Fused-Orth-ALL TBM
Figure 2: Comparison of single-mode clustering for Fused-Orth-ALS, DTC CP-Kmeans and TBM
5.1.2 Multi-mode clustering
To validate the performance multi-mode clustering, we extend the simulation settings above by imposing clustering structures on all three factor matrices A, B, C, similarly to (3). The corresponding choice of cluster means are: (u, 0.5,), (-0.5 (u), (0, -u), (-4, 0) for mode 1, (0, p)) (-,, 0), (0.1μ, -p), (-μ, 10.14) for mode and ( (), (-,, μ), (14, -1), (-H, for mode 3. Accordingly, each mode has four clusters, i.e., $1 82 $3
of same size: [d./4]. d22/4], |ds/ / respectively. For example, for the first mode, there are |d//4 rows of the factor matrix taking value (μ, 0.5j), [d//4] rows taking value (-0.5 p), [d1/4] rows taking
value (0, -4) and the rest rows taking value 0). We fix dj = d2 = 20, dg = 40, μ and vary E (1, 2, 5, Since the detailed comparison between Fused-Orth-AL.S and DTC has been illustrated above, we omit DTC (originally designed for dynamic single mode clustering of tensors). It is clear that Fused-Orth-ALS algorithm outperforms
Comparison results can be found in Figure 3
the other two methods, especially under the noisy case when ar is large.
M Mode>
Mode 2
Mode 3
0.3
0.3-
0.3-
Clustering emor
Clustering emor
Clustering emor
0.2
0.1
0.1
0.1-
2
sigma
sigma
sigma
CP -means Fused-Orth-ALSS TBM
Figure 3: Comparison of multi-modes clustering for Fused-Orth Oatt-OrtthhSS.t CP-Kmeans and TBM
5.2 Real datasets
Experiments are conducted on two real datasets: the brain node structural connectivity from Human
Connectome Project (HCP) | [25 and political relationships between nations 13 Details for choosing rank K and number of clusters are provided supplementary material.
The HCP dataset is formatted as ) E 1 and consists of brain connectivity among 68 brain nodes for 136 individuals. Each entry takes on ordinal value {0, 1,2] which indicates the strength level of connectivity {low, moderate, high} between different brain nodes. Clustering results in Table 2 capture spatial connectivity between hemispheres of brain (the first character in the node name indicates the left right hemisphere and the number in the parenthesis indicates the node count with same name). Cluster and mainly capture connectivity either left or right hemisphere while
at hattp://wre www.anccoometttmmpooeet......... Dataset is wttpapat.aittttr/aaaoaoooooononn
org
Dataset is available hattp://wwe wwwrleskeep.p.oooooo com/
code/ wwwp://p/cccccccccctttttccccccct

---

cluster III represents the cross-section connection between left and right hemisphere. rsupramarginal and Lsupramarginal are picked as smaller clusters those regions are known to play a critical role in visual word recognition and reading.
Table 2: Clustering result for 68 brain nodes in HCP dataset
Cluster
Brain Nodes
II
III
..cueeus.r.ppraahpppccaapp
IV
rsupraamrrinal(4)
V
Isupramarginal(4)
The Nations dataset is formatted as y E 141x14455 consisting of 56 political relationships of 14 countries between 1950-1965. Each entry represents the presence or absence of a political action, such as 'treaties', 'send tourists to between different nations. Since 78.9% entries are zero in this dataset, we include an &o penalized tensor block model (denoted The clustering result in for Stseed-Ott-AAa algorithm assigned the nations into three clusters, one representing western-bloc countries (Cluster 1), one communist bloc (Cluster III), one neutral bloc (Cluster II). This result is consistent with the structure of political environment after world war II. addition, the goodness-of-fit (proportion of variance explained) for five methods is provided in Table 4
"TBM--passe" )) for comparison.
For the nations dataset, Fused-Orth-ALL Orth-ALS algorithm shows the highest variance explained.
Table 3: Clustering result for 14 nations in Nations dataset
Cluster
Characteristic
Country
Cluster
Western
Brazil, Netherlands, UK, USA Burma, Egypt, Israel, India, Indonesia China, Cuba, Poland, USSR
Cluster II
Neutral
Cluster Communist
Table 4: Comparison of goodness-of-fit for HCP and nations dataset
Fused-Orth-ALS.
DTC
TBM
TBM-Sparse
CP-Kmeans
HCP
0.921
0.921
0.925
0.921
Nations
0.522
0.458
0.439
0.433
0.324
6 Conclusions
This paper studies the properties and convergence rates for a novel multiway tensor clustering algo- rithm, Frused-OrthAAL. The underlying CP decomposition structure imposes low-rank constraints on the tensor observation, and graph similarity based fusion regularization encourages smoothness on factor matrices, thus leading to automatic clustering. The Fused-Orth- sparse and dense data tensors, and achieves clustering consistency even under model misspecification. An open problem developing theory for rank selection consistency. Another interesting topic is on clustering consistency if there are missing entries in the data. We leave these directions for further study. There are foresecable negative social impacts of this work.
algorithm can handle

---

References
[1] Genevera Allen. Sparse higher-order principal components analysis. In Artificial Intelligence and Statistics. pages 27-36, 2012.
[2] Taylor Amold and Ryan Tibshirani. Efficient implementations of the generalized lasso dual path algorithm. Journal of Computational and Graphical Statistics, 25(1):1-27, 2016.
Xuan Annic Qu, Junhui Wang, and Xiaotong Shen. group-speciffe recommender system. Journal of the American Statistical Association, 1344519) 1353, 2017.
[4] Norris Bruce BPS Murthi, Ram C Rao. dynamic model for digital advertising: The effects of creative format, message content, and targeting on engagement. Journal of marketing research, 54(2): 202-218, 2017.
[5] Gary K. Chen, Eric C. Chi, John M.O. Ranola, and Kenneth Lange. Convex clustering: An attractive alternative to hierarchical clustering. PLoS Computational Biology, 11()ee1004225 2015. doi: 10.1371/ Prrnmal.pbii.00000
[6] Eric C. Chi and Kenneth Lange. Splitting methods for convex clustering. Journal of Computational and Graphical Statistics, 24(4):99-1013, 2015. doi:
[7] Eric C. Chi, Genevera I. Allen, and Richard G. Baraniuk. Convex biclustering. Biometrics, 73(1):10- 19, 2017. doi: 10.1111/biom. 12540. URL http:ttpp///dddo org/ 100 1111/bicm. 12540
Eric C Chi, Brian J Gaines. Will Wei Sun, Hua Zhou, and Jian Yang. Provable convex co-clustering of tensors, Journal of Machine Learning Research, 21(214), 2020.
Victoria Hore, Ana Viñuela, Alfonso Buil, Julian Knight, Mark I McCarthy, Kerrin Small, and Jonathan Marchini. Tensor decomposition for multiple-tissue gene expression experiments. Nature genetics, 48(9): 1094, 2016.
[10] Furong Huang, Sergiy Matusevych, Anima Anandkumar, Nikos Karampatziakis, and Paul Mineiro. Distributed latent dirichlet allocation via tensor factorization, In NIPS Optimization Workshop, volume 1, 2014.
[11] Lawrence Hubert and Phipps Arabie. Comparing partitions. Journal of classification, 2(1):193-18,, 1985.
[12] U Kang, Evangelos Papalexakis, Abhay Harpale, and Christos Faloutsos. Gigatensor: scaling tensor analy- sis up by 100 times-al algorithms and discoveries. In Proceedings of the 18th ACM SIGKDD international conference Knowledge discovery and data mining, pages 316 324,,3244555 2012.
[13] Charles Kemp, Joshua B Tenenbaum, Thomas L Griffiths, Takeshi Yamada, and Naonori Ueda. Learning systems of concepts with an infinite relational model. In AAAI, volume 3, page 5, 2006.
[14] Tamara G Kolda and Brett W Bader. Tensor decompositions and applications. SIAM review, 51(3):455-500 500, 2009.
[15] Jean Kossaifi, Yannis Panagakis, Anima Anandkumar, and Maja Pantic. Tensorly: Tensor learning in python. The Journal of Machine Learning Research, 20(1):925 930, 2019.
[16] Tianqi Liu, Ming Yuan, and Hongyu 7hao. Characterizing spatiotemporal transcriptome of human brain via low rank tensor decomposition. arXiv preprint urXiv:1702.07449 2017.
[17] Marina Meilä. Comparing clusterings- an information based distance. Journal of multivariate analysis, 98(5):873 895, 2007.
[18] Maximilian Nickel, Volker Tresp, and Hans-Peter Kriegel. A three-way model for collective learning on multi-relational data. In /cml, volume 11. pages 809-816, 2011.
[19] Vatsal Sharan and Gregory Valiant. Orthogonalized als: theoretically principled tensor decomposition algorithm for practical use. In Proceedings of the 34th International Conference on Machine Learning- Volume 70, pages
JMLR. org, 2017.
[20] Yiyuan She et al. Sparse regression with exact clustering. Electronic Journal of Statistics, 4: 4:1055-1096. 2010.
[21] Will Wei Sun and Lexin Li. Dynamic tensor clustering. Journal of the American Statistical Association. pages 1-28, 2019.
10

---

[22] Kean Ming Tan and Daniela Witten. Statistical properties of convex clustering. Electronic journal of statistics, 9(2)23324 2015.
[23] Robert Tibshirani, Guenther Walther, and Trevor Hastie. Estimating the number of clusters in a data set the gap statistic. Journal of the Royal Statistical Society: Series (Statistical Methodology), 63(2):
41 1-23, 2001.
[24] Ryota Tomioka and Taiji Suzuki. Spectral norm random tensors. arXiv preprint arXiv:1407.187. 2014.
[25] David Van Essen, Stephen M Smith, Deanna Barch, Timothy EJ Behrens, Yacoub, Kamil Ugurbil. Wu-Minn HCP Consortium, et al. The wu-minn human connectome project: an overview. Neuroimage, 80: 62-79, 2013.
[26] Miaoyan Wang and Yuchen Zeng. Multiway clustering via tensor block models. In Advances in Neural Information Processing Systems, pages 713-723, 2019.
[27] Miaoyan Wang, Jonathan Fischer, Yun S Song, et al. Thyee-way clustering of multi-tissun multi-individul gene expression data using semi-nonnegative tensor decomposition. The Annals of Applied Statistics, 13 12):1103-1127. 2019.
[28] Yanqing Zhang, Xuan Bi, Niansheng Tang, and Annic Qu. Dynamic tensor recommender systems. arXiv preprint ar N0v:2003.00565 2020.
[29] Yunzhang Zhu, An augmented admm algorithm with application to the generalized lasso problem. Journal of Computational and Graphical Statistics, 26( 26(1):195-20) 2017.
Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's
contributions and scope? [Yes] See Section [
(b) Did you describe the limitations of your work? [Yes] See Section 6
(c) Did you discuss any potential negative societal impacts of your work? [Yes] See
Section
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [Yes] See Section 4
(b) Did you include complete proof' of all theoretical results? [Yes] See supplementary materials.
3. If you ran experiments...
(a) Did you include the code, data, and instructions needed to reproduce the main experi- mental results (either in the supplemental material or as a URL)? [Yes] See Section and supplementary material
(b) Did you specify all the training details (c.g., data splits, hyperparameters, how they were chosen)? [Yes] See Section 5 and supplementary material
(c) Did you report error bars (e.g., with respect the random seed after running experi- ments multiple times)? [Yes] See Section
(d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)? [Yes] See supplementary material.
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
(a) If your work uses existing assets, did you cite the creators? [Yes] See Section
(b) Did you mention the license of the assets? [Yes] See Section 3
(c) Did you include any new assets either in the supplemental material or as URL? [Yes]
See supplementary material.
(d) Did you discuss whether and how consent was obtained from people whose data you're using/curating? [Yes] See Section

---

(c) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [N/A]
5. If you used crowdsoureing conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, if applicable? [N/A]
(b) Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, applicable? [N/A]
(c) Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation? [N/A]
12


Reviews:
[{"review_id": "OnYlF6nvpIt", "reviewer": "Reviewer_c3Lt", "summary": "This paper proposes a new formulation and a new ALS based algorithm for tensor clustering. Upper bound of the theoretical convergence rate of this algorithm is provided, and experimental results show that this algorithm yields more accurate clustering for multiple cases compared to existing methods.", "questions": "\nLine 70: spectrum norm formulation is wrong, I guess the authors mean max rather than min.\n\nLine 79: There is one redundant ) at the end of this line.", "limitations": "There are no negative social impacts of this work.", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\n- The problem of tensor clustering is important\n- Detailed theoretical analysis of the algorithm, and better experimental results compared to existing work.\n\nWeaknesses:\n\n- As is declared in the paper: \"First, the orthogonalization step is performed before each iteration of ALS, which allows for the avoidance of local optima and more rapid convergence to the true factors\". There is no detailed analysis in the paper to show the advantage of this orthogonalization step. Additional experiments that compare the algorithm with orthogonalization and the one without orthogonalization would be helpful.\n- The formulation (1) includes multiple hyperparameters. The choice of these hyperparameters in the paper is presented in 3.2 However, questions still remain as to 1) is the formulation robust to these hyperparameters, and 2) is there other reasonable ways to choose these hyperparameters, and discussions over them would be helpful.\n\n-------\nAfter rebuttal: I thank authors for the detailed response. My concerns are addressed.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "JaR9eozPdWx", "reviewer": "Reviewer_cTv9", "summary": "This paper tackles the tensor multi-mode clustering problem. The authors model the tensor to cluster as the sum of two terms, the signal which admits a rank-K CANDECOMP/PARAFAC (CP) decomposition, and some noise. They establish a link between the CP decomposition factors and the means of the clusters from the corresponding modes to propose the core optimization problem described in the paper. They aim to approximate as closely as possible the tensor by a rank-K CP decomposition, with a generalized Lasso penalty that limits the number of unique rows in the CP decomposition factors. An Alternating Least Squares (ALS) strategy is used to solve the optimization problem one factor at a time. An orthogonalization step is added before each ALS step, and a Fuse operator is applied after the ALS step.\nThe algorithm is presented with convergence guarantees, and the proposed approach is compared with other tensor mode clustering methods on both simulated and real datasets. The authors also explain how to choose every parameter of the proposed algorithm.", "questions": "I start with an essential question determining my acceptance of the paper.\nPlease correct me if I am wrong, but it seems that the proof of Theorem 1 bites its tail. Indeed, Theorem 1 gives an upper bound on the distance between the true and the estimated factors. However, its proof is based on Lemma 3, which assumes this same upper bound. Furthermore, it is written in Lemma 3 that this assumption can be made \"without loss of generality\". I may be missing something here, so please enlighten me.\n\nI will now ask some questions, give suggestions to improve readability and point out mistakes or typos.\n- Could it be possible to assess the validity of the A1-A4 assumptions on a given dataset? Were they all valid on the simulated datasets?\n- General remark: it would help the reader if different indices were used for the rows and the columns. For instance, use $k$ for columns and $i$ for rows: $\\mathbf A_{:k}$ and $\\mathbf A_{i:}$.\n- Line 71: the definition of the spectral norm should be $ \\Vert \\mathcal Y \\Vert = \\text{max}_{\\Vert \\mathbf u_1 \\Vert_2 = \\Vert \\mathbf u_2 \\Vert_2 = \\Vert \\mathbf u_3 \\Vert_2 = 1} \\mathcal Y(\\mathbf u_1, \\mathbf u_2, \\mathbf u_3)$.\n- Line 79: \"rank-K CP decomposition\" and no closing parenthesis at the end of the line.\n- Line 80: \"unit ~~one~~ norm\".\n- Line 121: it should be [...] $+ \\lambda \\Vert \\boldsymbol \\Delta \\mathbf u \\Vert_1$. $\\boldsymbol \\Delta$ is also implicitly defined as specific matrices per mode were previously introduced.\n- Line 169: $\\gamma$ is used for the signal ratio, there is no possible confusion with the regularization parameters, but maybe another symbol should be preferred.\n- Line 170: in A4, constant $M$ has not been introduced before.\n- Line 172: \"on **the** rank $K$\".\n- Line 191: observation made that large $\\gamma$ results in lower error bounds seems wrong according to the bound given in Theorem 1.\n- Line 200: what does the condition on $w_\\text{min}$ means? Is it a high signal-to-noise ratio condition? Was it verified in the simulations?\n- Lines 256-257: the motivation for adding orthogonality was, according to line 119, to get more rapid convergence, but here, a concern is addressed that this may increase the number of iterations. Talking about verification in lines 256-257 would be more logical.\n\n- Appendix line 65, on the second line of the equation defining $\\Delta$ (maybe another symbol would be a good choice not to get confused with $\\boldsymbol \\Delta$): why does this become an upper bound with $\\mathbf A_{: i}$ replaced with $(\\mathbf A_{: i} + \\hat \\xi_i)$?\n- Appendix line 68 (and after): $II_{11}'$ is used but never introduced. It happens after with different quantities.\n- Appendix line 71, on the first line of the equation: $\\mathcal Y^*$ should be replaced with $w_l$.\n- Appendix line 73, on the first of the equation: $w_l$ is missing.\n- Appendix line 79: indices are wrong. It should be $j_1$ for the $\\mathbf A$ terms and $j_2$ for the $\\mathbf B$ terms. I get the idea of the majorization, but it is not well-written: the sums have been removed, but there are still terms in $j_1$ (and normally $j_2$) appearing in the upper bound.\n- Appendix line 83, equation 11: the term $2 w_\\text{max} \\rho$ supposes that $j_1 \\neq j_2$, but we could have the equality.\n- Appendix line 84, the whole bound should be multiplied by $((K - 1) \\Delta)^2$ and the last term is $(K - 2) \\rho^2 w_\\text{max}$.\n- Appendix line 146: $\\Vert . \\Vert_2$.\n- Appendix line 149 (and after): $\\alpha^2$ turned into $\\alpha$.\n- Appendix line 155: I do not get how the upper bound is simplified.\n- Appendix line 156: I do not get the bounds on $| 1 / \\kappa |$, and I do not get how the upper bound is simplified.", "limitations": "The authors identified that their convergence guarantees were valid for the real rank $K$ and that more work is needed to derive results for an estimated rank $K$. The convergence results show that if the number of samples goes to infinity and the tensor is noiseless, the true factors will be retrieved. However, it does not say anything about the algorithm reaching a stationary point. Is it possible to show, for example, that the objective function decreases monotonously?", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "### Strengths\n1. The paper is clear and well written. The progression is logical and the choices made are well justified.\n\n2. The proposed approach is compared favorably with concurrent methods on simulated and real datasets.\n\n3. Theoretical convergence properties of the given algorithm are derived.\n\n### Weaknesses\n4. I have concerns about the proof of Theorem 1. See the section **Questions** below.\n\n5. No code is given to reproduce the results.\n\n**Update after rebuttal**: the authors have provided code (although I did not have time to go through it) and addressed my concerns regarding the proof of Theorem 1. I think that the clarity of the proof can still be improved, but the authors are working on it, and I have good hope that they will manage to do so in the given amount of time.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "I2VW2HAMJ3K", "reviewer": "Reviewer_UQxe", "summary": "This work proposes an algorithm for simultaneous tensor decomposition and clustering using the alternating least squares approach. The proposed method is accompanied by a complexity and convergence analysis. Experimental studies on both synthetic and real-world datasets demonstrate the effectiveness of the proposed method in clustering multiway array data.", "questions": "1. Why is a CP model chosen over other tensor decomposition models? Please elaborate.\n2. Is the CP assumption relevant on real-world datasets?\n3. The paper claims that employing orthogonalization in ALS avoids local minima and speeds up convergence. Although empirical evidence of this general statement it provided in previous works, does this hold in this setting where ALS is employed for clustering? Please provide intuition.\n\n", "limitations": "The paper discusses a couple of limitations of this work in section 6 and aims to address them in future work. ", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n1. Paper is well-written and easy to follow.\n2. The proposed method is well motivated.\n3. The proposed method is presented as a concise pseudocode and it is easy to follow.\n4. Complexity analysis shows the proposed method's cost is similar to that of ordinary ALS, showing that it is practically applicable.\n5. Section 3.2 on how to choose regularization weights shows that these vital hyperparameters are chosen systematically from the data and not heuristically.\n\nWeakness:\n1. Missing detailed intuition behind the reason for choosing a CP model as opposed to a Tucker model for example.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept