Splitting Methods for Convex Bi-Clustering and Co-Clustering

TL;DR

This paper introduces three efficient operator-splitting algorithms for convex bi-clustering and co-clustering, improving computational efficiency and providing theoretical analysis for large structured data clustering tasks.

Contribution

It presents novel operator-splitting methods, including a Generalized ADMM, tailored for convex co-clustering, with demonstrated efficiency and theoretical complexity analysis.

Findings

Generalized ADMM outperforms other methods on large problems

Theoretical complexity analysis supports experimental results

Operator-splitting methods are effective for convex co-clustering

Abstract

Co-Clustering, the problem of simultaneously identifying clusters across multiple aspects of a data set, is a natural generalization of clustering to higher-order structured data. Recent convex formulations of bi-clustering and tensor co-clustering, which shrink estimated centroids together using a convex fusion penalty, allow for global optimality guarantees and precise theoretical analysis, but their computational properties have been less well studied. In this note, we present three efficient operator-splitting methods for the convex co-clustering problem: a standard two-block ADMM, a Generalized ADMM which avoids an expensive tensor Sylvester equation in the primal update, and a three-block ADMM based on the operator splitting scheme of Davis and Yin. Theoretical complexity analysis suggests, and experimental evidence confirms, that the Generalized ADMM is far more efficient for large problems.

Tables

Table 1. Table 1 : Performance of the ADMM, Generalized ADMM, Davis-Yin Splitting, and COBRA on the presidential speech and TCGA data sets. For COBRA, the numbers in parentheses are the total number of sub-problem iterations taken. The † † \dagger indicates that the method was stopped after failing to converge in 10,000 iterations.

Method	Presidential Speeches			TCGA Breast Cancer
	Iter/sec	Total Time (s)	Total Iterations	Iter/sec	Total Time (s)	Total Iterations
Generalized ADMM	$663$	3.98	2639	$2.58$	1318	3411
AB Accelerated	$645$	4.69	3025	$2.55$	835	2135
Davis-Yin / AMA	$657$	†15.23	†10000	$2.60$	†3850	†10000
AB Accelerated	$653$	†15.31	†10000	$2.58$	†3879	†10000
ADMM	$189$	8.12	1536	$0.84$	1521	1277
AB Accelerated	$188$	7.08	1334	$0.84$	1324	1107
COBRA	$5.06$	37.20	188 (29581 total)	$0.054$	5994	325 (16519 total)
AB Accelerated	$6.89$	27.15	187 (21362 total)	$0.065$	4953	323 (13549 total)

Equations

$$\hat{\mathbf{U}}=\underset{\mathbf{U}\in {\mathbb{R}}^{n\times p}}{\mathrm{argmin}}\frac{1}{2}\Vert \mathbf{X}-\mathbf{U}{\Vert }_F^2+\lambda \sum _{\begin{array}{c}i,j=1\\ i<j\end{array}}^n{w}_{ij}\Vert {\mathbf{U}}_{i\cdot }-{\mathbf{U}}_{j\cdot }{\Vert }_q$$

$$\hat{\mathbf{U}}$$

$$+\lambda (\sum _{\begin{array}{c}i,j=1\\ i<j\end{array}}^n{w}_{ij}\Vert {\mathbf{U}}_{i\cdot }-{\mathbf{U}}_{j\cdot }{\Vert }_q+\sum _{\begin{array}{c}k,l=1\\ k<l\end{array}}^p{\tilde{w}}_{kl}\Vert {\mathbf{U}}_{\cdot k}-{\mathbf{U}}_{\cdot l}{\Vert }_q)$$

$$\hat{\mathcal{U}}=\underset{\mathcal{U}}{\mathrm{argmin}}\frac{1}{2}\Vert \mathcal{X}-\mathcal{U}{\Vert }_F^2+\lambda \sum _{j=1}^J\sum _{\begin{array}{c}k,l=1\\ k<l\end{array}}^{n_J}{w}_{k,l}^j\Vert {\mathcal{U}}^{j/k}-{\mathcal{U}}^{j/l}{\Vert }_q$$

$$\underset{(\mathbf{u},\mathbf{v})\in {\mathcal{H}}_1\times {\mathcal{H}}_2}{\mathrm{argmin}}f(\mathbf{u})+g(\mathbf{v})\text{ subject to }{\mathfrak{L}}_1\mathbf{u}+{\mathfrak{L}}_2\mathbf{v}=\mathbf{b}$$

$${\mathcal{L}}_{\rho }(\mathbf{u},\mathbf{v},\mathbf{z})$$

$$+\frac{\rho }{2}\Vert {\mathfrak{L}}_1\mathbf{u}+{\mathfrak{L}}_2\mathbf{v}-\mathbf{b}{\Vert }_{{\mathcal{H}}_{\ast }}^2.$$

$${\mathbf{u}}^{(k+1)}$$

$${\mathbf{v}}^{(k+1)}$$

$${\mathbf{z}}^{(k+1)}$$

$${\mathbf{u}}^{(k+1)}$$

$${\mathbf{v}}^{(k+1)}$$

$${\mathbf{z}}^{(k+1)}$$

$$\underset{(\mathbf{u},\mathbf{v},\mathbf{w})\in {\mathcal{H}}_1\times {\mathcal{H}}_2\times {\mathcal{H}}_3}{\mathrm{argmin}}$$

$${\mathfrak{L}}_1\mathbf{u}+{\mathfrak{L}}_2\mathbf{v}+{\mathfrak{L}}_3\mathbf{w}=\mathbf{b}.$$

$${\mathcal{L}}_{\rho }(\mathbf{u},\mathbf{v},\mathbf{w},\mathbf{z})=$$

$$+\frac{\rho }{2}\Vert {\mathfrak{L}}_1\mathbf{u}+{\mathfrak{L}}_2\mathbf{v}+{\mathfrak{L}}_3\mathbf{w}-\mathbf{b}-{\rho }^{-1}\mathbf{z}{\Vert }_{{\mathcal{H}}_{\ast }}^2$$

$${\mathbf{u}}^{(k+1)}$$

$${\mathbf{v}}^{(k+1)}$$

$${\mathbf{w}}^{(k+1)}$$

$${\mathbf{z}}^{(k+1)}$$

$$\hat{\mathbf{U}}$$

$${\mathbf{U}}^{(k+1)}$$

$$+\frac{\rho }{2}{\Vert \mathbf{U}{\mathbf{D}}_{\text{col}}-{\mathbf{V}}_{\text{col}}^{(k)}+{\rho }^{-1}{\mathbf{Z}}_{\text{col}}^{(k)}\Vert }_F^2$$

$$\left(\begin{array}{c}{\mathbf{V}}_{\text{row}}^{(k+1)}\\ {\mathbf{V}}_{\text{col}}^{(k+1)}\end{array}\right)$$

$$\left(\begin{array}{c}{\mathbf{Z}}_{\text{row}}^{(k+1)}\\ {\mathbf{Z}}_{\text{col}}^{(k+1)}\end{array}\right)$$

$$\mathbf{X}+\rho {\mathbf{D}}_{\text{row}}^T({\mathbf{V}}_{\text{row}}^{(k)}-{\rho }^{-1}{\mathbf{Z}}_{\text{row}}^{(k)})+\rho ({\mathbf{V}}_{\text{col}}^{(k)}-{\rho }^{-1}{\mathbf{Z}}_{\text{col}}^{(k)}){\mathbf{D}}_{\text{col}}^T$$

$$=\mathbf{U}+\rho {\mathbf{D}}_{\text{row}}^T{\mathbf{D}}_{\text{row}}\mathbf{U}+\rho \mathbf{U}{\mathbf{D}}_{\text{col}}{\mathbf{D}}_{\text{col}}^T.$$

$$\Vert {\mathfrak{L}}_1{\Vert }_{{\mathcal{H}}_1\to {\mathcal{H}}_{\ast }}\le {\sigma }_{\max }({\mathbf{D}}_{\text{row}})+{\sigma }_{\max }({\mathbf{D}}_{\text{col}})$$

$${\mathbf{U}}^{(k+1)}$$

$$+\rho ({\mathbf{V}}_{\text{col}}^{(k)}-{\rho }^{-1}{\mathbf{Z}}_{\text{col}}^{(k)}-{\mathbf{U}}^{(k)}{\mathbf{D}}_{\text{col}}){\mathbf{D}}_{\text{col}}^T)/(1+\alpha )$$

$${\mathbf{U}}^{(k+1)}$$

$${\mathbf{Z}}_{\text{row}}^{(k+1)}$$

$${\mathbf{Z}}_{\text{col}}^{(k+1)}$$

$$\hat{\mathcal{U}}=\underset{\mathcal{U}}{\mathrm{argmin}}\frac{1}{2}\Vert \mathcal{X}-\mathcal{U}{\Vert }_F^2+\lambda \sum _{j=1}^J\Vert \mathcal{U}{\times }_j{\mathcal{D}}_j{\Vert }_{j,q}$$

$${\mathcal{V}}_j^{(k+1)}$$

$${\mathcal{Z}}_j^{(k+1)}$$

$$\mathcal{X}+\rho \sum _{j=1}^J({\mathcal{V}}_j^{(k)}-{\rho }^{-1}{\mathcal{Z}}_j^{(k)}){\times }_j{\mathcal{D}}_j^{\mathcal{T}}=\mathcal{U}+\rho \sum _{j=1}^J\mathcal{U}{\times }_j{\mathcal{D}}_j{\times }_j{\mathcal{D}}_j^{\mathcal{T}}.$$

$${\mathcal{U}}^{(k+1)}$$

Full text

**Index Terms— **convex clustering, convex bi-clustering, optimization, splitting-methods, ADMM

1 Introduction

1.1 Convex Clustering, Bi-Clustering, and Co-Clustering

In recent years, there as been a surge of interest in the expression of classical statistical problems as penalized estimation problems. These convex re-formulations have many noteworthy advantages over classical methods, including improved statistical performance, analytical tractability, and efficient and scalable computation. In particular, a convex formulation of clustering [1, 2, 3] has sparked much recent research, as it provides theoretical and computational guarantees not available for classical clustering methods [4, 5, 6]. Convex clustering combines a squared Frobenius norm loss term, which encourages the estimated centroids to remain near the original data, with a convex fusion penalty, typically the $\ell_{q}$-norm of the row-wise differences, which shrinks the estimated centroids together, inducing clustering behavior in the solution:

[TABLE]

Observations $i$ and $j$ are assigned to the same cluster if $\hat{\bm{U}}_{i\cdot}=\hat{\bm{U}}_{j\cdot}$.

For highly structured data, it is often useful to cluster both the rows and the columns of the data. For example, in genomics one may wish to simultaneously estimate disease subtypes (row-wise clusters of patients) while simultaneously identifying genetic pathways (column-wise clusters of genes). While these clusterings can be performed independently, it is often more efficient to jointly estimate the two clusterings, a problem known as bi-clustering. Building on (1), [7] propose a convex formulation of bi-clustering, given by the following optimization problem:

[TABLE]

for fixed $\lambda\in\mathbb{R}_{\geq 0}$, $\bm{w}\in\mathbb{R}^{\binom{n}{2}}_{\geq 0}$, and $\tilde{\bm{w}}\in\mathbb{R}^{\binom{p}{2}}_{\geq 0}$. As with standard convex clustering, the fusion penalties fuse elements of $\bm{U}$ together, but here the rows and columns are both fused together, resulting in a characteristic checkerboard pattern in $\hat{\bm{U}}$.

Extending this, [8] propose a convex formulation of co-clustering, the problem of jointly clustering along each mode of a tensor (data array). Their estimator is defined by the following optimization problem, where both $\mathcal{X}$ and $\mathcal{U}$ are order $J$ tensors of dimension $n_{1}\times n_{2}\times\cdots\times n_{J}$ and $\mathcal{X}^{j/i}$ denotes the $i$th slice of $\mathcal{X}$ along the $j$th mode:

[TABLE]

for fixed $\lambda\in\mathbb{R}_{\geq 0}$ and $\bm{w}^{j}\in\mathbb{R}^{\binom{n_{j}}{2}}_{\geq 0},j=1,\dots,J$, and where $\|\mathcal{X}\|_{q}$ denotes the $\ell_{q}$-norm of the vectorization of $\mathcal{X}$. While many algorithms have been proposed to solve Problem (1), there has been relatively little work on efficient algorithms to solve Problems (2) and (3). Operator splitting methods have been shown to be among the most efficient algorithms for convex clustering [9] and, as we will show below, they are highly-efficient for bi-clustering and co-clustering as well.

1.2 Operator Splitting Methods

Many problems in statistical learning, compressive sensing, and sparse coding can be cast in a ‘‘loss + penalty’’ form, where a typically smooth loss function, measuring the in-sample accuracy, is combined with a structured penalty function (regularizer) which induces a simplified structure in the resulting estimate to improve performance. The popularity of this two-term structure has spurred a resurgence of interest in so-called operator splitting methods, which break the problem into simpler subproblems which can be solved efficiently [10]. Among the most popular of these methods is the Alternating Direction Method of Multipliers (ADMM), which is typically quite simple to implement, has attractive convergence properties, is easily parallelizable and distributed, and has been extended to stochastic and accelerated variants [11, 12, 13].

The ADMM can be used to solve problems of the form

[TABLE]

where $\bm{u}$, $\bm{v}$ take values in arbitrary (real, finite-dimensional) Hilbert spaces $\mathcal{H}_{1}$, $\mathcal{H}_{2}$, $\mathfrak{L}_{i}:\mathcal{H}_{i}\to\mathcal{H}_{*},i=1,2$ are linear operators from $\mathcal{H}_{i}$ to a common range space $\mathcal{H}_{*}$, $\bm{b}$ is a fixed element of $\mathcal{H}_{*}$, and $f,g$ are closed, proper, and convex functions from $\mathcal{H}_{1}$ and $\mathcal{H}_{2}$ to $\overline{\mathbb{R}}=\mathbb{R}\cup\{\infty\}$. The ADMM proceeds by Gauss-Seidel updates of the augmented Lagrangian

[TABLE]

The ADMM steps are given by:

[TABLE]

While this level of generality is not typically required in applications, we will see below that judicious use of the general formulation is key to developing efficient ADMMs for the co-clustering problem.

While the simplified sub-problems are often much easier to solve than the original problem, they may still be impractically expensive. [14] consider a generalized ADMM, which augments the subproblems (\theparentequation.1-\theparentequation.2) with positive-definite quadratic operators $\mathfrak{A}$ and $\mathfrak{B}$ to obtain the updates:

[TABLE]

If $\mathfrak{A}$ and $\mathfrak{B}$ are appropriately chosen, the generalized subproblems (\theparentequation.1-\theparentequation.2) may be easier to solve than their standard counterparts.

The wide adoption of the ADMM has lead to many attempts to extend the ADMM to objectives with three or more terms. Recently, [15, Algorithm 8 ] proposed an operator-splitting method for the three-block problem:

[TABLE]

Like the standard ADMM, the Davis-Yin scheme alternately updates the augmented Lagrangian in a Gauss-Seidel fashion:

[TABLE]

The Davis-Yin iterates bear similarities to the standard ADMM and the Alternating Minimization Algorithm [16], both of which are special cases:

[TABLE]

Like the AMA, the Davis-Yin three block ADMM requires strong convexity of $f(\cdot)$. Unlike the standard ADMM, which allows $\rho$ to be arbitrarily chosen or even to vary over the course of the algorithm, both the AMA and Davis-Yin ADMM place additional constraints on $\rho$ to ensure convergence.

2 Operator-Splitting Methods for Convex Bi-Clustering

In this section, we present three operator-splitting methods for convex co-clustering. For simplicity of exposition, we focus only on the bi-clustering (second-order tensor) case in this section and suppress the fusion weights $\bm{w},\tilde{\bm{w}}$. The bi-clustering problem (2) then becomes

[TABLE]

where $\bm{D}_{\text{row}}$ and $\bm{D}_{\text{col}}$ are directed difference matrices implied by the (non-zero) fusion weights and $\|\cdot\|_{\text{row},q}$ and $\|\cdot\|_{\text{col},q}$ are the sums of the row- and column-wise $\ell_{q}$-norms of a matrix. The extension of our methods to to tensors of arbitrary order is natural and discussed in more detail in Section 5.

2.1 Alternating Direction Method of Multipliers

Because the convex bi-clustering problem (9) has three terms, it is not immediately obvious how a standard two-block ADMM can be applied to this problem. However, if we take $\mathcal{H}_{1}=\mathbb{R}^{n\times p}$ to be the natural space of the primal variable and $\mathcal{H}_{2}=\mathcal{H}_{\text{row}}\times\mathcal{H}_{\text{col}}$ to be the Cartesian product of the row- and column-edge difference spaces, we can obtain the ADMM updates:

[TABLE]

The $\bm{V}$- and $\bm{Z}$-updates are straightforward to derive and parallel those of the ADMM for the convex clustering problem [9, 17, Appendix A.1]. Additionally, we note that, due to the separable Cartesian structure of $\mathcal{H}_{2}$, the $\bm{V}$ and $\bm{Z}$-updates can be applied separately, and potentially in parallel, for the row and column blocks. As we will see, this separable structure holds for all algorithms considered in this paper and for tensors of arbitrary rank.

The $\bm{U}$-update is somewhat more difficult as it requires us to solve

[TABLE]

This is a Sylvester equation with coefficient matrices $\frac{1}{2}\bm{I}+\rho\bm{D}_{\text{row}}^{T}\bm{D}_{\text{row}}$ and $\frac{1}{2}\bm{I}+\rho\bm{D}_{\text{col}}\bm{D}_{\text{col}}^{T}$ and can be solved using standard numerical algorithms [18, 19]. Typically, the most expensive step in solving a Sylvester equation is either a Schur or Hessenberg decomposition of the coefficient matrices, both of which scale cubically with the size of the coefficient matrices. Since the coefficient matrices are fixed by the problem structure and do not vary between iterations, this factorization can be cached and amortized over the ADMM iterates, though, as we will see below, it is often more efficient to modify the $\bm{U}$-update to avoid the Sylvester equation entirely.

2.2 Generalized ADMM

To avoid the Sylvester equation in the $\bm{U}$-update, we turn to the Generalized ADMM framework of [14] and augment the $\bm{U}$-update with the quadratic operator $\mathfrak{A}(\bm{U})=\frac{1}{2}(\alpha\|\bm{U}\|_{F}^{2}-\rho\|\mathfrak{L}_{1}\bm{U}\|_{\mathcal{H}_{*}}^{2})$. The constant $\alpha$ must be chosen to ensure $\mathfrak{A}$ is positive-definite so that sub-problem (\theparentequation.1) remains convex. The smallest valid $\alpha$ depends on $\rho$ and the operator norm of $\mathfrak{L}_{1}$ which can be bounded above by

[TABLE]

where $\sigma_{\max}(\cdot)$ is the maximum singular value of a matrix. To avoid expensive singular value calculations, an upper bound based on the maximum degree of any vertex in the graph can be used instead [20, 9, Section 4.1].

With this augmentation, the $\bm{U}$-update (\theparentequation.1) becomes:

[TABLE]

We note that this can be interpreted as a weighted average of a standard update and the previous estimate of $\bm{U}$. This update does not require solving a Sylvester equation and can be performed in quadratic, rather than cubic, time. As we will see in Section 3, this update significantly decreases the computational cost of an iteration without significantly slowing the per iteration convergence rate, resulting in much improved performance.

2.3 Davis-Yin Splitting

The three-block structure of the convex bi-clustering problem (9) naturally lends itself to the Davis-Yin three-block ADMM, which has provable convergence due to the strongly convex squared Frobenius norm loss. We take $\mathcal{H}_{1}$ as before, but now split the row- and column-edge differences into two spaces ($\mathcal{H}_{2}=\mathcal{H}_{\text{row}}$ and $\mathcal{H}_{3}=\mathcal{H}_{\text{col}}$). The joint range space, which also contains the dual variable, is still the Cartesian product $\mathcal{H}_{*}=\mathcal{H}_{2}\times\mathcal{H}_{3}$. The Davis-Yin iterates have the same $\bm{V}$- and $\bm{Z}$-updates as before, but yield a simpler $\bm{U}$-update:

[TABLE]

Unlike the ADMM, the Davis-Yin updates require $\rho$ to be less than $1/\|\mathcal{L}_{1}\|_{\mathcal{H}_{1}\to\mathcal{H}_{*}}$. The same bound used for $\alpha$ in the Generalized ADMM can be used here.

We note that, due to the ‘‘Cartesian’’ structure of the linear operators $\mathcal{L}_{2}\bm{V}=(-\bm{V}_{\text{row}},\bm{0}_{\mathcal{H}_{\text{col}}})$ and $\mathcal{L}_{3}\bm{W}=(\bm{0}_{\mathcal{H}_{\text{row}}},-\bm{V}_{\text{col}})$, the Davis-Yin updates have the same separable structure as the ADMM and Generalized ADMM updates. In fact, it is not hard to see that the Davis-Yin updates can also be derived as updates of the Alternating Minimization Algorithm (AMA) [16]. The AMA is known to be equivalent to proximal gradient on the dual problem, which [8] recommended to solve Problem (3).

As [9] note for the convex clustering problem, Moreau’s decomposition [21] can be used to eliminate the $\bm{V}$-update by replacing the proximal operator with a projection onto the dual norm ball, giving the simplified $\bm{Z}$-updates

[TABLE]

This simplification lowers the computational cost of the Davis-Yin updates, but it does not improve their per iteration convergence which, as we will see in the next section, is not competive with the ADMM variants.

3 Simulation Studies

In this section, we compare the performance of the three algorithms considered on the presidential speeches data set ($n=44$, $p=75$) of [17] and log-transformed Level III RPKM gene expression levels from 438 breast cancer samples collected by The Cancer Genome Atlas (TCGA) ($n=438$, $p=353$) [22]. In addition to the three algorithms discussed above, we also compare to the COBRA algorithm of [7], an application of the Dykstra-Like Proximal Algorithm of [23] to Problem (2), which works by solving alternating convex clustering problems (1) on the rows and columns of $\bm{U}$ until convergence. For each method, we compare a standard version and a version using the acceleration techniques proposed in [13]. (For COBRA, acceleration was applied to both the row and column sub-problems.)

We show results for the rotationally-invariant $q=2$ penalty with $\lambda=1\times 10^{4}$ for the presidents data and $\lambda=1\times 10^{6}$ for the breast cancer data, though our results are similar for other penalty functions and values of $\lambda$. For the ADMM variants, we fixed $\rho=1$. The Davis-Yin step size and the Generalized ADMM coefficient ($\alpha$) were both set to twice the maximum degree of the row- or column-indicidence graph. The default sparse Gaussian kernel weights of the clustRviz package [24] are used for each data set.

As can be seen in Figure 1, the COBRA algorithm clearly has the fastest per iteration convergence, rapidly converging to the optimal solution. (The slower convergence near optimality appears to be an artifact of the inexact solution of the sub-problems rather than an inherent behavior of the DLPA.) The ADMM and Generalized ADMM exhibit relatively fast convergence, while the Davis-Yin iterates are the slowest to converge. The standard ADMM appears to have slightly faster convergence than the Generalized ADMM near initialization, but the two methods appear to eventually attain the same convergence rate. Acceleration provides a minor, but consistent, improvement for all methods except COBRA, where the acceleration is applied within each sub-problem.

The apparent advantage of COBRA disappears, however, when we instead consider total elapsed time, as shown in Figure 2, as the COBRA updates require solving a pair of expensive convex clustering sub-problems at each step. On a wall-clock basis, the Generalized ADMM performs the best for both problems, followed by the standard ADMM, COBRA, and Davis-Yin splitting in that order.

Practically, the Generalized ADMM and Davis-Yin updates have essentially the same cost and are both three to four times faster per iteration than the standard ADMM with caching. Precise numerical results are given in Table 1. These results suggest that our Generalized ADMM scheme attains the best of both worlds, achieving both the high per iteration convergence rate of a standard ADMM with the low per iteration computational cost of the Davis-Yin / AMA updates, consistent with our theoretical analysis in the next section.

4 Complexity Analysis

Both the ADMM and AMA are known to exhibit $O(1/K)$ convergence in general, [25, 26, 13], but the ADMM obtains a superior convergence rate of $o(1/K)$ for strongly convex problems [27]. This rate can be further improved to linear convergence if certain full-rank conditions are fulfilled by the problem-specific $\bm{D}_{\text{row}}$ and $\bm{D}_{\text{col}}$ matrices [14]. As [14] show, the Generalized ADMM achieves essentially the same convergence rate as the standard ADMM. Consequently, theory suggests that the two-block ADMM variants have superior per iteration convergence, consistent with our experimental results.

The computational cost of the Generalized ADMM and Davis-Yin $\bm{U}$-updates are essentially the same, both being dominated by matrix multiplications giving $\mathcal{O}(np\max(\#\text{rows},\#\text{columns}))$ complexity, where $\#\text{rows}$ and $\#\text{columns}$ are the number of non-zero row and column fusion weights respectively. For sparse weighting schemes, we typically have $\#\text{rows}=\mathcal{O}(n)$ and $\#\text{columns}=\mathcal{O}(p)$, yielding an overall per iteration complexity of $\mathcal{O}(\max(np^{2},n^{2}p))$ for the $\bm{U}$-update. A naive implementation of the ADMM update requires $\mathcal{O}((\#\text{rows})^{3}+(\#\text{columns})^{3})$ to solve the Sylvester equation, but if the initial factorization is cached, the per iteration complexity is again reduced to $\mathcal{O}(\max(np^{2},n^{2}p))$ under a sparse weight scheme.

To the best of our knowledge, a convergence rate for the DLPA, on which COBRA is based, has not been established. Experimentally, the DLPA exhibits very rapid convergence, consistent with known rates for Dykstra’s alternating projections algorithm. Despite this, the high per iteration cost makes it impractical on larger problems.

5 Extensions to General Co-Clustering

In Section 2, we restricted our attention to co-clustering tensors of order two, i.e., bi-clustering. In this section, we show how the ADMM, Generalized ADMM, and Davis-Yin (AMA) algorithms can be extended to the general problem of co-clustering order-$K$ tensors (3). As before, we can re-express the co-clustering problem using directed difference tensors $\mathcal{D}_{j}$ as

[TABLE]

where $\|\mathcal{X}\|_{j,q}$ is the sum of the $\ell_{q}$-norms of the vectorization of each slice of $\mathcal{X}$ along the $j$th mode. For all three methods, the $\mathcal{V}$ and $\mathcal{Z}$-updates are straightforward extensions of the bi-clustering case:

[TABLE]

As with the bi-clustering case, we have separate $\mathcal{V}$ and $\mathcal{Z}$ variables for each term in the fusion penalty (i.e., each mode of the tensor), each of which can be updated in parallel. For the standard ADMM, the $\mathcal{U}$-update requires solving the following tensor Sylvester equation, where $(\cdot)^{\mathcal{T}}$ denotes a transpose, with the direction implied by context:

[TABLE]

As above, the $\mathcal{U}$-update can be solved explicitly for our other methods:

[TABLE]

6 Discussion

We have introduced three operator-splitting methods for solving the convex bi-clustering and co-clustering problems: a standard ADMM, a Generalized ADMM, and a three-block ADMM based on Davis-Yin splitting. The Davis-Yin three-block ADMM was found to be equivalent to the AMA which is in turn equivalent to the dual projected gradient recommended by [8]. The standard ADMM achieves the fastest per iteration convergence, but requires solving a Sylvester equation at each step. The Generalized ADMM avoids the expensive Sylvester equation while still maintaining the rapid convergence of the ADMM. Unlike [9], we do not find that the AMA performs well, perhaps due to the more complex structure of the bi- and co-clustering problems. On the TCGA breast cancer data, the Generalized ADMM was able to find a better solution in three and a half minutes than the AMA could in over an hour. All three of our methods can be extended to co-cluster tensors of arbitrary order, with the advantages of the Generalized ADMM becoming even more pronounced as it avoids a tensor Sylvester equation.

Our ADMM and Generalized ADMM methods provide a unified and efficient computational approach to convex co-clustering of large tensors of arbitrary order, but several interesting methodological questions remain unanswered. In particular, the conditions under which the co-clustering solution path is purely agglomerative, enabling the construction of dendrogram-based visualizations of the kind proposed for convex clustering by [17], remain unknown. While not discussed here, our algorithms can naturally be extended for missing data handling and cross-validation, using a straight-forward extension of the Majorization-Minimization (MM) imputation scheme proposed by [7]. Our algorithms make convex bi- and co-clustering practical for large structured data sets and we anticipate they will encourage adoption and additional study of these useful techniques.

7 References

Appendix A Step Size Selection

The Generalized ADMM requires that the quadratic operator $\mathfrak{A}(\bm{U})=\alpha\|\bm{U}\|_{F}^{2}-\rho\|\mathfrak{L}_{1}\bm{U}\|_{\mathcal{H}^{*}}^{2}$ be positive-definite, which holds if $\alpha/\rho$ is greater than the operator norm of $\mathfrak{L}_{1}$. Similarly, the Davis-Yin ADMM requires that the step size $\rho^{-1}$ be greater the operator norm of $\mathfrak{L}_{1}$. (See the related discussion in Section 4.2 of [9].) While this operator norm is somewhat inconvenient to compute, an upper bound is easily found:

[TABLE]

where $\sigma_{\max}(\cdot)$ is the maximum singular value of a matrix. To avoid a potentially expensive SVD for large problems, we can instead take advantage of the special structure of $\bm{D}_{\text{row}}$ and $\bm{D}_{\text{col}}$ using known results about the eigenstructure of a graph Laplacian [20] which imply

[TABLE]

where the first maximum is taken over the pairs of connected vertices and the second is taken over all vertices. This bound can be computed with minimal effort and, in our experience, is sufficiently tight for most problems.

Appendix B Detailed Derivations

In this section, we provide detailed derivations of the ADMM, Generalized ADMM, and Davis-Yin algorithms for convex clustering discussed above. Unless otherwise stated, our notation is the same as that used in the main body of the paper. Given the large number of Hilbert spaces, norms, and inner products used in these derivations, we err on the side of explicitness rather than concision in our notation.

B.1 ADMM

The ADMM updates are

[TABLE]

Rescaling the dual variable by $\rho^{-1}$ ($\bm{Z}\to\rho^{-1}\bm{Z}$), these can be simplified to

[TABLE]

For the convex bi-clustering problem (9), we take:

• •

$\mathcal{H}_{1}=\mathbb{R}^{n\times p}$ equipped with the Frobenius norm and inner product

• •

$\mathcal{H}_{\text{row}}=\mathbb{R}^{|\mathcal{E}_{\text{row}}|\times p}$ and $\mathcal{H}_{\text{col}}=\mathbb{R}^{n\times|\mathcal{E}_{\text{col}}|}$ both equipped with the Frobenius norm and inner product

• •

$\mathcal{H}_{2}=\mathcal{H}_{*}=\mathcal{H}_{\text{row}}\times\mathcal{H}_{\text{col}}$ equipped with the inner product

[TABLE]

and the norm

[TABLE]

• •

$\mathfrak{L}_{1}:\mathcal{H}_{1}\to\mathcal{H}_{*}$ given by $\mathfrak{L}_{1}\bm{U}=(\bm{D}_{\text{row}}\bm{U},\bm{U}\bm{D}_{\text{col}})$

• •

$\mathfrak{L}_{2}:\mathcal{H}_{2}\to\mathcal{H}_{*}$ given by $\mathfrak{L}_{2}(\bm{V}_{\text{row}},\bm{V}_{\text{col}})=-(\bm{V}_{\text{row}},\bm{V}_{\text{col}})$

• •

$\bm{b}=\bm{0}_{\mathcal{H}_{*}}=(\bm{0}_{\mathcal{H}_{\text{row}}},\bm{0}_{\mathcal{H}_{\text{col}}})\in\mathcal{H}_{*}$

• •

$f(\bm{U})=\frac{1}{2}\|\bm{X}-\bm{U}\|_{F}^{2}$

• •

$g((\bm{V}_{\text{row}},\bm{V}_{\text{col}}))=\lambda\|\bm{V}_{\text{row}}\|_{q,1}+\lambda\|\bm{V}_{\text{col}}\|_{1,q}$

First, we consider the $\bm{U}$-update:

[TABLE]

using the separability of the squared $\mathcal{H}_{*}$ norm ($\|\cdot\|_{\mathcal{H}_{*}}^{2}$). This is fully smooth and so we take the gradient with respect to $\bm{U}$ to obtain the stationarity conditions:

[TABLE]

using the identity111See Equation (119) in the Matrix Cookbook: https://www.math.uwaterloo.ca/~hwolkowi/matrixcookbook.pdf.

[TABLE]

Solving for $\bm{U}$, we obtain

[TABLE]

The $\bm{V}$-updates are straight-forward to derive using the separable structure of $\mathcal{H}_{2}=\mathcal{H}_{\text{row}}\times\mathcal{H}_{\text{col}}$:

[TABLE]

If fusion weights are included, then they are reflected in the row- and column-wise proximal operators. Finally, the $\bm{Z}$-updates are standard:

[TABLE]

If we unscale $\bm{Z}\to\rho^{-1}\bm{Z}$, we recover the updates given in the main body of this note.

B.2 Generalized ADMM

To avoid the Sylvester equation, we recommend use of a Generalized ADMM, where the $\bm{U}$-update is augmented with the quadratic operator $\mathfrak{A}(\bm{U}-\bm{U}^{(k)})=\frac{\alpha}{2}\|(\bm{U}-\bm{U}^{(k)})\|_{F}^{2}-\frac{\rho}{2}\|\mathfrak{L}_{1}(\bm{U}-\bm{U}^{(k)})\|_{F}^{2}$ where $\alpha$ is fixed sufficiently large to ensure positive-definiteness. This augmentation gives:

[TABLE]

As before, we take gradients with respect to $\bm{U}$ to obtain the stationarity condition:

[TABLE]

Simplifying, we obtain:

[TABLE]

The $\bm{V}$ and $\bm{Z}$-updates are unchanged.

B.3 Davis-Yin

The Davis-Yin updates are given by:

[TABLE]

For the convex bi-clustering problem (9), we take:

• •

$\mathcal{H}_{1}=\mathbb{R}^{n\times p}$ equipped with the Frobenius norm and inner product

• •

$\mathcal{H}_{2}=\mathcal{H}_{\text{row}}=\mathbb{R}^{|\mathcal{E}_{\text{row}}|\times p}$ equipped with the Frobenius norm and inner product

• •

$\mathcal{H}_{3}=\mathcal{H}_{\text{col}}=\mathbb{R}^{n\times|\mathcal{E}_{\text{col}}|}$ equipped with the Frobenius norm and inner product

• •

$\mathcal{H}_{*}=\mathcal{H}_{2}\times\mathcal{H}_{3}$ equipped with the inner product

[TABLE]

and the norm

[TABLE]

• •

$\mathfrak{L}_{1}:\mathcal{H}_{1}\to\mathcal{H}_{*}$ given by $\mathfrak{L}_{1}\bm{U}=(\bm{D}_{\text{row}}\bm{U},\bm{U}\bm{D}_{\text{col}})$

• •

$\mathfrak{L}_{2}:\mathcal{H}_{2}\to\mathcal{H}_{*}$ given by $\mathfrak{L}_{2}\bm{V}=(-\bm{V},\bm{0}_{\mathcal{H}_{3}})$

• •

$\mathfrak{L}_{3}:\mathcal{H}_{3}\to\mathcal{H}_{*}$ given by $\mathfrak{L}_{3}\bm{W}=(\bm{0}_{\mathcal{H}_{2}},-\bm{W})$

• •

$\bm{b}=\bm{0}_{\mathcal{H}_{*}}=(\bm{0}_{\mathcal{H}_{2}},\bm{0}_{\mathcal{H}_{3}})\in\mathcal{H}_{*}$

• •

$f_{1}(\bm{U})=\frac{1}{2}\|\bm{X}-\bm{U}\|_{F}^{2}$

• •

$f_{2}(\bm{V})=\lambda\|\bm{V}\|_{\text{row},q}$

• •

$f_{3}(\bm{W})=\lambda\|\bm{W}\|_{\text{col},q}$

First, we consider the $\bm{U}$-update:

[TABLE]

This is fully smooth, so we take the gradient with respect to $\bm{U}$, using the identity222See Equation (101) of the Matrix Cookbook: https://www.math.uwaterloo.ca/~hwolkowi/matrixcookbook.pdf.

[TABLE]

to obtain the stationarity condition

[TABLE]

Solving for $\bm{U}$, we obtain the Davis-Yin $\bm{U}$-update:

[TABLE]

The $\bm{V}_{\text{row}}$-update (i.e., the $\bm{V}$-update above) is given by

[TABLE]

Using the separability of the squared $\mathcal{H}_{*}$-norm, we see that the $\mathcal{H}_{3}$ component of the second term does not depend on $\bm{V}$ and hence the above problem is equivalent to

[TABLE]

A similar analysis gives the $\bm{V}_{\text{col}}$-update (i.e., the $\bm{W}$-update above):

[TABLE]

Finally, the $\bm{Z}$-update is given by:

[TABLE]