###### Abstract

Biclustering is a powerful data mining technique that allows simultaneously
clustering rows (observations) and columns (features) in a matrix-format data
set, which can provide results in a checkerboard-like pattern for
visualization and exploratory analysis in a wide array of domains. Multiple
biclustering algorithms have been developed in the past two decades, among
which the convex biclustering can guarantee a global optimum by formulating in
as a convex optimization problem. On the other hand, the application of
biclustering has not progressed in parallel with the algorithm techniques. For
example, biclustering for increasingly popular microbiome research data is
under-applied possibly due to its compositional constraints for each sample.
In this manuscript, we propose a new convex biclustering algorithm, called the
bi-ADMM, under general setups based on the ADMM algorithm, which is free of
extra smoothing steps to visualize informative biclusters required by existing
convex biclustering algorithms. Furthermore, we tailor it to the algorithm
named biC-ADMM specifically to tackle compositional constraints confronted in
microbiome data. The key step of our methods utilizes the Sylvester Equation
to derive the ADMM algorithm, which is new to the clustering research. The
effectiveness of the proposed methods is examined through a variety of
numerical experiments and a microbiome data application.

A New Algorithm for Convex Biclustering and
Its Extension to the Compositional Data

arg

Key words: Compositional data; Convex biclustering; Fused lasso; Microbiome
data; Sylvester Equation

##  1 Introduction

Biclustering is a powerful data mining technique that allows simultaneously
clustering rows (observations) and columns (features) in a matrix-format data
set. It is a 2-way extension of clustering analysis which aims to assign
observations into a number of clusters such that observations in the same
group are similar to each other. Biclustering can provide results in a
checkerboard-like pattern, which can be used for visualization and exploratory
analysis in a wide array of domains. For example, Cheng and Church, (2000)
introduced the concept of biclustering for the first time to gene expression
data, aiming to overcome some problems of traditional clustering methods such
as information loss during oversimplied similarity and grouping computation.
Instead, biclustering identified genes expression patterns under a subset of
all the conditions/samples, which provided a better representation for genes
with multiple functions or regulated by many factors. Furthermore, in the big
data era, electronic health records contain abundant information that can be
transformed into disease phenotypes. These phenotype data can be organized
into a matrix, with individuals as rows and phenotype features as columns,
from which biclustering can discover a subgroup of patients from a subset of
phenotypes (Hripcsak and Albers,, 2012). Results from biclustering can shed
lights on the downstream analyses to interpret the identified biclusters and
to demonstrate bicluster associations coupled with other statistical
evaluations tool.

In the past two decades, multiple biclustering algorithms have been developed.
The most popular approach for genomic data is the clustered dendrogram
(Madeira and Oliveira,, 2004; Busygin et al.,, 2008), which actually performs
the hierarchical clustering (Friedman et al.,, 2001) on both subjects and
genes. However, dendrogram is less ideal for generating reproducible results
because it greedily fuses observations or features to minimize a specific
objective function such as the within-clustering variance. As a consequence,
this algorithm may return a locally optimal solution with respect to the
criterion, and dendrogram is not stable to small perturbations of the data due
to its greedy algorithm. Later on, more advanced biclustering methods have
been proposed, which are either based on the singular value decomposition
(SVD) or graph cuts. However, Chi et al., (2017) pointed out that none of
these methods addressed the fundamental issues that underlay the clustered
dendrogram. Alternatively, they formulated the biclustering problem as a
convex optimization problem, and solved it with an iterative algorithm based
on the splitting methods for convex clustering.

Existing biclustering methods have several ways to define biclusters: 1)
defining a submatrix with a large average as a bicluster (Shabalin et al.,,
2009); 2) defining a bicluster by coherent low-rank patterns (Lee et al.,,
2010; Li,, 2020); 3) assuming the elements of a bicluster which shares the
same underlying mean (Flynn and Perry,, 2012; Tan and Witten,, 2014; Chi et
al.,, 2017). To be specific, this manuscript adopts the third way as it is
motivated by Chi et al., (2017); Wang et al., (2018).

On the other hand, application of biclustering has not progressed in parallel
with the algorithmic development. Recently, Xie et al., (2018) provided a
comprehensive review of biclustering applications in biological and biomedical
data, and pointed out there is a need for developing supporting computational
techniques and providing the guidance on selecting the appropriate
biclustering tools in a specific study. Data generated from different
biotechnologies have their own properties. For example, 16S and shotgun
sequencing are two widely used sequencing techniques in microbiome studies.
The raw microbiome data are typically organized into large matrices with rows
representing samples, and columns containing observed counts of clustered
sequences commonly known as operational taxonomic units (OTUs). Because the
observed counts are not comparable across samples due to their constrained
total by the sequencing depth (Gloor et al.,, 2017), normalizing counts toward
their totals into the relative abundance is a commonly used metric in
microbiome studies, which imposes a compositional constraint for data from
each sample. As reviewed by Xie et al., (2018), there are limited applications
of biclustering in microbiome studies, except one example. Falony et al.,
(2016) identified sample subsets with specific taxonomic signatures using a
biclustering approach called FABIA, proposed by Hochreiter et al., (2010).

In this manuscript, we propose a new convex biclustering algorithm, named bi-
ADMM, under general setups based on the alternating direction method of
multipliers (ADMM) algorithm. Furthermore, we tailor it to the algorithm,
named biC-ADMM, specifically to tackle the compositional constraints
confronted in microbiome data. As far as we have known, it is the first
biclustering algorithm that can handle data with compositional constraints.
Compared with Chi et al., (2017), the novelty of this paper is four-fold.
First, we solve the convex biclustering problem via the ADMM algorithm
directly. As Chi et al., (2017) mentioned, although an ADMM algorithm was
feasible, they solved it in an alternative way as a convex bicluster ring
algorithm (COBRA) via iteratively using Dykstra-like proximal algorithm (DLPA)
combined with the alternating minimization algorithm (AMA) introduced in Chi
and Lange, (2015) as the convex clustering solver. In the meanwhile, their
methods require extra smoothing steps to facilitate the visualization of
estimtaed biclusters, seeing Figure 1(b) in Chi et al., (2017). The key step
of our methods utilizes the Sylvester Equation for the first time to solve the
challenging estimation problem involving constraints on both rows and columns
of a matrix. Second, the proposed algorithms can naturally incorporate
additional constraints for rows or columns. For example, the modified
algorithms for the microbiome data with compositional constraint for data rows
are provided. Third, the proposed algorithms use two tuning parameters to
control the resulting numbers of row and column clusters separately, which
includes the single parameter setting used by Chi et al., (2017) to controls
the biclustering path as a special case. As a result, our methods offer users
an opportunity to fix clusters for one dimenstion and monitor the clustering
path of the other one. Fourth, the fused-lasso penalty terms in the proposed
methods accommodate $L_{q}$-norm, $q=1,2,\infty$, and thus, are more flexible
than the COBRA which requires the $L_{2}$-norm. Note that our methods are not
only theoretically sound, but also practically promising. The superior
performance of our methods is demonstrated in extensive simulated examples and
the application to analyze a murine gut microbiome dataset (Livanos et al.,,
2016).

We demonstrate the superior performance of the proposed algorithm under the
simulation setting in Section 4.1 with $\sigma=4$, where $\sigma$ is the
standard deviation of the random noise for data generation. To ensure a fair
comparison, we adopt a special case of our methods using the objective
function (10) with a single tuning parameter, which is also used by COBRA. In
this setting, there are $n=80$ samples from 4 clusters and $p=40$ features
from 4 clusters. Figure 1 compares the performance of COBRA without the
smoothing step (indicated as COBRA on the left panel) and our method
(indicated as bi-ADMM on the right panel), by visualizing the estimate for
each method with the same set of weights and the same value for the single
tuning parameter, $\gamma=6.72$. The heat maps show that the proposed bi-ADMM
algorithm can directly visualize the estimate with a clear checkerboard-like
pattern without any further procedure.

Figure 1: The heat maps of estimated matrices based on the COBRA (left) and
the bi-ADMM (right) algorithms with the same tuning parameter $\gamma=6.72$,
respectively. Rows indicate 80 samples with 4 clusters and columns indicate 40
features with 4 clusters. The row(column) side color bars indicate the true
row(column) clusters.

The rest of the paper is organized as follows. We propose the main methods and
algorithms in Section 2. Some practical consideraions for applying the
proposed algorithms are discussed in Section 3. The performance of the
proposed methods via simulation studies are evaluated in Section 4, and a
microbiome application is presented in Section 5. We conclude the paper with a
brief summary in Section 6, and defer all technical proofs to the Appendix.

##  2 Convex Biclustering

In this section, we develop algorithms for a general biclustering problem
first, and then extend them to incorporate the compositional constraints, a
unique feature of the microbiome data.

###  2.1 Specification of the General Model

We consider a general biclustering problem as follows. Let
$\mathbf{X}\in\mathbb{R}^{n\times p}$ be a data matrix with $n$ observations
$X_{i\cdot}=(X_{i1},\ldots,X_{ip})^{\rm T}$ with $p$ features, $i=1,\cdots,n$.
The underlying assumption for applying a biclustering method is that the $n$
observations belong to $K$ unknown and non-overlapping classes,
$C_{1},\ldots,C_{K}$, and the $p$ features belong to $R$ unknown and non-
overlapping classes, $D_{1},\ldots,D_{R}$. Such checkerboard-like structure of
a data matrix motivates researchers to use the fused-lasso penalty (Tibshirani
et al.,, 2005) to fuse rows and columns simultaneously.

To facilitate further derivations, we rewrite the data matrix $\mathbf{X}$ in
feature-level as column vectors $\mathbf{x}_{j}$ so that
$\mathbf{X}=(\mathbf{x}_{1},\cdots,\mathbf{x}_{p})$, where
$\mathbf{x}_{j}=(X_{1j},\cdots,X_{nj}){{}^{\rm T}}$, $j=1,\ldots,p$. Similarly
we denote $\mathbf{A}$ in the feature-level consisting of column vectors
$\mathbf{A}=(\mathbf{a}_{1},\cdots,\mathbf{a}_{p})$, where
$\mathbf{a}_{j}=(A_{1j},\cdots,A_{nj}){{}^{\rm T}}$, and as
$(A_{1\cdot},\ldots,A_{n\cdot})^{\rm T}$ in the observation-level. Define
${\cal E}_{1}=\\{l=(l_{1},l_{2}):1\leq l_{1}<l_{2}\leq n\\}$ and ${\cal
E}_{2}=\\{k=(k_{1},k_{2}):0\leq k_{1}<k_{2}\leq p\\}$. Then denote $|{\cal
E}_{1}|$ and $|{\cal E}_{2}|$ as the numbers of components of ${\cal E}_{1}$
and ${\cal E}_{2}$, respectively.

We formulate the convex biclustering problem as the following minimization
problem with tuning parameters $\gamma_{1}$ and $\gamma_{2}$:

$\displaystyle\min_{\mathbf{A}\in\mathbb{R}^{p\times n}}$
$\displaystyle{1\over
2}\sum_{i=1}^{n}\|X_{i\cdot}-A_{i\cdot}\|_{2}^{2}+\gamma_{1}\sum_{l\in{\cal
E}_{1}}w_{l}\|A_{l_{1}\cdot}-A_{l_{2}\cdot}\|_{q}+\gamma_{2}\sum_{k\in{\cal
E}_{2}}u_{k}\|\mathbf{a}_{k_{1}}-\mathbf{a}_{k_{2}}\|_{q},$ (1)

where the weights $w_{l}\geq 0$ and $u_{k}\geq 0$, and $\|\cdot\|_{q}$ is the
$L_{q}$-norm of a vector with $q\in\\{1,2,\infty\\}$. Chi et al., (2017) and
Wang et al., (2018) suggested
$w_{l}=\iota_{l_{1},l_{2}}^{m_{1}}\exp(-\phi\|X_{l_{1}\cdot}-X_{l_{2}\cdot}\|_{2}^{2})$,
where $\iota_{l_{1},l_{2}}^{m_{1}}$ is 1 if observation $l_{2}$ is among
$l_{1}$’s $m_{1}$ nearest neighbors or vice verse, and 0 otherwise. Similarly,
we suggest using
$u_{k}=\iota_{k_{1},k_{2}}^{m_{2}}\exp(-\phi\|\mathbf{x}_{k_{1}}-\mathbf{x}_{k_{2}}\|_{2}^{2})$.

By incorporating the fused-lasso penalty in the second and third terms of (1),
the above formulation encourages that some of the rows or columns of the
solution $\widehat{{\bf A}}$ are identical. If
$\widehat{A}_{l_{1}\cdot}=\widehat{A}_{l_{2}\cdot}$, then observation $l_{1}$
and observation $l_{2}$ will be assigned to the same observation (row)
cluster; if $\widehat{\mathbf{a}}_{k_{1}}=\widehat{\mathbf{a}}_{k_{2}}$, then
feature $k_{1}$ and feature $k_{2}$ will be assigned to the same feature
(column) cluster. The tuning parameter $\gamma_{1}$ in (1) controls the number
of unique rows of $\widehat{{\bf A}}$, that is, the number of estimated
observation clusters, while the tuning parameter $\gamma_{2}$ controls the
number of estimated feature clusters. When $\gamma_{1}=\gamma_{2}=0$,
$\widehat{{\bf A}}={{\bf X}}$, and thus, each combination of feature and
observation by itself is a bicluster. As $\gamma_{1}$ and $\gamma_{2}$
increases, some of the rows and columns of $\widehat{{\bf A}}$ become
identical, respectively, which demonstrates a fusing process. For sufficiently
large $\gamma_{1}$ and $\gamma_{2}$, all rows or columns of $\widehat{{\bf
A}}$ will be identical, implying that all observations or features are
estimated to form a single bicluster. With such a strictly convex objective
function in (1), the solution $\widehat{{\bf A}}$ from convex biclustering is
unique for each given pair of $\gamma_{1}$ and $\gamma_{2}$.

###  2.2 bi-ADMM: a Biclustering Algorithm for the General Model

In order to implement the ADMM algorithm, we rewrite the above problem into
the following constrained optimization problem:

$\displaystyle\min_{\mathbf{A}\in\mathbb{R}^{p\times n}}$
$\displaystyle{1\over
2}\sum_{i=1}^{n}\|X_{i\cdot}-A_{i\cdot}\|_{2}^{2}+\gamma_{1}\sum_{l\in{\cal
E}_{1}}w_{l}\|\mathbf{v}_{l}\|_{q}+\gamma_{2}\sum_{k\in{\cal
E}_{2}}u_{k}\|\mathbf{z}_{k}\|_{q}$ (2) $\displaystyle{\rm s.t.}$
$\displaystyle A_{l_{1}\cdot}-A_{l_{2}\cdot}-\mathbf{v}_{l}={\bf 0},\ \forall
l\in{\cal E}_{1}$
$\displaystyle\mathbf{a}_{k_{1}}-\mathbf{a}_{k_{2}}-\mathbf{z}_{k}={\bf 0},\
\forall k\in{\cal E}_{2},$

which is also equivalent to an augmented Lagrangian problem given by

$\displaystyle\mathcal{L}_{\nu_{1},\nu_{2}}(\mathbf{A},\mathbf{V},\mathbf{Z},{\mbox{\boldmath$\Lambda$}}_{1},{\mbox{\boldmath$\Lambda$}}_{2})$
$\displaystyle=$ $\displaystyle{1\over
2}\sum_{i=1}^{n}\|X_{i\cdot}-A_{i\cdot}\|_{2}^{2}+\gamma_{1}\sum_{l\in{\cal
E}_{1}}w_{l}\|\mathbf{v}_{l}\|_{q}+\gamma_{2}\sum_{k\in{\cal
E}_{2}}u_{k}\|\mathbf{z}_{k}\|_{q}$ (3) $\displaystyle+\sum_{l\in{\cal
E}_{1}}\langle{\mbox{\boldmath$\lambda$}}_{1l},\mathbf{v}_{l}-A_{l_{1}\cdot}+A_{l_{2}\cdot}\rangle+{\nu_{1}\over
2}\sum_{l\in{\cal
E}_{1}}\|\mathbf{v}_{l}-A_{l_{1}\cdot}+A_{l_{2}\cdot}\|_{2}^{2}$
$\displaystyle+\sum_{k\in{\cal
E}_{2}}\langle{\mbox{\boldmath$\lambda$}}_{2k},\mathbf{z}_{k}-\mathbf{a}_{k_{1}}+\mathbf{a}_{k_{2}}\rangle+{\nu_{2}\over
2}\sum_{k\in{\cal
E}_{2}}\|\mathbf{z}_{k}-\mathbf{a}_{k_{1}}+\mathbf{a}_{k_{2}}\|_{2}^{2},$

where $\nu_{1}$ and $\nu_{2}$ are nonnegative constants,
$\mathbf{V}=(\mathbf{v}_{1},\ldots,\mathbf{v}_{|{{\cal E}}_{1}|})$,
$\mathbf{Z}=(\mathbf{z}_{1},\ldots,\mathbf{z}_{|{{\cal E}}_{2}|})$,
${\mbox{\boldmath$\Lambda$}}_{1}=({\mbox{\boldmath$\lambda$}}_{11},\ldots,{\mbox{\boldmath$\lambda$}}_{1|{{\cal
E}}_{1}|})$ and
${\mbox{\boldmath$\Lambda$}}_{2}=({\mbox{\boldmath$\lambda$}}_{21},\ldots,{\mbox{\boldmath$\lambda$}}_{2|{{\cal
E}}_{2}|})$. This optimization problem is challenging in dealing with fusing
both observation-level and feature-level vectors in the objective function
simultaneously.

The bi-ADMM minimizes the augmented Lagrangian problem by alternatively
solving one block of variables at a time. Specifically, bi-ADMM solves

$\displaystyle\mathbf{A}^{m+1}$ $\displaystyle=$ $\displaystyle\mathop{\rm
argmin}_{\mathbf{A}}{\cal
L}_{\nu_{1},\nu_{2}}(\mathbf{A},\mathbf{V}^{m},\mathbf{Z}^{m},{\mbox{\boldmath$\Lambda$}}_{1}^{m},{\mbox{\boldmath$\Lambda$}}_{2}^{m}),$
$\displaystyle\mathbf{V}^{m+1}$ $\displaystyle=$ $\displaystyle\mathop{\rm
argmin}_{\mathbf{V}}{\cal
L}_{\nu_{1},\nu_{2}}(\mathbf{A}^{m+1},\mathbf{V},\mathbf{Z}^{m},{\mbox{\boldmath$\Lambda$}}_{1}^{m},{\mbox{\boldmath$\Lambda$}}_{2}^{m}),$
$\displaystyle\mathbf{Z}^{m+1}$ $\displaystyle=$ $\displaystyle\mathop{\rm
argmin}_{\mathbf{Z}}{\cal
L}_{\nu_{1},\nu_{2}}(\mathbf{A}^{m+1},\mathbf{V}^{m+1},\mathbf{Z},{\mbox{\boldmath$\Lambda$}}_{1}^{m},{\mbox{\boldmath$\Lambda$}}_{2}^{m}),$
(4) $\displaystyle{\mbox{\boldmath$\lambda$}}_{1l}^{m+1}$ $\displaystyle=$
$\displaystyle{\mbox{\boldmath$\lambda$}}_{1l}^{m}+\nu_{1}(\mathbf{v}_{l}^{m+1}-A_{l_{1}\cdot}^{m+1}+A_{l_{2}\cdot}^{m+1}),\
l\in{\cal E}_{1},$ $\displaystyle{\mbox{\boldmath$\lambda$}}_{2k}^{m+1}$
$\displaystyle=$
$\displaystyle{\mbox{\boldmath$\lambda$}}_{2k}^{m}+\nu_{2}(\mathbf{z}_{k}^{m+1}-\mathbf{a}_{k_{1}}^{m+1}+\mathbf{a}_{k_{2}}^{m+1}),\
k\in{\cal E}_{2}.$

Next, we develop the detailed updating implementations for
$\mathbf{A},\mathbf{V},\mathbf{Z},{\mbox{\boldmath$\Lambda$}}_{1}$ and
${\mbox{\boldmath$\Lambda$}}_{2}$ in three steps. A summary of the bi-ADMM
algorithm is shown in Algorithm 1.

Step 1: update $\mathbf{A}$. To update $\mathbf{A}$, we need to minimize

$\displaystyle f(\mathbf{A})$ $\displaystyle=$ $\displaystyle{1\over
2}\sum_{i=1}^{n}\|X_{i\cdot}-A_{i\cdot}\|_{2}^{2}+{\nu_{1}\over
2}\sum_{l\in{\cal
E}_{1}}\|\widetilde{\mathbf{v}}_{l}-A_{l_{1}\cdot}+A_{l_{2}\cdot}\|_{2}^{2}+{\nu_{2}\over
2}\sum_{k\in{\cal
E}_{2}}\|\widetilde{\mathbf{z}}_{k}-\mathbf{a}_{k_{1}}+\mathbf{a}_{k_{2}}\|_{2}^{2},$
(5)

where
$\widetilde{\mathbf{v}}_{1}=\mathbf{v}_{l}+{1\over\nu_{1}}{\mbox{\boldmath$\lambda$}}_{1l}$
and
$\widetilde{\mathbf{z}}_{k}=\mathbf{z}_{k}+{1\over\nu_{2}}{\mbox{\boldmath$\lambda$}}_{2k}$.

The optimization problem is challenging because it involves both rows and
columns of a single matrix. To tackle this difficulty, the following key lemma
associating $(\ref{eqn:bi_update_A})$ with a Sylvester Equations is proposed.
The proof of Lemma 1 is provided in Appendix.

######  Lemma 1.

Let ${\bf I}_{n}$ be an $n\times n$ identity matrix, ${\bf 1}_{n}$ be an
$n$-dimensional vector with each component being 1, $\mathbf{e}_{i}$ be an
$n$-dimensional vector with all components being 0, except the $i$-th
component being 1, and $\mathbf{e}_{j}^{*}$ is a $p$-dimensional vector with
the $j$-th element as 1 and 0 otherwise. In addition, define

$\displaystyle\mathbf{M}$ $\displaystyle=$
$\displaystyle\mathbf{I}_{n}+\nu_{1}\sum_{l\in{\cal
E}_{1}}(\mathbf{e}_{l_{1}}-\mathbf{e}_{l_{2}})(\mathbf{e}_{l_{1}}-\mathbf{e}_{l_{2}})^{\rm
T}$ $\displaystyle\mathbf{N}$ $\displaystyle=$
$\displaystyle\nu_{2}\sum_{k\in{\cal
E}_{2}}(\mathbf{e}_{k_{1}}^{*}-\mathbf{e}_{k_{2}}^{*})(\mathbf{e}_{k_{1}}^{*}-\mathbf{e}_{k_{2}}^{*})^{\rm
T}$ $\displaystyle\mathbf{G}$ $\displaystyle=$
$\displaystyle\mathbf{X}+\sum_{l\in{\cal
E}_{1}}(\mathbf{e}_{l_{1}}-\mathbf{e}_{l_{2}})({\mbox{\boldmath$\lambda$}}_{1l}+\nu_{1}\mathbf{v}_{l})^{\rm
T}+\sum_{k\in{\cal
E}_{2}}({\mbox{\boldmath$\lambda$}}_{2k}+\nu_{2}\mathbf{z}_{k})(\mathbf{e}_{k_{1}}^{*}-\mathbf{e}_{k_{2}}^{*})^{\rm
T}.$

Then, minimizing $(\ref{eqn:bi_update_A})$ is equivalent to solving the
following Sylvester Equation

$\displaystyle\mathbf{M}\mathbf{A}+\mathbf{A}\mathbf{N}=\mathbf{G}.$ (6)

If the edge sets ${\cal E}_{1}$ and ${\cal E}_{2}$ contain all possible edges,
it is straightforward to verify

$\displaystyle\sum_{l\in{\cal
E}_{1}}(\mathbf{e}_{l_{1}}-\mathbf{e}_{l_{2}})(\mathbf{e}_{l_{1}}-\mathbf{e}_{l_{2}})^{\rm
T}$ $\displaystyle=$ $\displaystyle n\mathbf{I}_{n}-{\bf 1}_{n}{\bf
1}_{n}^{\rm T}$ $\displaystyle\sum_{k\in{\cal
E}_{2}}(\mathbf{e}_{k_{1}}^{*}-\mathbf{e}_{k_{2}}^{*})(\mathbf{e}_{k_{1}}^{*}-\mathbf{e}_{k_{2}}^{*})^{\rm
T}$ $\displaystyle=$ $\displaystyle p\mathbf{I}_{p}-{\bf 1}_{p}{\bf
1}_{p}^{\rm T}.$

Then, $\mathbf{M}=(1+n\nu_{1})\mathbf{I}_{n}-\nu_{1}{\bf 1}_{n}{\bf
1}_{n}^{\rm T}$ and $\mathbf{N}=p\nu_{2}\mathbf{I}_{p}-\nu_{2}{\bf 1}_{p}{\bf
1}_{p}^{\rm T}$.

The type of equation (6) is called Sylvester Equation, which is very important
in control theory and many other branches of engineering. Its theoretical
solution is based on eigenvector and eigenvalue decomposition (Jameson,,
1968), but it is computationally expensive. The Bartels-Stewart algorithm
(Bartels and Stewart,, 1972) is the standard numerical solution that
transforms the Sylvester Equation into a triangular system with the Schur
decomposition and then solves it with forward or backward substitutions. In
this manuscript and its accompanied R package, we implement a modified
Bartels-Stewart algorithm proposed by Sorensen et al., (2003), which is more
efficient.

Step 2: Update $\mathbf{V}$ and $\mathbf{Z}$. For any $\sigma>0$ and norm
$\Omega(\cdot)$, we define a proximal map,

$\displaystyle\textrm{prox}_{\sigma\Omega}({\bf u})=\mathop{\rm argmin}_{{\bf
v}}\left[\sigma\Omega({\bf v})+{1\over 2}\|{\bf u}-{\bf v}\|_{2}^{2}\right].$

We refer the readers to Table 1 in Chi and Lange, (2015) for the solutions to
the proximal map of $L_{q}$-norm for $q=1,2$ and $\infty$. It is clear that
the vectors ${\mathbf{v}}_{l}$ and ${\mathbf{z}}_{k}$ are separable in the
objective function (3), and thus ${\mathbf{v}}_{l}$ and ${\mathbf{z}}_{k}$ can
be solved via proximal maps:

$\displaystyle\mathbf{v}_{l}$ $\displaystyle=$ $\displaystyle\mathop{\rm
argmin}_{\mathbf{v}_{l}}{1\over
2}\|\mathbf{v}_{l}-(A_{l_{1}\cdot}-A_{l_{2}\cdot}-\nu_{1}^{-1}{\mbox{\boldmath$\lambda$}}_{1l})\|_{2}^{2}+{\gamma_{1}w_{l}\over\nu_{1}}\|\mathbf{v}_{l}\|_{q}$
$\displaystyle=$
$\displaystyle\textrm{prox}_{\sigma_{1l}\|\cdot\|_{q}}(A_{l_{1}\cdot}-A_{l_{2}\cdot}-\nu_{1}^{-1}{\mbox{\boldmath$\lambda$}}_{1l})$
$\displaystyle\mathbf{z}_{k}$ $\displaystyle=$ $\displaystyle\mathop{\rm
argmin}_{\mathbf{z}_{k}}{1\over
2}\|\mathbf{z}_{k}-(\mathbf{a}_{k_{1}}-\mathbf{a}_{k_{2}}-\nu_{2}^{-1}{\mbox{\boldmath$\lambda$}}_{2k})\|_{2}^{2}+{\gamma_{2}u_{k}\over\nu_{2}}\|\mathbf{z}_{k}\|_{q}$
$\displaystyle=$
$\displaystyle\textrm{prox}_{\sigma_{2k}\|\cdot\|_{q}}(\mathbf{a}_{k_{1}}-\mathbf{a}_{k_{2}}-\nu_{2}^{-1}{\mbox{\boldmath$\lambda$}}_{2k}),$

where $\sigma_{1l}=\gamma_{1}w_{l}/\nu_{1}$ and
$\sigma_{2k}=\gamma_{2}u_{k}/\nu_{2}$.

Step 3: Update ${\mbox{\boldmath$\Lambda$}}_{1}$ and
${\mbox{\boldmath$\Lambda$}}_{2}$. Finally, ${\mbox{\boldmath$\lambda$}}_{1l}$
and ${\mbox{\boldmath$\lambda$}}_{2k}$ can be updated by
${\mbox{\boldmath$\lambda$}}_{1l}={\mbox{\boldmath$\lambda$}}_{1l}+\nu_{1}(\mathbf{v}_{l}-A_{l_{1}\cdot}+A_{l_{2}\cdot})$
and
${\mbox{\boldmath$\lambda$}}_{2k}={\mbox{\boldmath$\lambda$}}_{2k}+\nu_{2}(\mathbf{z}_{k}-\mathbf{a}_{k_{1}}+\mathbf{a}_{k_{2}})$.

Algorithm 1 bi-ADMM

  1. 1.

Initialize $\mathbf{V}^{0},\mathbf{Z}^{0},{\mbox{\boldmath$\Lambda$}}_{1}^{0}$
and ${\mbox{\boldmath$\Lambda$}}_{2}^{0}$. Calculate

$\displaystyle\mathbf{M}$ $\displaystyle=$
$\displaystyle\mathbf{I}_{n}+\nu_{1}\sum_{l\in{\cal
E}_{1}}(\mathbf{e}_{l_{1}}-\mathbf{e}_{l_{2}})(\mathbf{e}_{l_{1}}-\mathbf{e}_{l_{2}})^{\rm
T}$ $\displaystyle\mathbf{N}$ $\displaystyle=$
$\displaystyle\nu_{2}\sum_{k\in{\cal
E}_{2}}(\mathbf{e}_{k_{1}}^{*}-\mathbf{e}_{k_{2}}^{*})(\mathbf{e}_{k_{1}}^{*}-\mathbf{e}_{k_{2}}^{*})^{\rm
T}$

For $m=1,2,\ldots$

  2. 2.

Solve the Sylvester Equation
$\mathbf{M}\mathbf{A}+\mathbf{A}\mathbf{N}=\mathbf{G}^{m-1}$ to obtain
$\mathbf{A}^{m}$, where

$\displaystyle\mathbf{G}^{m-1}=\mathbf{X}+\sum_{l\in{\cal
E}_{1}}(\mathbf{e}_{l_{1}}-\mathbf{e}_{l_{2}})({\mbox{\boldmath$\lambda$}}_{1l}^{m-1}+\nu_{1}\mathbf{v}_{l}^{m-1})^{\rm
T}+\sum_{k\in{\cal
E}_{2}}({\mbox{\boldmath$\lambda$}}_{2k}^{m-1}+\nu_{2}\mathbf{z}_{k}^{m-1})(\mathbf{e}_{k_{1}}^{*}-\mathbf{e}_{k_{2}}^{*})^{\rm
T}.$

  3. 3.

For $l\in{\cal E}_{1}$, do

$\displaystyle\mathbf{v}_{l}^{m}=\textrm{prox}_{\sigma_{1l}\|\cdot\|_{q}}(A_{l_{1}\cdot}^{m}-A_{l_{2}\cdot}^{m}-\nu_{1}^{-1}{\mbox{\boldmath$\lambda$}}_{1l}^{m-1}).$

  4. 4.

For $k\in{\cal E}_{2}$, do

$\displaystyle\mathbf{z}_{l}^{m}=\textrm{prox}_{\sigma_{2k}\|\cdot\|_{q}}(\mathbf{a}_{k_{1}}^{m}-\mathbf{a}_{k_{2}}^{m}-\nu_{2}^{-1}{\mbox{\boldmath$\lambda$}}_{2k}^{m-1}).$

  5. 5.

For $l\in{\cal E}_{1}$ and $k\in{\cal E}_{2}$, do

$\displaystyle{\mbox{\boldmath$\lambda$}}_{1l}^{m}$ $\displaystyle=$
$\displaystyle{\mbox{\boldmath$\lambda$}}_{1l}^{m-1}+\nu_{1}(\mathbf{v}_{l}^{m}-A_{l_{1}\cdot}^{m}+A_{l_{2}\cdot}^{m})$
$\displaystyle{\mbox{\boldmath$\lambda$}}_{2k}^{m}$ $\displaystyle=$
$\displaystyle{\mbox{\boldmath$\lambda$}}_{2k}^{m-1}+\nu_{2}(\mathbf{z}_{k}^{m}-\mathbf{a}_{k_{1}}^{m}+\mathbf{a}_{k_{2}}^{m})$

  6. 6.

Repeat Steps 2-5 until convergence.

###  2.3 biC-ADMM: Biclustering Algorithm for Data with Compositional
Constraints

For microbiome data, due to its data-generating mechanism, the counts of all
taxa need to be normalized into the relative abundances such that their
summation for each sample is 1. This type of data is also called compositional
data in mathematics. Our newly proposed clustering algorithm, biC-ADMM, aims
to incorporate the compositional constraints into the optimization. To be
specific, the data from each sample satisfy $X_{i\cdot}^{\rm T}{\bf
1}_{p}=1,i=1,\ldots,n$, and then it is natural to require the estimate
$\mathbf{A}$ to satisfy the same constraints, i.e., $A_{i\cdot}^{\rm T}{\bf
1}_{p}=1$.

Because the bi-ADMM solves an optimization problem with constraints, it can be
easily extended to the biC-ADMM algorithm by incorporating the compositional
constraints. With the additional compositional constraints into (2), the
augmented Lagrangian problem for compositional data is given by

$\displaystyle\mathcal{L}_{\nu_{1},\nu_{2},\nu_{3}}(\mathbf{A},\mathbf{V},\mathbf{Z},{\mbox{\boldmath$\Lambda$}}_{1},{\mbox{\boldmath$\Lambda$}}_{2},{\mbox{\boldmath$\lambda$}}_{3})$
$\displaystyle=$ $\displaystyle{1\over
2}\sum_{i=1}^{n}\|X_{i\cdot}-A_{i\cdot}\|_{2}^{2}+\gamma_{1}\sum_{l\in{\cal
E}_{1}}w_{l}\|\mathbf{v}_{l}\|_{q}+\gamma_{2}\sum_{k\in{\cal
E}_{2}}u_{k}\|\mathbf{z}_{k}\|_{q}$ (7) $\displaystyle+\sum_{l\in{\cal
E}_{1}}\langle{\mbox{\boldmath$\lambda$}}_{1l},\mathbf{v}_{l}-A_{l_{1}\cdot}+A_{l_{2}\cdot}\rangle+{\nu_{1}\over
2}\sum_{l\in{\cal
E}_{1}}\|\mathbf{v}_{l}-A_{l_{1}\cdot}+A_{l_{2}\cdot}\|_{2}^{2}$
$\displaystyle+\sum_{k\in{\cal
E}_{2}}\langle{\mbox{\boldmath$\lambda$}}_{2k},\mathbf{z}_{k}-\mathbf{a}_{k_{1}}+\mathbf{a}_{k_{2}}\rangle+{\nu_{2}\over
2}\sum_{k\in{\cal
E}_{2}}\|\mathbf{z}_{k}-\mathbf{a}_{k_{1}}+\mathbf{a}_{k_{2}}\|_{2}^{2}$
$\displaystyle+\langle{\mbox{\boldmath$\lambda$}}_{3},{\bf
1}_{n}-\mathbf{A}{\bf 1}_{p}\rangle+{\nu_{3}\over 2}\|{\bf
1}_{n}-\mathbf{A}{\bf 1}_{p}\|_{2}^{2},$

where $\nu_{3}$ is a nonnegative constant and
${\mbox{\boldmath$\lambda$}}_{3}=(\lambda_{31},\ldots,\lambda_{3n})^{\rm T}$.

As a variant of the bi-ADMM for biclustering compositional data, the biC-ADMM
algorithm is similar with the bi-ADMM detailed in Section 2.2. The key step
for updating $\mathbf{A}$ of the biC-ADMM is summarized in the following Lemma
2, which minimizes:

$\displaystyle f(\mathbf{a})={1\over
2}\|\mathbf{x}-\mathbf{a}\|_{2}^{2}+{\nu_{1}\over
2}\|\mathbf{B}\mathbf{P}\mathbf{a}-\widetilde{\mathbf{v}}\|_{2}^{2}+{\nu_{2}\over
2}\|\mathbf{C}\mathbf{a}-\widetilde{\mathbf{z}}\|_{2}^{2}+{\nu_{3}\over
2}\|\mathbf{D}\mathbf{a}-{\bf s}\|_{2}^{2},$ (8)

with $\mathbf{x}=\textrm{vec}(\mathbf{X})$,
$\mathbf{a}=\textrm{vec}(\mathbf{A})$, the vectorization of the matrices
$\mathbf{X}$ and $\mathbf{A}$, and

$\displaystyle\mathbf{B}^{\rm T}=\left(\mathbf{B}^{\rm
T}_{1},\ldots,\mathbf{B}^{\rm T}_{|{\cal
E}_{1}|}\right),\quad\widetilde{\mathbf{v}}^{\rm
T}=\left(\widetilde{\mathbf{v}}_{1}^{\rm
T},\ldots,\widetilde{\mathbf{v}}_{|{\cal E}_{1}|}^{\rm T}\right),$
$\displaystyle\mathbf{C}^{\rm T}=\left(\mathbf{C}^{\rm
T}_{1},\ldots,\mathbf{C}^{\rm T}_{|{\cal
E}_{2}|}\right),\quad\widetilde{\mathbf{z}}^{\rm
T}=\left(\widetilde{\mathbf{z}}_{1}^{\rm
T},\ldots,\widetilde{\mathbf{z}}_{|{\cal E}_{2}|}^{\rm T}\right),$
$\displaystyle\mathbf{D}^{\rm T}=\left(\mathbf{D}^{\rm
T}_{1},\ldots,\mathbf{D}^{\rm T}_{n}\right),\quad{\bf
s}=\left(s_{1},\ldots,s_{n}\right)^{\rm T}={\bf
1}_{n}+{\mbox{\boldmath$\lambda$}}_{3}/\nu_{3},$

where $\mathbf{B}_{l}=(\mathbf{e}_{l_{1}}-\mathbf{e}_{l2})^{\rm
T}\otimes\mathbf{I}_{p},\mathbf{C}_{k}=(\mathbf{e}_{k_{1}}^{*}-\mathbf{e}_{k2}^{*})^{\rm
T}\otimes\mathbf{I}_{n},\mathbf{D}_{i}={\bf 1}_{p}^{\rm
T}\otimes\mathbf{e}_{i}^{\rm T}$, and $\mathbf{P}$ is a permutation matrix
such that $\textrm{vec}(\mathbf{A}^{\rm
T})=\mathbf{P}\textrm{vec}(\mathbf{A}),l=1,\ldots,|{\cal
E}_{1}|,k=1,\ldots,|{\cal E}_{2}|,i=1,\ldots,n$, respectively.

######  Lemma 2.

Define

$\displaystyle\mathbf{M}$ $\displaystyle=$
$\displaystyle\mathbf{I}_{n}+\nu_{1}\sum_{l\in{\cal
E}_{1}}(\mathbf{e}_{l_{1}}-\mathbf{e}_{l_{2}})(\mathbf{e}_{l_{1}}-\mathbf{e}_{l_{2}})^{\rm
T}$ $\displaystyle\mathbf{N}$ $\displaystyle=$
$\displaystyle\nu_{2}\sum_{k\in{\cal
E}_{2}}(\mathbf{e}_{k_{1}}^{*}-\mathbf{e}_{k_{2}}^{*})(\mathbf{e}_{k_{1}}^{*}-\mathbf{e}_{k_{2}}^{*})^{\rm
T}+\nu_{3}{\bf 1}_{p}{\bf 1}_{p}^{\rm T}$ $\displaystyle\mathbf{G}$
$\displaystyle=$ $\displaystyle\mathbf{X}+\sum_{l\in{\cal
E}_{1}}(\mathbf{e}_{l_{1}}-\mathbf{e}_{l_{2}})({\mbox{\boldmath$\lambda$}}_{1l}+\nu_{1}\mathbf{v}_{l})^{\rm
T}+\sum_{k\in{\cal
E}_{2}}({\mbox{\boldmath$\lambda$}}_{2k}+\nu_{2}\mathbf{z}_{k})(\mathbf{e}_{k_{1}}^{*}-\mathbf{e}_{k_{2}}^{*})^{\rm
T}+\nu_{3}\mathbf{s}{\bf 1}_{p}^{\rm T}.$

Then, minimizing $(\ref{eqn:biC_update_A})$ is equivalent to solving the
following Sylvester Equation

$\displaystyle\mathbf{M}\mathbf{A}+\mathbf{A}\mathbf{N}=\mathbf{G}.$ (9)

Proof of Lemma 2 is provided in the Appendix, and the biC-ADMM algorithm is
summarized in Algorithm 2.

Algorithm 2 biC-ADMM

  1. 1.

Initialize
$\mathbf{V}^{0},\mathbf{Z}^{0},{\mbox{\boldmath$\Lambda$}}_{1}^{0},{\mbox{\boldmath$\Lambda$}}_{2}^{0}$
and ${\mbox{\boldmath$\lambda$}}_{3}^{0}$. Calculate

$\displaystyle\mathbf{M}$ $\displaystyle=$
$\displaystyle\mathbf{I}_{n}+\nu_{1}\sum_{l\in{\cal
E}_{1}}(\mathbf{e}_{l_{1}}-\mathbf{e}_{l_{2}})(\mathbf{e}_{l_{1}}-\mathbf{e}_{l_{2}})^{\rm
T}$ $\displaystyle\mathbf{N}$ $\displaystyle=$
$\displaystyle\nu_{2}\sum_{k\in{\cal
E}_{2}}(\mathbf{e}_{k_{1}}^{*}-\mathbf{e}_{k_{2}}^{*})(\mathbf{e}_{k_{1}}^{*}-\mathbf{e}_{k_{2}}^{*})^{\rm
T}+\nu_{3}{\bf 1}_{p}{\bf 1}_{p}^{\rm T}$

For $m=1,2,\ldots$

  2. 2.

Solve the Sylvester Equation
$\mathbf{M}\mathbf{A}+\mathbf{A}\mathbf{N}=\mathbf{G}^{m-1}$ to obtain
$\mathbf{A}^{m}$, where

$\displaystyle\mathbf{G}^{m-1}=\mathbf{X}+\sum_{l\in{\cal
E}_{1}}(\mathbf{e}_{l_{1}}-\mathbf{e}_{l_{2}})({\mbox{\boldmath$\lambda$}}_{1l}^{m-1}+\nu_{1}\mathbf{v}_{l}^{m-1})^{\rm
T}+\sum_{k\in{\cal
E}_{2}}({\mbox{\boldmath$\lambda$}}_{2k}^{m-1}+\nu_{2}\mathbf{z}_{k}^{m-1})(\mathbf{e}_{k_{1}}^{*}-\mathbf{e}_{k_{2}}^{*})^{\rm
T}+\nu_{3}\mathbf{s}^{m-1}{\bf 1}_{p}^{\rm T},$

where $\mathbf{s}^{m-1}={\bf
1}_{n}+{\mbox{\boldmath$\lambda$}}_{3}^{m-1}/\nu_{3}$.

  3. 3.

For $l\in{\cal E}_{1}$, do

$\displaystyle\mathbf{v}_{l}^{m}=\textrm{prox}_{\sigma_{1l}\|\cdot\|_{q}}(A_{l_{1}\cdot}^{m}-A_{l_{2}\cdot}^{m}-\nu_{1}^{-1}{\mbox{\boldmath$\lambda$}}_{1l}^{m-1}).$

  4. 4.

For $k\in{\cal E}_{2}$, do

$\displaystyle\mathbf{z}_{l}^{m}=\textrm{prox}_{\sigma_{2k}\|\cdot\|_{q}}(\mathbf{a}_{k_{1}}^{m}-\mathbf{a}_{k_{2}}^{m}-\nu_{2}^{-1}{\mbox{\boldmath$\lambda$}}_{2k}^{m-1}).$

  5. 5.

For $l\in{\cal E}_{1}$ and $k\in{\cal E}_{2}$, do

$\displaystyle{\mbox{\boldmath$\lambda$}}_{1l}^{m}$ $\displaystyle=$
$\displaystyle{\mbox{\boldmath$\lambda$}}_{1l}^{m-1}+\nu_{1}(\mathbf{v}_{l}^{m}-A_{l_{1}\cdot}^{m}+A_{l_{2}\cdot}^{m})$
$\displaystyle{\mbox{\boldmath$\lambda$}}_{2k}^{m}$ $\displaystyle=$
$\displaystyle{\mbox{\boldmath$\lambda$}}_{2k}^{m-1}+\nu_{2}(\mathbf{z}_{k}^{m}-\mathbf{a}_{k_{1}}^{m}+\mathbf{a}_{k_{2}}^{m})$
$\displaystyle{\mbox{\boldmath$\lambda$}}_{3}^{m}$ $\displaystyle=$
$\displaystyle{\mbox{\boldmath$\lambda$}}_{3}^{m-1}+\nu_{3}({\bf
1}_{n}-\mathbf{A}^{m}{\bf 1}_{p})$

  6. 6.

Repeat Steps 2-5 until convergence.

##  3 Implementation of bi-ADMM and biC-ADMM in Practical Settings

In this section, we discuss a few practical issues for implementing the
proposed algorithms, including algorithmic convergence and selection of tuning
parameters. The weights $w_{l}$ and $u_{k}$ are suggested in Section 2.1.

###  3.1 Algorithmic Convergence

The proposed convex biclustering algorithms are derived by applying standard
techniques of the ADMM method. For convex clustering problems, Chi and Lange,
(2015) and Wang et al., (2018) as well as the references therein provided
sufficient conditions for the convergence of their proposed methods. Chi et
al., (2017) discussed the convergence of their algorithm COBRA for convex
biclustering problem.

The convergence of our bi-ADMM and biC-ADMM algorithms follows similar
arguments. Note that the only difference between the objective function in (1)
and its counterpart in Chi and Lange, (2015) is an additional fused-lasso
penalty term for feature-level vectors, which is strongly convex. According to
Chi and Lange, (2015), one can show that, under mild regularization
conditions, the convergence of the bi-ADMM and biC-ADMM algrithms is
guaranteed for any $\nu_{j}>0,j=1,2,3$.

###  3.2 Selection of Tuning Parameters

This subsection discusses the methods for selecting tuning parameters
$\gamma_{1}$ and $\gamma_{2}$. Recall that $\gamma_{1}$ controls the number of
estimated observation-level clusters and $\gamma_{2}$ controls the number of
estimated feature-level clusters. The usage of two separate tuning parameters
to control the numbers of row and column clusters separately provide great
flexibility in applications.

Generally speaking, a two-dimensional grid search can provide an optimal pair
of tuning parameters under some criterion if the computational resource is
sufficient. Sometimes, if the data matrix is in a large scale, we may adopt a
similar strategy suggested in Chi et al., (2017) by using a single tuning
parameter $\gamma$ for properly rescaled penalty terms. Using a single tuning
parameter can reduce the computational burden, but row and columns have to be
clustered in a proportional manner.

To be specific, we rewrite the convex biclustering problem as the following
minimization problem with a single tuning parameter $\gamma$:

$\displaystyle\min_{\mathbf{A}\in\mathbb{R}^{p\times n}}$
$\displaystyle{1\over
2}\sum_{i=1}^{n}\|X_{i\cdot}-A_{i\cdot}\|_{2}^{2}+\gamma\left\\{\sum_{l\in{\cal
E}_{1}}w_{l}\|A_{l_{1}\cdot}-A_{l_{2}\cdot}\|_{q}+\sum_{k\in{\cal
E}_{2}}u_{k}\|\mathbf{a}_{k_{1}}-\mathbf{a}_{k_{2}}\|_{q}\right\\}.$ (10)

The key to the validity of this formulation is that the two penalty terms
$\sum_{l\in{\cal E}_{1}}w_{l}\|A_{l_{1}\cdot}-A_{l_{2}\cdot}\|_{q}$ and
$\sum_{k\in{\cal E}_{2}}u_{k}\|\mathbf{a}_{k_{1}}-\mathbf{a}_{k_{2}}\|_{q}$
should be on the same scale. Note that row vectors of $\mathbf{A}$ are in
$\mathbb{R}^{p}$ while column vectors of $\mathbf{A}$ are in $\mathbb{R}^{n}$.
Thus, we choose row weights $w_{l}$ to sum to $1/\sqrt{p}$ and the column
weights $u_{k}$ to sum to $1/\sqrt{n}$.

Next, we discuss two criteria for tuning parameter selection in a data-driven
manner. One is spiritually similar with cross-validation and the other is
based on stability selection. Chi et al., (2017) proposed a hold-out
validation method for convex biclustering by randomly selecting a hold-out of
elements in the data matrix and assessing the quality of predicting the hold-
out set with an estimated model based on the rest elements. Alternatively,
Wang et al., (2018) proposed to use stability selection in Fang and Wang,
(2012) to tune two parameters for sparse convex clustering. For convex
biclustering, we can apply stability section in a similar way to tune
$\gamma_{1}$ and $\gamma_{2}$. To be specific, for any given $\gamma_{1}$ and
$\gamma_{2}$, two biclustering results can be produced via (1) based on two
sets of bootstrapped samples, and then the stability measurement (Fang and
Wang,, 2012) can be computed to measure the agreement between the two
biclustering results. A repetition of this procedure for $50$ times is
required to enhance the robustness of the stability selection method. Finally,
the optimal parameters are selected as the one achieving maximum averaged
stability value.

##  4 Simulation Studies

In this section, we conduct simulation studies to thoroughly evaluate the
performance of the proposed bi-ADMM under general biclustering setups using
$L_{2}$\- and $L_{1}$-norms, and the biC-ADMM using $L_{2}$-norm under
scenarios with compositional constrains motivated by a real microbiome dataset
Livanos et al., (2016), and compare them with the competing method COBRA. We
do not illustrate the proposed algorithms with $L_{\infty}$-norm because extra
computation burden is needed to solve the proximal map with simplex algorithms
in Step 2 of Algorithms 1 and 2. For each setting, we run 100 repetitions. In
all simulation studies, the adjusted RAND index (ARI) (Hubert and Arabie,,
1985) is used to measure the agreement between the estimated clustering result
and the underlying true clustering assignment. The ARI ranges from 0 to 1 with
a higher value indicating better performance. Note that the underlying cluster
labels are known or pre-specified in simulation studies, and thus, it is
feasible to evaluate how well the candidate methods can perform if they are
tuned by maximizing the ARI.

To ensure fair comparisons with the COBRA, we adopt the formulation shown in
(10), using a single tuning parameter. For each repetition, an optimal tuning
parameter $\gamma$ is chosen by maximizing the ARI on a separate validation
data set for candidate algorithms.

###  4.1 General Setups

To generate datasets with checkerboard cluster structures, we consider various
sizes of datasets with $n$ observations (row) and $p$ features (column), where
$(n,p)$ are set as $(50,40),(100,80),(200,160)$, respectively. The number of
row and column clusters varies as well. Given the dataset size, $K=4,8,12$, or
$16$ row clusters and $R=4,8$, or $16$ column clusters are generated
accordingly. The corresponding total number of biclusters are $M\doteq K\times
R=16,32,64,96,128$, or $256$.

$X_{ij}$ is generated as follows. First, we assign cluster indexes to the
observations (rows) by sampling from a set $\\{1,...,K\\}$ uniformly, and the
cluster indexes are assigned to features (columns) following a similar
procedure. Thus, each $X_{ij}$ belongs to one of those $M$ biclusters. Second,
random samples for each bicluster are generated from a normal distribution,
$X_{ij}\text{ i.i.d. }\sim{\cal N}(\mu_{kr},\sigma^{2})$, i.e., samples from
the row cluster $k\in\\{1,\ldots,K\\}$ and column cluster
$r\in\\{1,\ldots,R\\}$ follow a normal distribution with mean $\mu_{kr}$ and
variance $\sigma^{2}$. $\mu_{kr}$ is chosen uniformly from a sequence
$\\{-10,-9,\ldots,9,10\\}$. We vary the variance $\sigma^{2}$ to change the
noise level, and the $\sigma$ is chosen as $2,4,6,8$, and $10$ in our
simulation settings, respectively. Both parameters $\nu_{1},\nu_{2}$ for the
bi-ADMM algorithm are set as 8.

Table 1: Simulation results for the bi-ADMM with $L_{2}$-norm, bi-ADMM with $L_{1}$-norm, and COBRA algorithms in terms of the ARI under various scenarios, respectively. The largest ARI(s) is bolded for each scenario. |  |  |  |  |  | bi-ADMM ($L_{2}$) | bi-ADMM ($L_{1}$) | COBRA
---|---|---|---|---|---|---|---|---
$n$ | $p$ | $K$ | $R$ | $M$ | $\sigma$ | mean | sd | mean | sd | mean | sd
50 | 40 | 4 | 4 | 16 | 2 | 1.00 | 0.03 | 1.00 | 0.02 | 0.99 | 0.05
|  |  |  |  | 4 | 0.96 | 0.09 | 0.96 | 0.09 | 0.88 | 0.15
|  |  |  |  | 6 | 0.85 | 0.16 | 0.83 | 0.18 | 0.66 | 0.19
|  |  |  |  | 8 | 0.65 | 0.21 | 0.63 | 0.19 | 0.42 | 0.17
|  | 8 | 4 | 32 | 2 | 0.99 | 0.02 | 0.99 | 0.03 | 0.97 | 0.06
|  |  |  |  | 4 | 0.92 | 0.08 | 0.92 | 0.08 | 0.82 | 0.12
|  |  |  |  | 6 | 0.76 | 0.12 | 0.77 | 0.11 | 0.57 | 0.15
|  |  |  |  | 8 | 0.54 | 0.16 | 0.55 | 0.14 | 0.33 | 0.12
100 | 80 | 4 | 4 | 16 | 2 | 1.00 | 0.02 | 1.00 | 0.03 | 1.00 | 0.02
|  |  |  |  | 4 | 0.97 | 0.06 | 0.94 | 0.10 | 0.91 | 0.11
|  |  |  |  | 6 | 0.86 | 0.15 | 0.75 | 0.19 | 0.84 | 0.20
|  |  |  |  | 8 | 0.67 | 0.18 | 0.54 | 0.19 | 0.65 | 0.20
|  | 8 | 4 | 32 | 4 | 0.91 | 0.09 | 0.91 | 0.09 | 0.93 | 0.08
|  |  |  |  | 6 | 0.73 | 0.12 | 0.71 | 0.20 | 0.75 | 0.13
|  |  |  |  | 8 | 0.53 | 0.12 | 0.51 | 0.14 | 0.53 | 0.14
|  | 8 | 8 | 64 | 4 | 0.98 | 0.05 | 0.96 | 0.06 | 0.98 | 0.05
|  |  |  |  | 6 | 0.82 | 0.17 | 0.78 | 0.13 | 0.81 | 0.13
|  |  |  |  | 8 | 0.45 | 0.14 | 0.48 | 0.13 | 0.30 | 0.19
200 | 160 | 8 | 8 | 64 | 4 | 1.00 | 0.02 | 0.99 | 0.03 | 1.00 | 0.02
|  |  |  |  | 6 | 0.95 | 0.07 | 0.92 | 0.09 | 0.95 | 0.07
|  |  |  |  | 8 | 0.80 | 0.13 | 0.71 | 0.15 | 0.77 | 0.14
|  |  |  |  | 10 | 0.51 | 0.17 | 0.42 | 0.16 | 0.46 | 0.16
|  | 12 | 8 | 256 | 4 | 1.00 | 0.02 | 1.00 | 0.02 | 1.00 | 0.01
|  |  |  |  | 6 | 0.95 | 0.05 | 0.94 | 0.06 | 0.97 | 0.05
|  |  |  |  | 8 | 0.74 | 0.12 | 0.72 | 0.13 | 0.73 | 0.15
|  |  |  |  | 10 | 0.46 | 0.14 | 0.40 | 0.13 | 0.39 | 0.13
|  | 16 | 8 | 128 | 4 | 0.99 | 0.03 | 0.99 | 0.02 | 1.00 | 0.01
|  |  |  |  | 6 | 0.94 | 0.05 | 0.91 | 0.07 | 0.77 | 0.11
|  |  |  |  | 8 | 0.66 | 0.11 | 0.69 | 0.11 | 0.69 | 0.12
|  |  |  |  | 10 | 0.44 | 0.12 | 0.37 | 0.12 | 0.35 | 0.11
|  | 16 | 16 | 256 | 4 | 1.00 | 0.003 | 1.00 | 0.006 | 1.00 | 0.003
|  |  |  |  | 6 | 0.98 | 0.03 | 0.95 | 0.05 | 0.96 | 0.04
|  |  |  |  | 8 | 0.77 | 0.12 | 0.65 | 0.12 | 0.62 | 0.17
|  |  |  |  | 10 | 0.40 | 0.15 | 0.23 | 0.08 | 0.28 | 0.10

We compare the performance of bi-ADMM with $L_{2}$-norm (bi-ADMM ($L_{2}$)),
bi-ADMM with $L_{1}$-norm (bi-ADMM ($L_{1}$)), and COBRA on identifying the
bicluster structures in terms of the ARI. The mean value and the standard
deviation of the ARI for each setting and algorithm across 100 replications
are presented in Table LABEL:tab:1. The largest mean ARI value for each
setting is bolded, and all tied values are bolded as well. Overall, we can see
that as $\sigma$ increases, the ARIs for all three methods decrease. These
observations are reasonable because a larger $\sigma$ implies a more
challenging task to identify correct bicluster structures due to a higher
noise level. In the meanwhile, given $\sigma$ and $M$, a larger sample size
implies better performance. As to the differences between three methods, it is
clear that bi-ADMM ($L_{2}$) outperform the others in most cases, while bi-
ADMM ($L_{1}$) leads in a few cases with narrow margins. COBRA works
comparably well in low-noise scenarios, but its performance deteriorates more
quickly than bi-ADMM counterparts as the noise level increases, even worse
than bi-ADMM ($L_{1}$) in some high-noise scenarios. The above results
indicate that the proposed bi-ADMM methods are more robust convex biclustering
algorithms. It is not surprising to see that bi-ADMM ($L_{1}$) is inferior to
the other two methods in most cases, because the spherical clusters are
generated under multivariate Gausiian distributions, which use $L_{2}$-norm
internally.

###  4.2 Compositional Data

In this section, we evaluate the biclustering performance of biC-ADMM with
$L_{2}$-norm and COBRA algorithms using simulated compositional data based on
a murine gut microbiome study (Livanos et al.,, 2016), which compared the gut
microbiome between the group with early-life pulsed therapeutic antibiotic
dosing (PAT) treatment and the group without any antibiotic treatment
(control). In accordance with prior microbiome data analysis (Koh et al.,,
2017; Hu et al.,, 2018), we first generate OTU counts for 100 subjects from
the Dirichlet-multinomial (DM) distribution, where the dispersion parameter
and proportion means were estimated from the control microbiome samples of the
murine study, and the total number of reads per sample is set to 10,000. As a
demonstration, a subset of the OTU count data consisting of 36 mice and 24
common taxa with average proportation $>0.01$ at age of 6 weeks is used in
estimating the parameters of the DM distribution. Here the rule $>0.01$ helps
to remove the extremely rare taxa and to facilitate the following
manipulation. Before we manually manipulate one half of the generated samples
to create biclusters, using the original data we sort the average proportions
of taxa in an ascending order and assign them into 3 groups: Group 1 has the
first six taxa with the smallest proportions, accounting for 3% overall
relative abundance on average; Group 2 has the last five taxa, accounting for
87% overall relative abundance on average; Group 3 has the remaining 13 taxa,
accounting for 10% overall relative abundance on average. We manipulate the
data in the following way. First, for each generated sample using the fitted
DM model, we calculate the ratio of the summed counts for Group 2 over the
summed counted for Group 1 (the averaged ratio is 126), and reduce the ratio
by 1,400 folds. By doing that, we increased the relative abundance for taxa in
Group 1 by roughly 93 folds from the original data. Second, we calculate the
re-scaled counts for Groups 1 and 2 using the adjusted ratio by fixing the
original total counts for Groups 1 and 2. Third, the count for each taxa in
Groups 1 and 2 is re-scaled according to the ratio of the re-scaled group
counts over the original ones. Lastly, the relative abundance for this sample
is obtained by standardizing it using the total counts. By doing these, the
relative abundance of the taxa in Group 3 remains unchanged, denoted as
“unchanged”. Group 1 has its relative abundance enlarged for its taxa, while
Group 2 has reduced relative abundance. Thus, Group 1 is denoted as
“enlarged”, and Group 2 is denoted as “shrunk”. The 50 untouched samples are
denoted as the “control” group, and the 50 manipulated samples are denoted as
the “treatment” group. Then, the resulting data set has two row clusters and
three column clusters.

Next, we evaluate the performance of COBRA and biC-ADMM algorithms on
identifying the bicluster structure in terms of the ARI. The biculster
structure with two row clusters and three column clusters may be artificial,
so the ARI only serves as a reference metric. Under the compositional setting,
a single tuning parameter $\gamma$ for the COBRA is tuned, using the default
setting provided by the R package “cvxbiclustr”, and we adopt a two-
dimensional grid search on $\gamma_{1}$ and $\gamma_{2}$ for the biC-ADMM
method to allow for more flexibility.

Table LABEL:tab:2 shows that the biC-ADMM is superior to the COBRA algorithm
in terms of the ARI. For one simulated data set, the best tuning parameters
are $\gamma_{1}=37.28$ and $\gamma_{2}=268.27$. Figure 2 presents the
corresponding heat map based on the estimated $\widehat{\mathbf{A}}$. We can
see that the biC-ADMM algorithm correctly identifies the “treatment” and
“control” groups for samples. For column clusters, the taxa from the same
groups are almost classified in the same algorithm-generated clusters. The
five taxa in the group “shrunk” are included in two generated clusters with
sizes as 1 and 4. The six taxa in the “enlarged” group are split into two
clusters also, with one mixed into a cluster primarily formed by taxa from the
“unchanged” group. Only one taxon from the “unchanged” group is merged with
taxa from the “shrunk” group.

Table 2: Simulations results for biC-ADMM and COBRA algorithms in terms of the ARI under the compositional setup. biC-ADMM | COBRA
---|---
mean | sd | mean | sd
0.54 | 0.10 | 0.48 | 0.09
Figure 2: Left panel: the heat map of the original simulated compositional
data with 100 samples in 2 row clusters and 24 taxa in 3 column clusters.
Right panel: the heat maps of $\widehat{\mathbf{A}}$ estimated by the biC-ADMM
algorithm. The row(column) side color bars indicate the “true” row(column)
clusters.

###  4.3 Computational Consideration

Although both the proposed Algorithm 1 and COBRA can target the same convex
objective function (10), as we can see from Figure 1, the resulting estimates
vary slightly. The dendrogram or an extra smoothing step can facilitate COBRA
to convey the bicluster pattern, while our method can clearly display the
pattern. One possible reason is that we directly derive the ADMM algorithm,
and COBRA adopts an alternative approximation algorithm.

Step 2 in both Algorithms 1 and 2 involves solving the Sylvester Equation,
which imposes extra computational burden compared with the COBRA based on the
DLPA and the AMA based convex clustering solver, which has less computational
burden compared with the ADMM based solver (Chi and Lange,, 2015; Wang et
al.,, 2018). As we implement the up-to-date algorithm to solve such type of
equation, it is interesting to compare the computational time of the proposed
methods to the COBRA on various data settings. In the following, we evaluate
the computational time in seconds of the COBRA and the bi-ADMM algorithm under
the $L_{2}$-norm. The COBRA is implemented with the R package “cvxbiclustr”
(Chi et al.,, 2017). We have developed an R package “biADMM” to implement our
convex biclustering algorithms, bi-ADMM and biC-ADMM. We also developed a
solver for the Sylvester Equation, which is based on a modified Bartels-
Stewart algorithm proposed by Sorensen et al., (2003). To improve the
efficiency of the bi-ADMM algorithm, we also develop a Python version, which
can be called back in R.

Table 3: The frequency table between true sample groups and clusters obtained
from the biC-ADMM algorithm.

|  $n=50$ |  |  |  |  |
---|---|---|---|---|---|---
$p=40$ |  |  |  |  |  |
$K=4$ |  |  |  |  |  |
|  $n=50$ |  |  |  |  |
$p=40$ |  |  |  |  |  |
$K=8$ |  |  |  |  |  |
|  $n=100$ |  |  |  |  |
$p=80$ |  |  |  |  |  |
$K=4$ |  |  |  |  |  |
|  $n=100$ |  |  |  |  |
$p=80$ |  |  |  |  |  |
$K=8$ |  |  |  |  |  |
|  $n=200$ |  |  |  |  |
$p=160$ |  |  |  |  |  |
$K=12$ |  |  |  |  |  |
|  $n=200$ |  |  |  |  |
$p=160$ |  |  |  |  |  |
$K=16$ |  |  |  |  |  |
| $R=8$ |  |  |  |  |
COBRA | 0.001 | 0.001 | 0.008 | 0.014 | 0.046 | 0.038
bi-ADMM | 0.791 | 0.788 | 7.786 | 7.765 | 104.208 | 104.386
bi-ADMM (Python) | 0.072 | 0.075 | 0.657 | 0.617 | 6.247 | 6.336

The data sets are generated with the simulation settings described in Section
4.1. We set $\sigma=2$ with a given tuning parameter $\gamma=0.1$ defined in
(10). Thirty replicates are simulated, and each algorithm is run for 10
iterations. Table LABEL:tab:3 compares the average running time across 30
replicates of COBRA, bi-ADMM and Python-version bi-ADMM under various
combinations of $n,p,K$ and $R$. The COBRA shows its superiority in program
efficiency, with the shortest running time and the slowest increase as the
scenario gets more complicated. The Python-version bi-ADMM program improves
the running time significantly compared with the R version. It is reasonable
to pay acceptably more running time for better biclustering performance. The
computer is equipped with a CPU i5-7267 (3.10 GHz) and 8G memory.

##  5 Application of biC-ADMM to a Gut Microbiome Data

Here, we apply the proposed biC-ADMM algorithm with two tuning parameters to
the murine microbiome dataset used in simulation studies to assess the
biclustering performance in compositional data. To be specific, we aim to
simultaneously cluster a subset at 13 weeks with 37 common taxa and 68
samples, consisting of 32 microbiome samples from the PAT group and 36 control
samples without retrieving their true treatment group information. The left
panel of Figure LABEL:fig:real_raw_best shows the heat map and dendrogram of
this data set. Since microbes in a community are usually dependent upon one
another and may react similarly with environmental changes, it is of great
interest to investigate whether there are multiple microbial taxa with the
relative abundance being altered concordantly with the antibiotic treatment.
In addition, the relative abundances of microbiome data are highly skewed,
with a few common taxa accounting for the majority of bacteria in the sample.
Microbiome studies usually need to group taxa at a higher rank to pull the
signal strength. Some other research found that microbial community can be
divided into sub-communities based on their functional reactions to the
perturbation (Holmes et al.,, 2012; Sankaran and Holmes,, 2019). All these
observations motivate the application of clustering algorithm to identify
data-driven taxa cluster patterns.

Following the commonly used pre-processing step in microbiome data analysis,
the OTU count data of each sample is normalized into relative abundances prior
to conducting biclustering, i.e., with relative abundances of 37 taxa summing
up to 1. In our biC-ADMM algorithm, we set $m_{1}=m_{2}=5$ and $\phi=1$ for
weight $w_{l}$ and $u_{k}$, and set $\nu_{1}=\nu_{2}=\nu_{3}=1$.

Then we apply the biC-ADMM algorithm to conduct convex biclustering for this
data, which is an unsupervised learning method, free of a priori grouping
knowledge. Although the true sample cluster labels are known, we mask them in
biclustering analysis and aim to identify taxa clusters associated with
potential sample subgroups. The tuning parameters are selected to maximize the
stability for taxa clusters based on stability selection in Section 3.2. The
resulting optimal tuning parameters are $\gamma_{1}=153.99$ and
$\gamma_{2}=82.54$. The bic-ADMM algorithm results into eight taxa clusters
and seven sample clusters. It is noted that more sample clusters are obtained
using biC-ADMM with an ARI of 0.40, but from Table LABEL:tab:real_cross we can
see the purity of biC-ADMM generated clusters is high with only five samples
are misclassified. The resulting taxa clusters indicate those taxa subgroups
which have the potential to discriminate sample subgroups.

Next, we illustrate how the solution $\widehat{\mathbf{A}}$ evolves as
$\gamma_{1}$ and $\gamma_{2}$ vary. Figure LABEL:fig:real_varying_gamma shows
snap shots of the biC-ADMM solution path of this data set, as the parameter
$\gamma_{1}$ increases by row and $\gamma_{2}$ increases by column. The row
path captures the whole range of behavior between under-smoothed estimates of
the mean structure for taxa clusters (small $\gamma_{2}$) to over-smoothed
estimates (large $\gamma_{2}$) with fixed sample clusters, while the row
column captures the whole range of behavior between under-smoothed estimates
for sample clusters (small $\gamma_{1}$) to over-smoothed estimates (large
$\gamma_{1}$) with fixed taxa clusters. In between the extremes in top-left
and bottom-right panels, rows and columns are fusing together as $\gamma_{1}$
and $\gamma_{2}$ increases. The entire solution path enable us to identify
those specific taxa that could discriminate sample groups. “Akkermansia”,
“Lactobacillus”, “Rikenellaceae—Other”, “Blautia” and “S24-7—Other” seem to
play an important role, and “S24-7—Other” accounts for the largest relative
abundance.

| biC-ADMM clusters
---|---
| 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8
Control | 1 | 2 | 0 | 0 | 1 | 1 | 25 | 6
PAT | 17 | 4 | 3 | 5 | 0 | 0 | 0 | 0

##  6 Summary

Biclustering is a powerful data mining technique which can provide results in
a checkerboard-like pattern to explore the feature patterns for subgroups,
while it is under-utilized in biological and biomedical data. Among
biclustering algorithms developed in the past two decades, the convex
biclustering algorithm is appealing due to its globally optimal solution.
However, existing convex biclustering solved the convex optimization problem
via an approximation to objective function of the ADMM algorithm, which
required extra smoothing steps to obtain informative biclusters. On the other
hand, existing biclustering methods cannot incorporate the compositional
constraints which are often needed for the increasingly popular microbiome
data.

In this manuscript, a new algorithm for the standard convex biclustering and
its extension under the compositional constraints are proposed to
simultaneously cluster observations and features, denoted as the bi-ADMM and
the biC-ADMM, respectively. Both algorithms directly apply the ADMM method
with solving a Sylvester Equation as a key step, which has a great potential
to solve other problems with both row and columns constraints for a data
matrix. The numerical results show that bi-ADMM and biC-ADMM are able to
provide clear checkerboard-like pattern directly without any further smoothing
step. Furthermore, the biC-ADMM can solve biclustering problems with
compositional constraints, and it can be extended to any linear constraint
easily, which has a great potential to be used in microbiome studies.
Furthermore, we have developed an R package “biADMM” to implement our convex
biclustering algorithms, bi-ADMM and biC-ADMM, which includes Python-version
functions to speed up the calculation.

Moreover, this work can motivate future research. Wang et al., (2018)
presented the Sparse Convex Clustering problem, seeking to cluster
observations and select features simultaneously. Feature selection can be
incorporated into the convex biclustering as well. In addition, Zhou et al.,
(2020) developed a very efficient smoothing proximal gradient algorithm
(Sproga) for convex clustering. Then, we can extend convex biclustering to
sparse convex biclustering and develop a highly efficient algorithms to solve
it, in order to conduct biclustering and feature selection simultaneously.
Furthermore, weights $w_{l}$ and $u_{k}$ serves important roles in convex
biclustering, but they suffer from the variable selection problem (Chakraborty
and Xu,, 2020). The proposed algorithms can be improved in this direction as
well.

## Acknowledgement

Hu and Li’s research is partially supported by a National Institute of Health
(NIH) grant (5R01DK110014).

## References

  * Bartels and Stewart, (1972) Bartels, R. H. and Stewart, G. W. (1972).  Solution of the matrix equation ax+ xb= c [f4].  Communications of the ACM, 15(9):820–826.
  * Busygin et al., (2008) Busygin, S., Prokopyev, O., and Pardalos, P. M. (2008).  Biclustering in data mining.  Computers & Operations Research, 35(9):2964–2987.
  * Chakraborty and Xu, (2020) Chakraborty, S. and Xu, J. (2020).  Biconvex clustering.
  * Cheng and Church, (2000) Cheng, Y. and Church, G. M. (2000).  Biclustering of expression data.  In Ismb, volume 8, pages 93–103.
  * Chi et al., (2017) Chi, E. C., Allen, G. I., and Baraniuk, R. G. (2017).  Convex biclustering.  Biometrics, 73(1):10–19.
  * Chi and Lange, (2015) Chi, E. C. and Lange, K. (2015).  Splitting Methods for Convex Clustering.  Journal of Computational and Graphical Statistics, 24(4):994–1013.
  * Falony et al., (2016) Falony, G., Joossens, M., Vieira-Silva, S., Wang, J., Darzi, Y., Faust, K., Kurilshikov, A., Bonder, M. J., Valles-Colomer, M., Vandeputte, D., et al. (2016).  Population-level analysis of gut microbiome variation.  Science, 352(6285):560–564.
  * Fang and Wang, (2012) Fang, Y. and Wang, J. (2012).  Selection of the number of clusters via the bootstrap method.  Computational Statistics & Data Analysis, 56:468–477.
  * Flynn and Perry, (2012) Flynn, C. J. and Perry, P. O. (2012).  Consistent biclustering.  arXiv preprint arXiv:1206.6927, 3.
  * Friedman et al., (2001) Friedman, J., Hastie, T., and Tibshirani, R. (2001).  The elements of statistical learning, volume 1.  Springer series in statistics New York.
  * Gloor et al., (2017) Gloor, G. B., Macklaim, J. M., Pawlowsky-Glahn, V., and Egozcue, J. J. (2017).  Microbiome datasets are compositional: and this is not optional.  Frontiers in microbiology, 8:2224.
  * Hochreiter et al., (2010) Hochreiter, S., Bodenhofer, U., Heusel, M., Mayr, A., Mitterecker, A., Kasim, A., Khamiakova, T., Van Sanden, S., Lin, D., Talloen, W., et al. (2010).  Fabia: factor analysis for bicluster acquisition.  Bioinformatics, 26(12):1520–1527.
  * Holmes et al., (2012) Holmes, I., Harris, K., and Quince, C. (2012).  Dirichlet multinomial mixtures: generative models for microbial metagenomics.  PloS one, 7(2):e30126.
  * Hripcsak and Albers, (2012) Hripcsak, G. and Albers, D. J. (2012).  Next-generation phenotyping of electronic health records.  Journal of the American Medical Informatics Association, 20(1):117–121.
  * Hu et al., (2018) Hu, J., Koh, H., He, L., Liu, M., Blaser, M. J., and Li, H. (2018).  A two-stage microbial association mapping framework with advanced fdr control.  Microbiome, 6(1):131.
  * Hubert and Arabie, (1985) Hubert, L. and Arabie, P. (1985).  Comparing partitions.  Journal of classification, 2(1):193–218.
  * Jameson, (1968) Jameson, A. (1968).  Solution of the equation ax+xb=c by inversion of an m*m or n*n matrix.  SIAM Journal on Applied Mathematics, 16(5):1020–1023.
  * Koh et al., (2017) Koh, H., Blaser, M. J., and Li, H. (2017).  A powerful microbiome-based association test and a microbial taxa discovery framework for comprehensive association mapping.  Microbiome, 5(1):45.
  * Lee et al., (2010) Lee, M., Shen, H., Huang, J. Z., and Marron, J. S. (2010).  Biclustering via Sparse Singular Value Decomposition.  Biometrics, 66(4):1087–1095.
  * Li, (2020) Li, G. (2020).  Generalized co-clustering analysis via regularized alternating least squares.  Computational statistics & data analysis, 150:106989.
  * Livanos et al., (2016) Livanos, A. E., Greiner, T. U., Vangay, P., Pathmasiri, W., Stewart, D., McRitchie, S., Li, H., Chung, J., Sohn, J., Kim, S., et al. (2016).  Antibiotic-mediated gut microbiome perturbation accelerates development of type 1 diabetes in mice.  Nature microbiology, 1(11):1–13.
  * Madeira and Oliveira, (2004) Madeira, S. C. and Oliveira, A. L. (2004).  Biclustering algorithms for biological data analysis: a survey.  IEEE/ACM Transactions on Computational Biology and Bioinformatics (TCBB), 1(1):24–45.
  * Sankaran and Holmes, (2019) Sankaran, K. and Holmes, S. P. (2019).  Latent variable modeling for the microbiome.  Biostatistics, 20(4):599–614.
  * Shabalin et al., (2009) Shabalin, A. A., Weigman, V. J., Perou, C. M., and Nobel, A. B. (2009).  Finding large average submatrices in high dimensional data.  The Annals of Applied Statistics, pages 985–1012.
  * Sorensen et al., (2003) Sorensen, D. C., Zhou, Y., et al. (2003).  Direct methods for matrix sylvester and lyapunov equations.  Journal of Applied Mathematics, 6(2003):277–303.
  * Tan and Witten, (2014) Tan, K. M. and Witten, D. M. (2014).  Sparse biclustering of transposable data.  Journal of Computational and Graphical Statistics, 23(4):985–1008.  PMID: 25364221.
  * Tibshirani et al., (2005) Tibshirani, R., Saunders, M., Rosset, S., Zhu, J., and Knight, K. (2005).  Sparsity and smoothness via the fused lasso.  Journal of the Royal Statistical Society: Series B, 67:1198–1232.
  * Wang et al., (2018) Wang, B., Zhang, Y., Sun, W. W., and Fang, Y. (2018).  Sparse convex clustering.  Journal of Computational and Graphical Statistics, 27:393–403.
  * Xie et al., (2018) Xie, J., Ma, A., Fennell, A., Ma, Q., and Zhao, J. (2018).  It is time to apply biclustering: a comprehensive review of biclustering applications in biological and biomedical data.  Briefings in bioinformatics, 1:16.
  * Zhou et al., (2020) Zhou, X., Du, C., and Cai, X. (2020).  An efficient smoothing proximal gradient algorithm for convex clustering.  arXiv preprint arXiv:2006.12592.

## Appendix

In Appendix, we provide proofs of Lemma 1 and Lemma 2.

Proof of Lemma 1

Denote ${\bf a}=\textrm{vec}({\bf A})$, a vectorization of the matrix ${\bf
A}$. According to the fact that $A_{i_{1}\cdot}-A_{i_{2}\cdot}={\bf A}^{\rm
T}({\bf e}_{i_{1}}-{\bf e}_{i_{2}})$,
$\mathbf{a}_{k_{1}}-\mathbf{a}_{k_{2}}=\mathbf{A}(\mathbf{e}_{k_{1}}^{*}-\mathbf{e}_{k2}^{*})$,
and the property of the tensor product $\textrm{vec}({\bf R}{\bf S}{\bf
T})=[{\bf T}^{\rm T}\otimes{\bf R}]\textrm{vec}({\bf S})$, solving the
minimization of $f({\bf A})$ is equivalent to minimize

$\displaystyle f(\mathbf{a})={1\over
2}\|\mathbf{x}-\mathbf{a}\|_{2}^{2}+{\nu_{1}\over 2}\sum_{l\in{\cal
E}_{l}}\|\mathbf{B}_{l}\mathbf{P}\mathbf{a}-\widetilde{\mathbf{v}}_{l}\|_{2}^{2}+{\nu_{2}\over
2}\sum_{k\in{\cal
E}_{2}}\|\mathbf{C}_{k}\mathbf{a}-\widetilde{\mathbf{z}}_{k}\|_{2}^{2},$

where $\mathbf{B}_{l},\mathbf{C}_{k}$ and $\mathbf{P}$ are defined in Section
2.3. Then it follows that

$\displaystyle f(\mathbf{a})={1\over
2}\|\mathbf{x}-\mathbf{a}\|_{2}^{2}+{\nu_{1}\over
2}\|\mathbf{B}\mathbf{P}\mathbf{a}-\widetilde{\mathbf{v}}\|_{2}^{2}+{\nu_{2}\over
2}\|\mathbf{C}\mathbf{a}-\widetilde{\mathbf{z}}\|_{2}^{2},$

where $\mathbf{B},\mathbf{C},\widetilde{\mathbf{v}}$ and
$\widetilde{\mathbf{z}}$ are defined in Section 2.3 as well.

Finally, the stationary equation can be obtained by

$\displaystyle(\mathbf{I}_{np}+\nu_{1}\mathbf{P}^{\rm T}\mathbf{B}^{\rm
T}\mathbf{B}\mathbf{P}+\nu_{2}\mathbf{C}^{\rm
T}\mathbf{C})\mathbf{a}=\mathbf{x}+\nu_{1}\mathbf{P}^{\rm T}\mathbf{B}^{\rm
T}\widetilde{\mathbf{v}}+\nu_{2}\mathbf{C}^{\rm T}\widetilde{\mathbf{z}}.$

This is a system of $np$ linear equations. By applying the properties for
Kronecker product such as $(\mathbf{S}\otimes\mathbf{T})^{\rm
T}=\mathbf{S}^{\rm T}\otimes\mathbf{T}^{\rm T}$ and
$(\mathbf{Q}\otimes\mathbf{R})(\mathbf{S}\otimes\mathbf{T})=(\mathbf{Q}\mathbf{S})\otimes(\mathbf{R}\mathbf{T})$,
it follows

$\displaystyle\nu_{2}\mathbf{C}^{\rm T}\mathbf{C}$ $\displaystyle=$
$\displaystyle\nu_{2}\sum_{k\in{\cal
E}_{2}}\left[\left((\mathbf{e}_{k_{1}}^{*}-\mathbf{e}_{k_{2}}^{*})(\mathbf{e}_{k_{1}}^{*}-\mathbf{e}_{k_{2}}^{*})^{\rm
T}\right)\right]\otimes\mathbf{I}_{n}$ $\displaystyle=$
$\displaystyle\left[\sum_{k\in{\cal
E}_{2}}\nu_{2}\left((\mathbf{e}_{k_{1}}^{*}-\mathbf{e}_{k_{2}}^{*})(\mathbf{e}_{k_{1}}^{*}-\mathbf{e}_{k_{2}}^{*})^{\rm
T}\right)\right]\otimes\mathbf{I}_{n}$ $\displaystyle\nu_{2}\mathbf{C}^{\rm
T}\widetilde{\mathbf{z}}$ $\displaystyle=$
$\displaystyle\nu_{2}\sum_{k\in{\cal
E}_{2}}[(\mathbf{e}_{k_{1}}^{*}-\mathbf{e}_{k_{2}}^{*})\otimes\mathbf{I}_{n}]\widetilde{\mathbf{z}}_{k}$
$\displaystyle=$ $\displaystyle\sum_{k\in{\cal
E}_{2}}\left(\nu_{2}(\mathbf{e}_{k_{1}}^{*}-\mathbf{e}_{k_{2}}^{*})\otimes\mathbf{I}_{n}\right)\widetilde{\mathbf{z}}_{k}$
$\displaystyle\nu_{3}\mathbf{D}^{\rm T}\mathbf{D}$ $\displaystyle=$
$\displaystyle\nu_{3}\sum_{i=1}^{n}\mathbf{D}_{i}^{\rm
T}\mathbf{D}_{i}=\nu_{3}({\bf 1}_{p}{\bf 1}_{p}^{\rm T})\otimes\mathbf{I}_{n}$
$\displaystyle\nu_{3}\mathbf{D}^{\rm T}{\bf s}$ $\displaystyle=$
$\displaystyle\nu_{3}\sum_{i=1}^{n}({\bf
1}_{p}\otimes\mathbf{e}_{i})s_{i}=\nu_{3}\sum_{i=1}^{n}s_{i}\textrm{vec}(\mathbf{e}_{i}{\bf
1}_{p}^{\rm T}).$

In addition, by applying Proposition 1 in Wang et al., (2018) for the
permutation matrix $\mathbf{P}$, we can prove

$\displaystyle\mathbf{I}_{np}+\nu_{1}\mathbf{P}^{\rm T}\mathbf{B}^{\rm
T}\mathbf{B}\mathbf{P}$ $\displaystyle=$
$\displaystyle\mathbf{I}_{p}\otimes\left[\mathbf{I}_{n}+\nu_{1}\sum_{l\in{\cal
E}_{1}}(\mathbf{e}_{l_{1}}-\mathbf{e}_{l_{2}})(\mathbf{e}_{l_{1}}-\mathbf{e}_{l_{2}})^{\rm
T}\right]$ $\displaystyle\nu_{1}\mathbf{P}^{\rm T}\mathbf{B}^{\rm
T}\widetilde{\mathbf{v}}$ $\displaystyle=$ $\displaystyle\sum_{l\in{\cal
E}_{1}}\left(\mathbf{I}_{p}\otimes\nu_{1}(\mathbf{e}_{l_{1}}-\mathbf{e}_{l_{2}})\right)\widetilde{\mathbf{v}}_{l}.$

Therefore, the system of equations is equivalent to

$\displaystyle(\mathbf{I}_{p}\otimes\mathbf{M})\textrm{vec}(\mathbf{A})+(\mathbf{N}\otimes\mathbf{I}_{n})\textrm{vec}(\mathbf{A})=\textrm{vec}(\mathbf{G}),$

which is equivalent to $\mathbf{M}\mathbf{A}+\mathbf{A}\mathbf{N}=\mathbf{G}$.
$\Box$

Proof of Lemma 2

Proof of Lemma 1 needs to be modified to incorporate the extra constraints.
The stationary equations for minimizing (8) can be obtained by

$\displaystyle(\mathbf{I}_{np}+\nu_{1}\mathbf{P}^{\rm T}\mathbf{B}^{\rm
T}\mathbf{B}\mathbf{P}+\nu_{2}\mathbf{C}^{\rm
T}\mathbf{C}+\nu_{3}\mathbf{D}^{\rm
T}\mathbf{D})\mathbf{a}=\mathbf{x}+\nu_{1}\mathbf{P}^{\rm T}\mathbf{B}^{\rm
T}\widetilde{\mathbf{v}}+\nu_{2}\mathbf{C}^{\rm
T}\widetilde{\mathbf{z}}+\nu_{3}\mathbf{D}^{\rm T}{\bf s}.$

With matrix techniques, it follows

$\displaystyle\nu_{3}\mathbf{D}^{\rm T}\mathbf{D}$ $\displaystyle=$
$\displaystyle\nu_{3}\sum_{i=1}^{n}\mathbf{D}_{i}^{\rm
T}\mathbf{D}_{i}=\nu_{3}({\bf 1}_{p}{\bf 1}_{p}^{\rm T})\otimes\mathbf{I}_{n}$
$\displaystyle\nu_{3}\mathbf{D}^{\rm T}{\bf s}$ $\displaystyle=$
$\displaystyle\nu_{3}\sum_{i=1}^{n}({\bf
1}_{p}\otimes\mathbf{e}_{i})s_{i}=\nu_{3}\sum_{i=1}^{n}s_{i}\textrm{vec}(\mathbf{e}_{i}{\bf
1}_{p}^{\rm T}).$

Therefore, the system of equations is equivalent to (I_p ⊗M)vec(A) + (N⊗I_n
)vec(A) = vec(G). $\Box$