# Estimation of Gaussian Bi-Clusters with General Block-Diagonal Covariance
Matrix and Applications

Anastasiia Livochka 111Department of Statistics and Actuarial Science, 200
University Ave W, Waterloo, ON, N2L 3G1, Canada Ryan Browne 222Department of
Statistics and Actuarial Science, 200 University Ave W, Waterloo, ON, N2L 3G1,
Canada. e<EMAIL_ADDRESS>Sanjeena Subedi333School of Mathematics
and Statistics, 4302 Herzberg Laboratories, Carleton University, 1125 Colonel
By Drive, Ottawa, ON, K1S 5B6, Canada e<EMAIL_ADDRESS>
###### Abstract

Bi-clustering is a technique that allows for the simultaneous clustering of
observations and features in a dataset. This technique is often used in
bioinformatics, text mining, and time series analysis. An important advantage
of biclustering algorithm is the ability to uncover multiple “views” (i.e.,
through rows and column groupings) in the data. Several Gaussian mixture model
based biclustering approach currently exist in the literature. However, they
impose severe restrictions on the structure of the covariance matrix. Here, we
propose a Gaussian mixture model-based bi-clustering approach that provides a
more flexible block-diagonal covariance structure. We show that the clustering
accuracy of the proposed model is comparable to other known techniques but our
approach provides a more flexible covariance structure and has substantially
lower computational time. We demonstrate the application of the proposed model
in bioinformatics and topic modelling.

Keywords: Bi-clustering, Gaussian mixture models

##  1 Introduction

Humans are particularly good at classifying objects into natural groups based
on some notion of similarity and the goal they want to achieve. However, as
the set of objects to organize and the number of features grows, we require
faster and more memory-efficient tools to complete the task. Cluster analysis
is an important area of unsupervised machine learning that deals with the
problem of grouping objects so that objects in one cluster are more similar to
each other than to objects in other clusters (Diday and Simon, 1976). The main
applications of cluster analysis include data exploration and pattern
recognition (Berkhin, 2006).

Cluster analysis has been widely studied for decades, and many clustering
algorithms have been developed. Interestingly, the great variety of the
available models can be attributed to the fact that precisely defining a
cluster presents a problem (Milligan and Cooper, 1987). Some algorithms (e.g.
DBSCAN, HDBSCAN) leverage density assumptions by defining clusters as
connected dense regions in the space (Khan et al., 2014; McInnes et al.,
2017). Other approaches view clusters as mixtures of certain distributions,
with Gaussian mixture model (GMM) being the most popular one (Reynolds, 2009).

Bi-clustering is a technique that allows for the simultaneous clustering of
observations and features in a dataset (Hartigan, 1972). An important
advantage of bi-clustering is the ability to uncover multiple “views” (i.e.,
through row and column groupings) in the data. This technique is often used in
bioinformatics, text mining, and time series analysis (Madeira and Oliveira,
2004; Madeira et al., 2008; Busygin et al., 2008). For example, Wang et al.
(2013) utilized biclustering algorithm on gene expression data and identified
subgroups of breast cancer tumour with similar clinical characteristics. In
topics modelling context, it’s been extensively used to interpret the
resulting clusters (Rugeles et al., 2017).

There is no single definition of a cluster and as such there is no single
definition of a bi-cluster. Therefore, a variety of bi-clustering models
exist. One of the first applications of bi-clustering to text documents is
based on the spectral graph (Dhillon, 2001) partitioning heuristic. Similarly,
Kluger et al. (2003) shows that the checkerboard structure of the data matrix
can be inferred from its eigenvectors. There are many greedy approaches
proposed for the co-clustering of gene expressions (Padilha and Campello,
2017; Cheng and Church, 2000; Bergmann et al., 2003). For instance, Ben-Dor et
al. (2002) introduces an approach that finds the hidden order-preserving sub-
matrices in the data matrix.

Alternate model-based biclustering approaches assume that the observations in
the same row cluster were generated from the same distribution. e.g.
$\mathcal{N}_{p}(\bm{\mu},\bm{\Sigma})$. In turn, column groups can be
inferred from the structure of the covariance matrix. Suppose we have a
$p$-dimensional mean vector $\bm{\mu}$ such that

$\bm{\mu}=\left[\bm{\mu}_{1},\bm{\mu}_{2},\bm{\mu}_{3}\right],$

where the dimesions of $\bm{\mu}_{i}$ is $c_{i}$, $i=1,\ldots,3$ and
$p=c_{1}+c_{2}+c_{3}$, and covariance covariance matrix, $\bm{\Sigma}$, can be
written as:

$\bm{\Sigma}=\,\mbox{diag}\left(\bm{\Sigma}_{1},\bm{\Sigma}_{2},\bm{\Sigma}_{3}\right)=\left[\begin{array}[]{ccc}\bm{\Sigma}_{1(c_{1}\times
c_{1})}&\mathbf{0}&\mathbf{0}\\\ \mathbf{0}&\bm{\Sigma}_{2(c_{2}\times
c_{2})}&\mathbf{0}\\\ \mathbf{0}&\mathbf{0}&\bm{\Sigma}_{3(c_{3}\times
c_{3})}\end{array}\right].$ (1)

Then ($\bm{\mu}_{1}$,
$\bm{\Sigma}_{1}),(\bm{\mu}_{2},\bm{\Sigma}_{2}),(\bm{\mu}_{3},\bm{\Sigma}_{3}$)
define three distinct column groupings. Note that column permutations may be
required to achieve a checkerboard-like structure as illustrated in 1. This
adds to the non-triviality of the estimation of block-diagonal covariance
matrix.

Figure 1: Illustration of the block diagonal matrices with off-diagonal
elements limited to positive (left) and general (right) structure.

Initial model-based bi-clustering approaches introduced additional latent
variables to indicate column clusters (Bhatia et al., 2014; Gu and Liu, 2008).
Gallaugher et al. (2022) proposed a Gaussian model-based clustering for high-
dimensional data. They partition columns of the dataset twice - once by means
and another by variances. Others (Martella et al., 2008; Wong et al., 2017; Tu
and Subedi, 2022) leveraged the bi-clustering framework rooted in the mixtures
of factor analyzers (McNicholas and Murphy, 2008) approach . The advantage of
the factor analyzer structure is that it decreases the number of estimated
parameters by assuming a latent low-dimensional representation of high-
dimensional variables. However, this representation restricts the covariance
matrix structure that can be recovered. For example, Martella et al. (2008)
and Wong et al. (2017) require the within block off-diagonal elements of the
covariance matrix to $1$. Tu and Subedi (2022) relaxed some constraints by
characterizing covariance to allow the all of the off-diagonal entries of
$i$-th covariance block to be equal to a positive constant. This is
illustrated in the left panel of Figure 1. In addition, the parameter
estimation utilized in Martella et al. (2008); Wong et al. (2017); Tu and
Subedi (2022) has limitations due to the complexity of finding column
groupings.

In this work, we aim to build on the model-based bi-clustering procedure by
allowing the covariance matrix of the underlying mixture components to have a
general block-diagonal structure as shown in the right panel of Figure 1. We
derive and compare three approaches to estimate the column groupings; greedy,
convex relaxation and hierarchical. Empirically the hierarchical approach is
shown to be the best approach in terms of accuracy and computational
efficiency. We also compare the hierarchical approach to readily available
covariance matrix estimators (Ledoit and Wolf, 2004; Friedman et al., 2008;
Broto et al., 2022). The readily available approaches do not satisfy our
desirable conditions for the covariance estimator, namely ($i$) enforcing the
block-diagonal structure on the covariance matrix; and ($ii$) allowing for
general covariance matrix structure inside the blocks. In the second half the
manuscript, we introduce the model-based bi-clustering procedure with the
general block-diagonal structure. We use simulated datasets to determine the
efficacy of the model. We show utility of our approach in bioinformatics and
topics modelling. First, we applied the proposed model and other existing bi-
clustering models on two high-dimensional genomics classification problems:
Alon (Alon et al., 1999) and Golub (Golub et al., 1999) datasets. Then, we
illustrate the proposed methodology on two datasets from topic modelling:
hotel ratings & reviews (Alam et al., 2016) and Disneyland ratings & reviews
(Chillar, 2022).

##  2 Estimating a Block-Diagonal Covariance Matrix

This section discusses parameterization of a block-diagonal matrix and its
estimation. We propose three methods for estimating the block-diagonal
covariance matrix (as illustrated on the 1 \- which is crucial for the GMM-
based bi-clustering procedure given section 3. For simplicity, we derive an
estimation for one component of GMM here. In Section 3, we extend the
techniques for multi-component mixture model. The main aim of this work is to
propose an approach such that: $(i)$ it imposes a block-diagonal structure
onto the covariance matrix; $(ii)$ it does not restrict the covariance matrix
structure inside the blocks.

###  2.1 Imposing Block-Diagonal Structure Onto a Covariance Matrix

Suppose we have a normally distributed data
$\mathbf{X}\sim\mathcal{N}_{p}(\bm{\mu},\bm{\Sigma})$. One way to compute
feature clusters, common for biclustering, is by estimating a block-diagonal
covariance matrix that is close to the original $\bm{\Sigma}$. In this work,
we propose to parameterize this block-diagonal covariance matrix as:

$\mathbf{D}_{1}^{\intercal}\bm{\Sigma}\mathbf{D}_{1}+\dots+\mathbf{D}_{K}^{\intercal}\bm{\Sigma}\mathbf{D}_{K}\qquad\mbox{where
}\qquad\mathbf{D}_{1}+\dots+\mathbf{D}_{K}=\mathbf{I}_{p}$ (2)

where, $K$ is the number of variable or column groups,
$\mathbf{D}_{k}=\,\mbox{diag}(\mathbf{d}_{k})$ and each element of
$\mathbf{d}_{k}$ is zero or one. In other words, the $\mathbf{D}_{k}$ are
binary diagonal matrices that sum up to the identity matrix. This allows for
each $\mathbf{D}_{k}$ to define a column or feature group $k$ using
$\mathbf{D}=\left[\mathbf{d}_{1},\ldots,\mathbf{d}_{K}\right]$. For example,
to get the general covariance structure in Figure 1, we set

$\mathbf{D}=\left[\begin{array}[]{ccc}\mathbf{d}_{1}&\mathbf{d}_{2}&\mathbf{d}_{3}\end{array}\right]=\left[\begin{array}[]{cccccccc}1&1&1&0&0&0&0&0\\\
0&0&0&1&1&1&0&0\\\ 0&0&0&0&0&0&1&1\end{array}\right]^{\intercal}.$ (3)

This proposed parametrization imposes a general block-diagonal structure on
the covariance matrix. i.e. we can obtain the general covariance structure
shown in the right panel of Figure 1.

If we have a sample $\\{\mathbf{x}_{1},\mathbf{x}_{2},\dots\mathbf{x}_{N}\\}$
generated from a $p$-dimensional multivariate Gaussian with mean, $\bm{\mu}$,
and a postive definite covariance matrix, $\bm{\Sigma}$. The scaled or
normalized log-likelihood is

$-\frac{1}{2}\log|\bm{\Sigma}|-\frac{1}{2N}\sum_{i=1}^{N}\,\mbox{tr}\left[\bm{\Sigma}^{-1}(\mathbf{x}_{i}-\hat{\bm{\mu}})(\mathbf{x}_{i}-\hat{\bm{\mu}})^{\intercal}\right].$

The maximum likelihood estimates (MLEs) of the parameters $\bm{\mu}$ and
$\bm{\Sigma}$ are the sample mean and sample covariance:

$\hat{\bm{\mu}}=\overline{\mathbf{x}}=\frac{\sum_{i=1}^{N}\mathbf{x}_{i}}{N},\qquad\mbox{and}\qquad\mathbf{S}=\frac{\sum_{i=1}^{N}(\mathbf{x}_{i}-\hat{\bm{\mu}})(\mathbf{x}_{i}-\hat{\bm{\mu}})^{\intercal}}{N}.$

Furthermore, we denote sample correlation matrix by $\mathbf{S}^{\star}$.

Since, the MLE for unconstrained covariance matrix is $\mathbf{S}$, the log-
likelihood function for the block diagonal covariance matrix as a function of
the feature groups, $\mathbf{D}=[\mathbf{d}_{1},\dots\mathbf{d}_{K}]$, is

$\mathcal{L}(\mathbf{d}_{1},\ldots,\mathbf{d}_{K})=-\frac{1}{2}\log\left|\sum_{k=1}^{K}\mathbf{D}_{k}\mathbf{S}\mathbf{D}_{k}\right|-\frac{1}{2}\,\mbox{tr}\left[\left(\sum_{k=1}^{K}\mathbf{D}_{k}\mathbf{S}\mathbf{D}_{k}\right)^{-1}\mathbf{S}\right]$
(4)

where $\mathbf{D}_{k}=\,\mbox{diag}(\mathbf{d}_{k})$. In the next three
sections, we discuss three approaches to estimate $\mathbf{D}$: greedy, convex
relaxation (numerical) and hierarchical.

###  2.2 Greedy Approach

Motivated by the approach in Tu and Subedi (2022), we decided to adopt a
greedy algorithm for estimating
$\mathbf{D}=[\mathbf{d}_{1},\dots\mathbf{d}_{K}]$ in (4). To initialize we
follow Tu and Subedi (2022) and set $\mathbf{D}$ to the first $K$ principal
components of the correlation matrix $\mathbf{S}^{\star}$. Then, in the
$j^{th}$ row of $\mathbf{D}$, we find the largest element, e.g. the $k^{th}$
element, set that element to one and then the remaining are set zero. This can
be interpreted as assigning the $j$-th row to the $k^{th}$-the column group.

After initialization, we iteratively reassign column clusters similar to the
initialization but instead we use the largest log-likelihood value (4).
Starting with variable $j$, we evaluate the equation (4) by considering this
feature or variable in each of the $K$ feature group and assign the variable
to the $k^{th}$ column group that maximizes the log-likelihood. Alternatively,
this is equivalent to considering each of the unit vectors, of length $K$ in
the row $j^{th}$ of $\mathbf{D}$. Mathematically, for the $j^{th}$ row of
$\mathbf{D}$, denoted by $\mathbf{D}_{j}$, we have

$\mathbf{D}_{j}=\operatorname*{arg\,max}_{\mathbf{a}\in\\{\mathbf{e}_{1},\ldots,\mathbf{e}_{K}\\}}\mathcal{L}(\mathbf{d}_{1},\ldots,\mathbf{d}_{j-1},\mathbf{a},\mathbf{d}_{j+1}\ldots,\mathbf{d}_{K})$

for $j=1,\ldots,p$ and $\mathbf{e}_{k}$ is a unit vector of length $K$ which
has a unit in the $k^{th}$ position and zero elsewhere. This means each
updated row is conditional on the rest of them.

###  2.3 Convex Relaxation Approach

We aim to find
$\mathbf{D}=\left[\mathbf{d}_{1},\mathbf{d}_{2},\dots,\mathbf{d}_{K}\right]$
that maximizes log-likelihood in (4). Keeping in mind that $\mathbf{d}_{k}$
should be binary by definition of (2), this is a non-linear integer problem,
which is known for its complexity (Hemmecke et al., 2009). To simplify, we
relaxed the $\mathbf{d}$ binary constraints and use quadratic programming to
solve the following:

$\max\mathcal{L}\quad\text{st.}\quad\forall k:\
\mathbf{d}_{k}>\mathbf{0}\quad\text{and}\quad\sum_{k=1}^{K}\mathbf{d}_{k}=\mathbf{1}$
(5)

where $\mathbf{d}_{k}>\mathbf{0}$ means that each element of $\mathbf{d}_{k}$
should be greater than zero. This relaxes the condition that the
$\mathbf{D}_{k}$ are binary diagonal matrices by allow each element of
$\mathbf{D}_{k}$ to be positive. A possible solution is that $\forall k:\
\mathbf{d}_{k}=1/K\times\mathbf{1}$ as it results in $\bm{\Sigma}=\mathbf{S}$.
To avoid the trivial solution, we added a penalty on the logarithm of the
determinant of $\mathbf{D}^{\intercal}\mathbf{D}$. Then in addition we will
use a penalty to push $\mathbf{d}_{k}$ towards binary values. In particular
the penalty is $\sum_{k=1}^{K}\mathbf{d}_{k}^{\intercal}\mathbf{d}_{k}-p$.
Finally, instead of applying a block diagonal structure to the covariance, we
apply it to the inverse of the covariance matrix also known as the precision
matrix. Note that the two parameterizations are equivalent because the inverse
of a block diagonal matrix is block diagonal matrix. However, using the
precision matrix will simply the derivatives.

The parameterized block diagonal precision matrix is

$\mathbf{D}_{1}\bm{\Sigma}^{-1}\mathbf{D}_{1}+\dots+\mathbf{D}_{K}\bm{\Sigma}^{-1}\mathbf{D}_{K}.$
(6)

We want to maximize the log-likelihood function with respect to the block-
diagonal structure of the covariance matrix: The MLE for $\bm{\Sigma}^{-1}$ is
$\mathbf{S}^{-1}$ so putting (4) in terms of the precision matrix we have

$\mathcal{L}^{\star}(\mathbf{d}_{1},\ldots,\mathbf{d}_{K})\frac{1}{2}\log\left|\sum_{k=1}^{K}\mathbf{D}_{k}\mathbf{S}^{-1}\mathbf{D}_{k}\right|-\frac{1}{2}\,\mbox{tr}\left[\left(\sum_{k=1}^{K}\mathbf{D}_{k}\mathbf{S}^{-1}\mathbf{D}_{k}\right)\mathbf{S}\right]$
(7)

We utilize the properties of the trace and element-wise product operators to
obtain

$\,\mbox{tr}(\mathbf{D}_{k}\mathbf{S}^{-1}\mathbf{D}_{k}\mathbf{S})=\mathbf{d}_{k}^{T}\left(\mathbf{S}^{-1}\odot\mathbf{S}\right)\mathbf{d}_{k}$
(8)

where $\odot$ is element-wise or Hadamard product.

Combing the penalties and the log-likelihood function we have

$\mathcal{F}^{\star}(\mathbf{d}_{1},\ldots,\mathbf{d}_{K})=\mathcal{L}^{\star}(\mathbf{d}_{1},\ldots,\mathbf{d}_{K})+\gamma\left(\sum_{k=1}^{K}\mathbf{d}_{k}^{\intercal}\mathbf{d}_{k}-m\right)-\lambda\log\left|\mathbf{D}^{\intercal}\mathbf{D}\right|$
(9)

for $\gamma,\lambda\in\mathbb{R}^{+}$. This optimization problem can be solved
numerically utilizing gradient descent. Gradient descent is an iterative,
first-order optimization method that finds local minima or maxima of the
differentiable function. We propose to estimate the columns of $\mathbf{D}$
iteratively conditioning on the values of all other column vectors.

To compute the partial derivatives of (9) with respect to $\mathbf{d}_{k}$, we
note that $\mathbf{D}_{k}=\,\mbox{diag}(\mathbf{d}_{k})$ and
$\partial\,\mbox{diag}(\mathbf{d}_{k})=\,\mbox{diag}(\partial\mathbf{d}_{k})$.
The partial derivative is

$\frac{\partial\mathcal{F}^{\star}}{\partial\mathbf{d}_{k}}=\left[\mathbf{T}\odot\mathbf{S}^{-1}\right]\mathbf{d}_{k}-\left(\mathbf{R}^{-1}\odot\mathbf{R}\right)\mathbf{d}_{k}+2\gamma\mathbf{d}_{k}-2\lambda\left[\mathbf{D}(\mathbf{D}^{\intercal}\mathbf{D})^{-1}\right]_{k}$
(10)

where
$\mathbf{T}=(\sum_{k=1}^{K}\mathbf{D}_{k}\mathbf{S}^{-1}\mathbf{D}_{k})^{-1}$
and $\left[\mathbf{A}\right]_{k}$ denotes the $k$-the column of the matrix
$\mathbf{A}$. To maximize the $\mathcal{F}^{\star}$, we take a step in the
direction of a positive gradient. Let us denote $\omega\in\mathbb{R}$ to be a
learning rate, then the proposed update for $\mathbf{d}_{k}$ is defined as:

$\mathbf{d}_{k}^{\ (t+1)}=\mathbf{d}_{k}^{\
(t)}+\omega\frac{\partial\mathcal{F}^{\star}}{\partial\mathbf{d}_{k}}.$ (11)

Through some experimentation, setting $\omega=0.1$ worked well in all the
simulations.

###  2.4 Hierarchical Clustering Approach

Hierarchical clustering refers to a widely utilized family of clustering
algorithms (Murtagh and Legendre, 2014; Beeferman and Berger, 2000; Erica et
al., 2018), that are based on an iterative procedure of either merging or
splitting nested clusters. Merging or splitting is also known as bottom-up and
top-down approaches correspondingly. Many readily available implementations
and generalizations of the algorithm ensure robustness on various input data
configurations (Balcan et al., 2014).

Here, we leverage the bottom-up hierarchical approach known as agglomerative
clustering to find the column clusters. However, instead of clustering the
transpose of the dataset, we propose to use the absolute values of
correlations as a feature representations for each variable. The absolute
value will remove the effect of positive and negative correlations. Then we
apply agglomerative clustering to the feature representations. To obtain the
variables features, if $\mathbf{S}^{\star}$ is theestimate of the correlation
matrix based on dataset and $\mathbf{S}$ is the MLE of the covariance matrix,
we calculate

$\mathbf{R}^{\star}=\left|\mathbf{S}^{\star}\right|=\left|\,\mbox{diag}(\mathbf{S})^{-\frac{1}{2}}\
\mathbf{S}\ \,\mbox{diag}(\mathbf{S})^{-\frac{1}{2}}\right|$ (12)

where $\left|\mathbf{A}\right|$ applies the absolute value element-wise to the
matrix $\mathbf{A}$ and $\,\mbox{diag}(\mathbf{A})$ creates a diagonal matrix
using the diagonal elements of the matrix $\mathbf{A}$. Treating
$\mathbf{R}^{\star}$ as feature vectors, we then use the $l_{2}$-norm to
construct a distance matrix. This distance results in variables with similar
correlation vectors to be close to each other.

The agglomerative clustering procedure starts by assigning every variable into
its own group. On every iteration, it merges two of the most similar groups,
as measured by a selected similarity/dissimilarity metric, which is known as
linkage criteria. The linkage criterion is a metric that defines the distance
between two clusters. There are three most widely used linkage criteria
(Vijaya et al., 2019): single linkage, complete linkage and average linkage.
single, complete and average linkage defines distance between two clusters as
the minimum, maximum or average distance among distances between the two
clusters, respectively. We propose utilizing average linkage criteria as this
method results in clusters with the highest cohesion as discussed in Sokal and
Michener (1958). Once the $K$ variable groups have been obtained, we construct
the corresponding block diagonal covariance matrix.

###  2.5 Comparison & Experiments

In this subsection, we aim to conduct experiments with simulated data to
establish which of the three proposed approaches for estimating the block-
diagonal matrix is the most accurate and efficient. Section 2.5.1 compares the
execution time and accuracy of the three approaches while varying the number
of observations in the dataset. We find that the hieracherical approach is the
the most accurate and efficient one.

The we compare the hierarchical approach to other known covariance estimators.
In the 2.5.2 we experiment with recovering two different structures of the
covariance matrix with our vs state-of-the-art approaches, including the
shrinkage estimator (Ledoit and Wolf, 2004), Graphical Lasso (Friedman et al.,
2008), factor-analyzer UCUU (Tu and Subedi, 2022) and MLE. We find that
methods that artificially enforce block-diagonal covariance structure are more
accurate and the proposed method is indeed more robust than the UCUU model (Tu
and Subedi, 2022) for matrices with negative covariances.

Finally, we investigate the limitations of the hierarchical estimator by
increasing the complexity the block-diagonal matrix $\bm{\Sigma}$ when the
number of columns group is known and unknown. We examine using the maximum
silhouette score (Rousseeuw, 1987) to determine the number of groups. We
discovered that the usage of the hierarchical algorithm is limited when
estimating high number of small blocks.

####  2.5.1 Accuracy and Efficiency of the three Proposed Estimation
Approaches

In this subsection, we investigate the accuracy and efficiency of the three
proposed approaches to the problem of recovering column clusters
$\mathbf{D}_{k}$. The experiment setup is similar to other proposed in Tu and
Subedi (2022). In this subsection, we generate $N$ observations from $p=8$
dimensional Guassian distribution with
$\bm{\mu}=(0,1,2,3,4,5,6,7)^{\intercal}$ and two different covariance
matrices. Figure 1 displays two options for the covariance matrix that are
considered in the experiments: $\bm{\Sigma}_{A}$ with only non-negative off-
diagonal entries; $\bm{\Sigma}_{B}$ with negative off-diagonal entries.
$\bm{\Sigma}_{A}$ and $\bm{\Sigma}_{B}$ are both block diagonal matrices
displayed in Figure 1. For case A, we have
$\bm{\Sigma}_{A}=\,\mbox{diag}(\bm{\Sigma}_{A1},\bm{\Sigma}_{A2},\bm{\Sigma}_{A3})$
where

$\bm{\Sigma}_{A1}=\left[\begin{array}[]{ccc}4.5&2&2\\\ 2&4.5&2\\\
2&2&4.5\end{array}\right],\quad\bm{\Sigma}_{A2}=\left[\begin{array}[]{ccc}4.5&2&2\\\
2&4.5&2\\\
2&2&4.5\end{array}\right],\quad\bm{\Sigma}_{A3}=\left[\begin{array}[]{cc}4.5&2\\\
2&4.5\\\ \end{array}\right],$

and for case B, we have
$\bm{\Sigma}_{B}=\,\mbox{diag}(\bm{\Sigma}_{B1},\bm{\Sigma}_{B2},\bm{\Sigma}_{B3})$
where

$\bm{\Sigma}_{B1}=\left[\begin{array}[]{ccc}4.5&-2&1\\\ -2&4.5&2\\\
1&2&4.5\end{array}\right],\quad\bm{\Sigma}_{B2}=\left[\begin{array}[]{ccc}4.5&-2&2\\\
-2&4.5&2\\\
2&2&4.5\end{array}\right],\quad\bm{\Sigma}_{B3}=\left[\begin{array}[]{cc}3&2\\\
2&4.5\\\ \end{array}\right].$

The choice of covariance matrices is motivated by restrictions for off-
diagonal elements to be strictly non-negative utilized commonly in the
literature (Wong et al., 2017; Tu and Subedi, 2022; Martella et al., 2008). To
goal here is to recover a covariance matrix of the general structure,
including the possibility of negative off-diagonal elements. Note that both
matrices in Figure 1 have a natural block-diagonal structure, which provides
us with ground truth for column cluster indicator matrices $\mathbf{D}_{k}$ as
shown previously in (3).

Figure 2 shows results with data generated from
$\mathcal{N}(\bm{\mu},\bm{\Sigma}_{A})$. We vary the number of data points,
$N=50,80,100,200,500,800,1600$, to understand the impact on both accuracy and
execution time of the proposed methods. For every size $N$, we run the
experiment $200$ times and set number of column clusters to three. The
accuracy is calculated as the total number of times when an algorithm returned
correct column groups $\mathbf{D}$ divided by the total number of attempts. 2
(right) demonstrates the average execution time in milliseconds over $200$
runs together with it’s standard deviation.

From Figure 2, it is clear that the hierarchical estimator for the column
groups $\mathbf{D}_{k}$ is both accurate and efficient. It gets to $100$%
starting from $100$ data points in the training set and demonstrates solid
accuracy $>90\%$ even with $N=50$. The numerical estimator benefits from the
increasing dataset size. However, the time it needs to converge grows as well.
In this experiment, the greedy estimator proved to be the least efficient and
the least accurate for larger sample sizes.

Figure 2:  Comparison of accuracy & execution time over $200$ runs for
estimating $D_{k}$ for block-diagonal $8\times 8$ covariance matrix with $3$
blocks defined by $\bm{\Sigma}_{A}$ which has only positive entries. (Left)
The accuracy of the proposed approaches is given a number of the data points
in the generated dataset $X$. (Right) Execution time of the proposed
approaches.

Figure 3:  Comparison of accuracy & execution time for estimating $D_{k}$ for
block-diagonal $8\times 8$ covariance matrix with $3$ blocks defined by
$\bm{\Sigma}_{B}$ which has some negative entries. (Left) The accuracy of the
proposed approaches given a number of the data points in the generated dataset
$X$. (Right) Execution time of the proposed approaches.

Figure 3 demonstrates the results of the analogous experiment with data
generated from $\mathcal{N}(\bm{\mu},\bm{\Sigma}_{B})$. Interestingly, the
execution time of both numerical and greedy estimators increases. While the
numerical estimator demonstrates degraded accuracy, the greedy one improves in
the presence of negative covariances or a greater variety of values within a
covariance block. At the same time, the hierarchical approach leads in terms
of both efficiency and accuracy, and it can correctly cluster columns even for
the matrix with negative covariances. Therefore, for the following
experiments, we proceed with the hierarchical approach.

####  2.5.2 Recovering Block-Diagonal Covariance Matrix & Comparison with the
Literature

In this section, we evaluate the accuracy of estimation of the block-diagonal
covariance matrix $\bm{\Sigma}$ using the mean absolute percentage error
(MAPE). For estimating $\bm{\Sigma}$ with an estimate $\hat{\bm{\Sigma}}$, we
compute MAPE % measure as following:

$100\times\frac{\sum_{i,j}|\sigma_{ij}-\hat{\sigma}_{ij}|}{\sum_{i,j}\sigma_{ij}}.$
(13)

We compare the hierarchical approach to the MLE and three other methods
utilized in the literature: shrinkage estimator (Ledoit and Wolf, 2004),
Graphical Lasso (Friedman et al., 2008) and UCUU model (Tu and Subedi, 2022).

Figure 4: Average estimations over ten runs of recovering block-diagonal
covariance matrix with only positive (top) / negative (bottom) off-diagonal
entries with five methods including MLE (empirical), Ledoit-Wolf, Graphical
Lasso, Hierarchical (proposed) & UCUU model by Tu and Subedi (2021). For
standard deviations, see A.

The first column of Figure 4 shows two block-diagonal covariance matrices used
to generate data-sets. The top row has a block-diagonal covariance matrix with
$4.5$ on the diagonal and $2$ on the off-diagonals within a block. The bottom
row has a block-diagonal covariance matrix with $4.5$ on the diagonal and $-1$
on the off-diagonals within a block. We use $N=300$, $p=12$ and
$\bm{\mu}=(1,2,\ldots,12)$. Every estimation of the covariance matrix is
replicated ten times with the average $\hat{\Sigma}$ for each estimator
depicted in Figure 4. For the hierarchical approach and UCUU model (Tu and
Subedi, 2022) we set the number of column groups to three.

Our analysis show that artificially enforcing block-diagonal covariance
structure is justified for the problems when such covariance structure is
expected. Note that hierarchical approach and UCUU model are at least two
times more accurate (in terms of MAPE) than others when the off-diagonal
elements in the blocks are positive. On the other hand, the proposed
hierarchical method is robust in the presence negative off-diagonal
covariances, while the performance of the UCUU model by Tu and Subedi (2022)
suffers dramatically as it is not designed for such case.

Note that, MAPE is only helpful for relative algorithm accuracy comparison but
not for absolute estimation quality assessment. It allows the comparisons of
the form “method X is 2 times more accurate than Y”. However, we do not deduce
that all methods are highly inaccurate for the covariance matrix with negative
off-diagonal entries because the denominator of equation 13 is small.
Therefore even small errors will result in large percentages.

###  2.6 Impact of the Covariance Block Structure on the Performance of
Hierarchical Estimator

In this section, we investigate the limitations of the hierarchical estimator
by increasing the complexity the block-diagonal matrix $\bm{\Sigma}$ when the
number of columns group is known and unknown. We measure the accuracy of the
estimated column groups compared the generated column groups. The complexity
of the variable groups $\mathbf{D}$ is influenced by the dimensionality of
number of features variables (i.e. columns) in the dataset, $p$; and the
number of blocks or column groups, $K$. We fix the number of rows in the
dataset to $N=50$ and vary $p$ and $K$ to understand its impact on the
estimator’s accuracy. This choice is motivated by common real-life scenarios
(e.g. genomics, natural language processing etc.) when the number of
observations is often limited compared to the number of potential features.

We evaluate the accuracy of the hierarchical estimator with three different
levels for the number of variables, $p$, and the number of blocks or variable
groups, $K$. They will be referenced as low, medium and high and summarized in
Table 1. The number of column groups at every level is proportional to number
of variables. For example, for the low-low setting $(p,K)=(24,3)$ while the
medium-low setting has $(p,K)=(96,4)$. For simplicity, the number of variables
within each block is set to $p/K$.

$(p,K)$ |  # of variables - $p$
---|---
# of column groups - $K$ | low | medium | high
low | $(24,3)$ | $(96,4)$ | $(384,8)$
medium | $(24,6)$ | $(96,8)$ | $(384,16)$
high | $(24,12)$ | $(96,16)$ | $(384,32)$
Table 1: Summary of different combinations of the number of variables, $p$,
and the number of variable groups, $K$, in $\bm{\Sigma}$. We set $K$
dynamically in relation to $p$ (e.g. low number of blocks is $3$ for the data
with $24$ variables and $12$ for the high-dimensional data with $384$ columns)
Figure 5: Examples of the covariance matrices $\bm{\Sigma}$ used to generate
data for the simulation study. These show when the number of columns is set to
medium $p=96$ and the number of blocks, $k$ varies from 6 (low), 8 (medium)
and 16 (high). Note the added noise.

For every combination of $p$ and $K$, we run the analysis $500$ times. In
every run, we sample the elements of $\bm{\mu}$ randomly from the uniform
$[0,1]$ distribution. For the block diagonal covariance matrix, each block,
$\bm{\Sigma}_{k}$, is set to $\mathbf{A}^{\intercal}\mathbf{A}$ where
$\mathbf{A}$ is $p/K$-dimensional square matrix and the entries are generated
from uniform $U(1,2)$. Once each block has been constructed, we add random
noise to the overall covariance matrix $\bm{\Sigma}$ to bring this simulation
closer to real-world scenarios. The random noise is equal to
$0.5\mathbf{E}^{\intercal}\mathbf{E}$ where $\mathbf{E}$ is $p$-dimensional
square matrix and the entries are generated from uniform $U(0,1)$.

Figure 6 demonstrates the accuracy of the hierarchical estimator for different
problem complexities in two scenarios: known vs an unknown actual number of
column clusters (blocks of the covariance matrix). In the second case, we pick
the number of groups that produces the maximum silhouette score (Rousseeuw,
1987) with the search range $K\pm 2\sqrt{p}$, where $K$ is the number of
clusters and $p$ is a dimensionality of the data.

Figure 6: The average accuracy of the hierarchical estimator for the number of
column groups $\mathbf{D}$ while number of variables and the number column
groups different complexities defined by the dimensionality of the
dataset(low-high) and the true number of blocks in the covariance matrix
$\bm{\Sigma}$ (number of column clusters). Tested for the scenario of known
(left) and unknown (right) $K$ clusters

Figure 6 shows that the proposed estimator is accurate for the covariance
structures in the low- and medium- setting data but has difficulty in the
high-dimensional ($p=384$) setting. Although this might be due the sample size
to be fixed at $50$). We see that accuracies for the known vs unknown number
of blocks are similar, which signifies the usefulness of the estimator even in
the case when the number of column clusters is unknown.

##  3 The Bi-Clustering Model with general block diagonal structure

This section extends the proposed block diagonal structure and estimation
procedure for model-based bi-clustering. We start by defining a probability
density function for a finite $G$-component finite Gaussian mixture model
(GMM),

$g(\mathbf{x}_{i}\mid\bm{\theta})=\sum_{g=1}^{G}\pi_{g}f\left(\mathbf{x}_{i}\mid\bm{\mu}_{g},\bm{\Sigma}_{g}\right)$
(14)

where $\pi_{g}>0$, $\sum_{g=1}^{G}\pi_{g}=1$ and $\bm{\theta}$ represents all
the model parameters. For bi-clustering model, we assume that
$\bm{\Sigma}_{g}$ has a block structure as given (2). The biggest difference
between the proposed model and existing models based on the factor-analyzer
model is in how the covariance matrices are constructed.

There is no closed-form solution for estimating the parameters of the Gaussian
mixture model. Parameters of such models are usually estimated by leveraging
the expectation-maximization (EM; Dempster et al., 1977) algorithm, or its
variants. The EM algorithm is an iterative algorithm that finds local maximum
likelihood estimates of the model parameters $\bm{\theta}$ when the data are
incomplete or are treated as incomplete. In the context of Gaussian mixture
model, the sample of $N$ observations from (14) are the observed data. We
define the component indicator variable $\mathbf{Z}$ for the row clusters such
that $Z_{ig}=1$ if observation $i$ belongs to group $g$ and 0 otherwise. In
the clustering context, the group memberships are unknown and thus, the
component indicator variable $Z_{ig}$ are treated as the unobserved latent
variable.

The EM algorithm alternates between the expectation (E-step) and maximization
(M-step) steps. In the E-step, we compute the expected value of the complete
data log-likelihood function using the current parameter values
$\bm{\theta}^{(t)}$. In the M-step, the algorithm finds the estimate of
$\bm{\theta}$ that maximizes the expected value of the log-likelihood from the
E-step. For a $G$-component GMM, the expected value of the complete data log-
likelihood function is

$\sum_{i=1}^{N}\sum_{g=1}^{G}\hat{z}_{ig}[\log\pi_{g}-\frac{1}{2}\log|\bm{\Sigma}_{g}|-\frac{1}{2}\,\mbox{tr}\left[\bm{\Sigma}_{g}^{-1}(\mathbf{x}_{i}-\bm{\mu}_{g})(\mathbf{x}_{i}-\bm{\mu}_{g})^{\intercal}\right],$

where expected values of the latent variable $z_{is}$ are

$\hat{z}_{ig}=P(Z_{ig}=1|\mathbf{x}_{i},\bm{\theta})=\frac{\pi_{g}f(\mathbf{x}_{i}|\bm{\mu}_{g},\bm{\Sigma}_{g})}{\sum_{h=1}^{G}\pi_{h}f(\mathbf{x}_{i}\mid\bm{\mu}_{h},\bm{\Sigma}_{h})}.$
(15)

For the M-step, the updates for the $\bm{\theta}$ are:

$\hat{\pi}_{g}=\frac{1}{N}\sum_{i=1}^{N}\hat{z}_{ig},\quad\hat{\bm{\mu}}_{g}=\frac{\sum_{i=1}^{N}\hat{z}_{ig}\mathbf{x}_{i}}{\sum_{i=1}^{N}\hat{z}_{ig}}.$
(16)

The maximum likelihood estimate of unconstrained $\bm{\Sigma}_{g}$ is given
by:

$\hat{\bm{\Sigma}}_{g}=\mathbf{S}_{g}=\frac{\sum_{i=1}^{N}\hat{z}_{ig}(\mathbf{x}_{i}-\bm{\mu}_{g})(\mathbf{x}_{i}-\bm{\mu}_{g})^{\intercal}}{\sum_{i=1}^{N}\hat{z}_{ig}}.$
(17)

To estimate a block-diagonal $\bm{\Sigma}_{g}$, we input the $\mathbf{S}_{g}$
to one of the three proposed column clustering estimators (e.g. greedy,
hierarchical, numerical). Then we compute the estimate of the block covariance
matrix for the $g^{th}$ component with $K$ blocks as:

$\sum_{k=1}^{K}\mathbf{D}_{kg}\mathbf{S}_{g}\mathbf{D}_{kg}.$ (18)

For initialization, we use the K-means algorithm (Forgy, 1965) to provide
initial assignment of the latent variables $z_{ig}$. We use Aiken’s
acceleration based criteria (Böhning et al., 1994) to determine convergence of
the EM algorithm using a tolerance of $10^{-4}$. To obtain row clustering
memberships, we take the argmax over the set
$\\{\hat{z}_{i1},\ldots,\hat{z}_{iG}\\}$ and the column groups are obtain by
examining at the matrices $\mathbf{D}_{k}$.

###  3.1 Estimating the Number of Row Clusters

Here, we study the estimation of the number of clusters by the proposed model
and it’s accuracy of selecting the number of row clusters given data generated
from a GMM with $G=3$ components. We consider two scenarios where generate
data from $\mathcal{N}(\bm{\mu}_{g},\bm{\Sigma}_{g})$ where the covariance
matrix has a block diagonal structure and a general structure (or non-block
diagonal structure).

We choose similar parameter values to the simulation from (Tu and Subedi,
2022). In the first scenario, we set

$\begin{split}\bm{\mu}_{1}&=\left(-5,-4,-3,-2,-1,0,1,2\right)^{\intercal},\\\
\bm{\mu}_{2}&=\left(0,1,2,3,4,5,6,7\right)^{\intercal},\\\
\bm{\mu}_{3}&=\left(5,6,7,8,9,10,11,12\right)^{\intercal}\end{split}$

and the covariance matrices are block diagonal. In particular we set
$\bm{\Sigma}_{1}=\,\mbox{diag}(\bm{\Sigma}_{11},\bm{\Sigma}_{12},\bm{\Sigma}_{13})$
where

$\bm{\Sigma}_{11}=\left[\begin{array}[]{ccc}2.5&0.5&0.5\\\ 0.5&3.5&0.5\\\
0.5&0.5&4.5\end{array}\right],\quad\bm{\Sigma}_{12}=\left[\begin{array}[]{ccc}2&1&1\\\
1&2&1\\\
1&1&2\end{array}\right],\quad\bm{\Sigma}_{13}=\left[\begin{array}[]{cc}3.5&3.0\\\
3.0&3.9\\\ \end{array}\right],$

and set
$\bm{\Sigma}_{2}=\,\mbox{diag}(\bm{\Sigma}_{21},\bm{\Sigma}_{22},\bm{\Sigma}_{23})$
where

$\bm{\Sigma}_{21}=\left[\begin{array}[]{ccc}4.5&-2&1\\\ -2&4.5&2\\\
1&2&4.5\end{array}\right],\quad\bm{\Sigma}_{23}=\left[\begin{array}[]{cc}4&2\\\
2&4\\\
\end{array}\right],\quad\bm{\Sigma}_{22}=\left[\begin{array}[]{ccc}4&3&3\\\
3&3.5&3\\\ 3&3&4\end{array}\right],$

and set $\bm{\Sigma}_{3}=\,\mbox{diag}(\bm{\Sigma}_{31},\bm{\Sigma}_{32})$
where

$\bm{\Sigma}_{31}=\left[\begin{array}[]{ccccc}2.1&2&2&2&2\\\ 2&2.5&2&2&2\\\
2&2&3&2&2\\\ 2&2&2&5&2\\\
2&2&2&2&4\end{array}\right],\quad\bm{\Sigma}_{32}=\left[\begin{array}[]{ccc}2&1&1\\\
1&3.5&1\\\ 1&1&2.5\end{array}\right].$

In the second scenario, we set $\bm{\Sigma}_{g}$ to be different random
positive-definite matrices generated by setting
$\bm{\Sigma}_{g}=\mathbf{A}_{g}^{T}\mathbf{A}_{g}$ where the elements of
$\mathbf{A}_{g}$ are generated independently from the uniform distribution on
the interval $(0.0,1.0)$. This comparison considers what happens when the data
is generated from a model that violates the bi-clustering framework or
assumption of block-diagonal covariance matrix.

In both scenarios, we set $\pi_{1}=\pi_{2}=\pi_{3}=1/3$ and we vary the number
of observations per component. Then we chose the number of row cluster by
iterating over the search space $\hat{G}\in\left\\{1\dots 5\right\\}$ and
applying the BIC (Stoica and Selen, 2004; Schwarz, 1978) criterion to select
the number of components. For both scenarios we replicate the experiment $100$
times. We compare performance of the proposed method to the performance of a
GMM.

Figure 7 shows the average estimated number of components with band showing
plus/minus the standard deviation for bi-clustering and GMM under the two
scenarios. Both panels show that the models select the true number of row
clusters more often as the size of input data grows. However, in the left
panel the bi-clustering model arrives to $100$% accuracy earlier. On the other
hand, the right panel shows that a general covariance structure causes issues
for the bi-clustering model which needs lots variable groups to achieve a good
fit of the data. This provides some evidence that in the presence of a block-
diagonal covariance matrix, the bi-clustering is efficient in determining the
number of row clusters. In addition, Figure 13 in B summarizes the accuracy of
identifying the correct number of groups.

Figure 7: Recovery of the number of groups when the covariance structure is
block diagonal (Left) and a random positive definite matrix (Right). Both
panels display the average and plus/minus the standard deviation for a GMM and
the proposed bi-clustering model when $G=3$ while varying the number of
observations within each component.

###  3.2 Bi-Clustering Applications

In this subsection, we demonstrate the applicability of the proposed
biclustering algorithms on benchmark datasets as well as high dimensional
datasets from bioinformatics and topics modelling. We examine the row
clustering performance of the proposed bi-clustering approach i.e., the
Gaussian mixture model with a block-diagonal covariance matrix estimated by
the hierarchical method. We compare the results of the proposed approach with
$K$-means (Forgy, 1965), the unrestricted Gaussian mixture model (Reynolds,
2009), and two other bi-clustering methods (UCUU by Tu and Subedi (2022) and
spectral co-clustering Dhillon (2001)). We use $K$-means, GMM and spectral co-
clustering implementations with their default parameters from the sklearn
library (Pedregosa et al., 2011). Additionally, we implemented the factor-
analyzer UCUU model by Tu and Subedi (2022) in Python.

####  3.2.1 Bi-clustering benchmark datasets

Here, we applied the algorithms on well-known low-dimensional benchmark
datasets: Wine (Cortez et al., 2009), Olive (Forina and Tiscornia, 1982) and
Ecoli (Nakai and Kanehisa, 1991). We scale the data to have the mean equal to
$0$ and variance equal to $1$ prior to clustering. We compare the clustering
performanace using the adjusted rand index (ARI) (Rand, 1971; Hubert and
Arabie, 1985; Zhang et al., 2012), accuracy and computational time. The
results are summarized in Table 2.

Table 2: Average adjusted rand index (ARI), accuracy and execution time over
10 different runs of all algorithms on the given datasets.

Datasets | Wine | Olive | Ecoli
---|---|---|---
| ARI | %, acc | time | ARI | %, acc | time | ARI | %, acc | time
|  |  | (in sec) |  |  | (in sec) |  |  | (in sec)
K-means | $0.897$ | $96.6\%$ | $0.015$ | $0.448$ | $76.5\%$ | $0.047$ | $0.509$ | $65.2\%$ | $0.036$
GMM | $0.831$ | $92.7\%$ | $\mathbf{0.009}$ | $\mathbf{0.601}$ | $78.1\%$ | $\mathbf{0.019}$ | $0.646$ | $74.4\%$ | 0.015
UCUU Model | $\mathbf{0.948}$ | $\mathbf{98.3}\%$ | $4.320$ | $0.517$ | $79.1\%$ | $19.62$ | $-$ | $-$ | $-$
Proposed | $0.945$ | $98.3\%$ | $0.120$ | $0.574$ | $\mathbf{80.4\%}$ | $0.169$ | $\mathbf{0.656}$ | $\mathbf{76.2}\%$ | $0.138$
Spectral co-clustering | $0.738$ | $90.9\%$ | $0.020$ | $0.237$ | $57.3\%$ | $0.042$ | $0.394$ | $56.4\%$ | $0.048$

The proposed algorithm has similar accuracy compared to known row clustering
algorithms and other bi-clustering approaches and even outperforms others on
the Ecoli dataset. However, is almost $10\times$ slower in terms of execution
time than other listed approaches except for the UCUU model. The UCUU model,
which is most similar algorithm to our proposed method, runs $400\times$
slower than $K$-means, GMM and spectral co-clustering. Thus, the proposed
approach has an advantage in terms of lower execution time compared to the bi-
clustering UCUU model, which is very similar in utility.

####  3.2.2 Bi-clustering in Bioinformatics

Biclustering is used in bioinformatics to simultenously cluster observations
and genes. Identifying groups of highly correlated genes that have different
correlation structure in different groups of individuals (i.e., disease vs
healthy individuals or between individuals with subtypes of diseases) can shed
light into underlying biological mechanisms of disease development.

Here, we compare the algorithms’ performance on the high-dimensional genomics
classification problems using two benchmark gene expression datasets Alon
(Alon et al., 1999) and Golub (Golub et al., 1999) available in the R package
plsgenomics(Boulesteix et al., 2018). We applied our approach on two versions
of the datasets: (i) a subset of 100 genes selected using ANOVA F-test (Feir-
Walsh and Toothaker, 1974) which from here on is referred to as Alon100 and
(ii) full datasets where all available genes are used. The ANOVA F-test
selects the set of $K=100$ features with highest ratio of explained to
unexplained variance between the given ground truth groups. We applied our
approach other competing approaches. Unfortunately, we were unable to run the
UCUU model by Tu and Subedi (2022) on the full set of features because of
$(i)$ increasing computational time; $(ii)$ authors’ note that this algorithm
is not designed for the scenario when the number of features is significantly
higher than the number of observations. Therefore, UCUU model was only applied
to Alon100.

The performance of the proposed algorithm along with comparisons with other
approaches in terms of ARI and % accuracy on Alon100 and
Alon${}_{\text{full}}$ datasets are provided in Table 3.

Table 3: ARI and accuracy of different algorithms on given high-dimensional datasets and their preprocessed versions with $100/150$ features selected with ANOVA F-testFeir-Walsh and Toothaker (1974).  Datasets |  Alon100 |  Golub150 |  Alon${}_{\text{full}}$ |  Golub${}_{\text{full}}$
---|---|---|---|---
| ARI | %, acc | ARI | %, acc | ARI | %, acc | ARI | %, acc
K-means | $0.510$ | $86.0\%$ | $0.892$ | $97.3\%$ | $-0.002$ | $55.3\%$ | $0.162$ | $70.8\%$
GMM | $0.545$ | $87.0\%$ | $\mathbf{0.935}$ | $\mathbf{98.4}\%$ | $0.001$ | $56.4\%$ | $0.161$ | $53.2\%$
UCUU Model | $\mathbf{0.609}$ | $\mathbf{89.2}\%$ | $0.892$ | $97.3\%$ | $-$ | $-$ | $-$ | $-$
Proposed | $0.541$ | $87.0\%$ | $0.892$ | $97.3\%$ | $-0.006$ | $54.8\%$ | $\mathbf{0.212}$ | $\mathbf{73.6\%}$
Spectral co-clustering | $0.541$ | $87.0\%$ | $0.892$ | $97.3\%$ | $\mathbf{0.040}$ | $\mathbf{62.9\%}$ | $0.185$ | $72.2\%$

In the case of Alon100, the UCUU model had a slightly higher ARI compared to
the proposed model. However, our proposed method takes substantially less time
to fit. Unfortunately, the performance of all methods degraded significantly
on the full Alon dataset, with the Spector co-clustering method achieving
maximum accuracy of $62.9$%.

In Figure 8, we visualize the true and estimated correlation structures in the
two estimated clusters of the Alon100 dataset. Note that column permutations
are done to achieve a checkerboard-like structure as illustrated. As can be
seen in Figure 8, the estimated column clusters are grouping variables that
have high correlation (either positive or negative). This can provide valuable
biological insight. For example, the variables X53586 (mRNA for integrin alpha
6), X54942 (ckshs2 - mRNA for Cks1 protein homologue) and X54941 (ckshs1 -
mRNA for Cks1 protein homologue) are all assigned to column cluster 5 in Group
1 (with majority of tumour samples) indicating that they have a high
correlation in tumour samples. Both integrin alpha 6 (Beaulieu, 2020) and Cks1
protein (Shapira et al., 2005) have previously been identified to play roles
in colorectal cancer. On the other hand, in Group 2 (with majority of healthy
samples), these variables were assigned to different column clusters: X53586
was assigned to column cluster 1, X54942 was assigned to column cluster 5 and
X54941 was assiged to column cluster 3, indicating that these variables have
little to no correlation in healthy samples.

(a) Observed correlation of predicted Group 1

(b) Estimated correlation of predicted Group 1

(c) Observed correlation of predicted Group 2

(d) Estimated correlation of predicted Group 2

Figure 8: True and estimated correlation matrices of the predicted groups from
the Alon100 dataset. Note that most tumour samples were classified as Group 1
and most healthy samples were classified as Group 2.

Golub dataset (Golub et al., 1999) comprises 2030 gene expressions of 72
patients, 47 of them are diagnosed with acute lymphoblastic leukemia (ALL) and
25 with acute myeloid leukemia (AML). Similar to the Alon dataset, we
construct Golub150 with subset of 150 features selected with ANOVA F-test
(Feir-Walsh and Toothaker, 1974). We compare the performance of the proposed
approach with competing approach on both the Golub150 and the full dataset.
Here, for the Golub150 dataset, the GMM has the largest ARI slightly higher
compared to our proposed model. However, due to the block diagonal covariance
structure, our proposed method is more parsimonious than GMM. Moreover, the
column clusters can provide valuable insights as illustrated below which the
GMM lacks. Similar to the full Alon data, the performance of all methods
degraded significantly on the full Golub dataset - in this case substantially
for the GMM, while our proposed method achieves maximum accuracy of $73.6$%
among the competiting models.

(a) Observed correlation of predicted Group 1

(b) Estimated correlation of predicted Group 1

(c) Observed correlation of predicted Group 2

(d) Estimated correlation of predicted Group 2

Figure 9: True and estimated correlation matrices of the predicted groups from
the Golub150 dataset. Note that all 11 ALL samples were classified as Group 1
and 26 out of 27 AML samples were classified as Group 2.

In Figure 9, we visualize the observed and estimated correlation structures in
the two predicted clusters of the Golub150 dataset. Column cluster 3 in Group
1 (i.e., the group comprising of all ALL samples) consists of 27 variables
which have high correlations (either positive or negative) in the ALL samples.
Out of the 27 variables, 24 of them were assigned to column cluster 1 in Group
2 (the group comprising of mostly AML samples). As evident from Figure 9, the
magnitude of the correlations among the variables in column cluster 1 of Group
2 is lower compared to that of column cluster 3 of Group 1. Interestingly,
pathway analysis using Reactome (Gillespie et al., 2022) identifies that 9 out
of the above 24 genes are known to be involved in pathways related to platelet
activation, signaling and aggregation. Platelets are specialized blood cells
which are primarily responsible for preventing bleeding but are also recently
found to promote tumor growth and metastasis (Yan and Jurasz, 2016).
Furthermore, cancer cells have been shown to induce platelet activation and
aggregation in several cancers, a phenomenon known as tumor cell–induced
platelet aggregation (Jurasz et al., 2004). Investigating the differences in
correlation structures among these genes between AML and ALL samples could
provide valuable biological insight.

####  3.2.3 Bi-clustering in Topic Modelling

Topic Modelling Wallach (2006) is an area of Natural Language Processing that
tries to find abstract human-interpretable topics in a collection of an
arbitrary number of documents. Essentially, the problem of grouping documents
into unsupervised topics can be solved by leveraging clustering techniques.
Human language is incredibly complex. There are multiple possible groupings
one can find in a set of documents, some having a stronger signals-no-noise
ratio than others. e.g., Trip advisor’s hotel reviews can be grouped by review
sentiment (recommend/do not recommend), hotel geographical location or more
generalized hotel features (has access to the beach, superb breakfast, polite
managers etc.). Each of the clustering can be desired depending on the
application. However, classical clustering algorithms will require additional
data post-processing to provide interpretations of the topics they produce
Sievert and Shirley (2014).

In this subsection, we apply the proposed bi-clustering algorithm to the topic
modelling problem. We test our algorithm on two real-world datasets of
customer reviews available on Kaggle. The first is Trip Advisor hotel Alam et
al. (2016) data that consists of 20491 unique free-text hotel reviews
accompanied by a numerical rating (from 1 to 5). The second one includes
reviews and ratings of three Disneyland (Chillar, 2022) branches. We use only
a subset of the data that reviews Disneyland Paris (13630 free-text
responses). The language of reviews in both datasets is English. There are no
reference labels available in this application. We aim to demonstrate that the
proposed algorithm as a tool for summarizing customer reviews data.

For better explainability of the resulted topics and column groupings, we
decided to apply the bi-clustering algorithm to the document-terms matrix
instead of document-embeddings Grootendorst (2022) data type. Testing the
proposed algorithm with semantic text embeddings is left for future work.

The entries of the document-terms matrix are term frequency-inverse document
frequency scores (TF-IDF) Aizawa (2003). Let us define $\text{tf}(t,d)$ as a
relative frequency of the term $t$ in some document $d$. Then, the inverse
document frequency of the term $t$ given a collection $D$ of $N$ documents is
computed as $\text{idf(t, D)}=\log\frac{N}{|{d\in D:t\in d}|}$.

To simplify, the term frequency measures how frequent a term is for a given
document. The inverse document frequency measures how rare this term is in the
overall collection of documents. We define a TF-IDF score of a term $t$ int he
document $d$ as $\text{tf-idf}(t,d)=\text{tf}(t,d)\text{idf(t, D)}$. The score
is maximized when a specific term is frequent in a given document but uncommon
for other documents. The TF-IDF score of the stop words is close to zero, as,
by definition, stop words have a high probability of being present in every
document in the collection.

Large collections of documents may include over a thousand unique words.
Therefore, before fitting the bi-clustering model, we select only $1000$
features with the highest average $TF$ scores. Additionally, we utilize
SelectKBest from the sklearn Pedregosa et al. (2011) Python library to select
$400$ terms that are the most predictive of the review rating.

Figure 10 and 11 demonstrate the column and row groupings for the Trip Advisor
and Disneyland Paris data. The heatmap entries are computed as the average TF-
IDF score of terms inside the column grouping $A$ row cluster $B$ divided by
average TF-IDF scores of $A$ among all row clusters and multiplied by $100$%.
This can be interpreted as how much ”more important” a column cluster $A$ is
for a row cluster $B$ than average.

From Figure 10, we infer that lower hotel ratings may be attributed to the
cigarette smell, rude managers, dirty conditions, problems with the check desk
and unexpected credit card charges. The highest ratings are correlated with
access to public transport or the hotel’s location within walking distance of
the city’s landmarks, friendly hotel staff and beautiful decorations of the
premises. Figure 11 demonstrates that both feature clusters of food (ice cream
and chips-burgers-fries) are highly associated with the lowest ratings of
Disneyland Paris. It may signal that the food quality is lower than expected.
At the same time, people are unhappy with closed attractions, rude staff,
pushing and perhaps the fact that most people in Paris speak exclusively
French.

Figure 10: Resulting column and row clusters from the Trip Advisor hotel
reviews. Entries of the heatmap represent how much “more important” a given
column cluster is for a given row cluster than average.  Figure 11: Resulting
column and row clusters from the Disneyland Paris reviews. Entries of the
heatmap represent how much “more important” a given column cluster is for a
given row cluster than average.

##  4 Conclusions

In this work, we have proposed a bi-clustering algorithm based on a finite
mixture of multivariate Gaussian models with block diagonal covariance
structure. We demonstrated that our approach has a comparable clustering
performance to the state-of-the-art algorithms in the field while providing
substantial improvement in the computational time. We also illustrate the
applicability of the proposed approach on real world datasets from various
fields.

## References

  * (1)
  * Aizawa (2003) Aizawa, A. (2003), ‘An information-theoretic perspective of tf–idf measures’, Information Processing & Management 39(1), 45–65.
  * Alam et al. (2016) Alam, M. H., Ryu, W.-J. and Lee, S. (2016), ‘Joint multi-grain topic sentiment: modeling semantic aspects for online reviews’, Information Sciences (339), 206–223.
  * Alon et al. (1999) Alon, U., Barkai, N., Notterman, D., Gish, K., Ybarra, S., Mack, D. and Levine, A. (1999), ‘Alon u, barkai n, notterman da, gish k, ybarra s, mack d, levine ajbroad patterns of gene expression revealed by clustering analysis of tumor and normal colon tissues probed by oligonucleotide arrays. proc natl acad sci usa 96: 6745-6750’, Proceedings of the National Academy of Sciences of the United States of America 96, 6745–50.
  * Balcan et al. (2014) Balcan, M.-F., Liang, Y. and Gupta, P. (2014), ‘Robust hierarchical clustering’, The Journal of Machine Learning Research 15(1), 3831–3871.
  * Beaulieu (2020) Beaulieu, J.-F. (2020), ‘Integrin $\alpha$6$\beta$4 in colorectal cancer: expression, regulation, functional alterations and use as a biomarker’, Cancers 12(1), 41.
  * Beeferman and Berger (2000) Beeferman, D. and Berger, A. (2000), Agglomerative clustering of a search engine query log, in ‘Proceedings of the sixth ACM SIGKDD international conference on Knowledge discovery and data mining’, pp. 407–416.
  * Ben-Dor et al. (2002) Ben-Dor, A., Chor, B., Karp, R. and Yakhini, Z. (2002), Discovering local structure in gene expression data: the order-preserving submatrix problem, in ‘Proceedings of the sixth annual international conference on Computational biology’, pp. 49–57.
  * Bergmann et al. (2003) Bergmann, S., Ihmels, J. and Barkai, N. (2003), ‘Iterative signature algorithm for the analysis of large-scale gene expression data’, Physical review E 67(3), 031902.
  * Berkhin (2006) Berkhin, P. (2006), A survey of clustering data mining techniques, in ‘Grouping multidimensional data’, Springer, pp. 25–71.
  * Bhatia et al. (2014) Bhatia, P., Iovleff, S. and Govaert, G. (2014), ‘blockcluster: An r package for model based co-clustering’, Journal of Statistical Software 76.
  * Böhning et al. (1994) Böhning, D., Dietz, E., Schaub, R., Schlattmann, P. and Lindsay, B. (1994), ‘The distribution of the likelihood ratio for mixtures of densities from the one-parameter exponential family’, Annals of the Institute of Statistical Mathematics 46, 373–388.
  * Boulesteix et al. (2018) Boulesteix, A.-L., Durif, G., Lambert-Lacroix, S., Peyre, J. and Strimmer., K. (2018), plsgenomics: PLS Analyses for Genomics.  R package version 1.5-2.
https://CRAN.R-project.org/package=plsgenomics

  * Broto et al. (2022) Broto, B., Bachoc, F., Clouvel, L. and Martinez, J.-M. (2022), ‘Block-diagonal covariance estimation and application to the shapley effects in sensitivity analysis’, SIAM/ASA Journal on Uncertainty Quantification 10(1), 379–403.
  * Busygin et al. (2008) Busygin, S., Prokopyev, O. and Pardalos, P. M. (2008), ‘Biclustering in data mining’, Computers & Operations Research 35(9), 2964–2987.
  * Cheng and Church (2000) Cheng, Y. and Church, G. M. (2000), Biclustering of expression data, in ‘Proc Int Conf Intell Syst Mol Biol.’, Vol. 8, pp. 93–103.
  * Chillar (2022) Chillar, A. (2022), ‘Tripadvisor disneyland reviews dataset’.
https://www.kaggle.com/datasets/arushchillar/disneyland-reviews

  * Cortez et al. (2009) Cortez, P., Cerdeira, A., Almeida, F., Matos, T. and Reis, J. (2009), ‘Modeling wine preferences by data mining from physicochemical properties’, Decision support systems 47(4), 547–553.
  * Dempster et al. (1977) Dempster, A. P., Laird, N. M. and Rubin, D. B. (1977), ‘Maximum likelihood from incomplete data via the em algorithm’, Journal of the Royal Statistical Society: Series B (Methodological) 39(1), 1–22.
  * Dhillon (2001) Dhillon, I. S. (2001), Co-clustering documents and words using bipartite spectral graph partitioning, in ‘Proceedings of the seventh ACM SIGKDD international conference on Knowledge discovery and data mining’, KDD ’01, Association for Computing Machinery, New York, NY, USA, pp. 269–274.
  * Diday and Simon (1976) Diday, E. and Simon, J. (1976), Clustering analysis, in ‘Digital pattern recognition’, Springer, pp. 47–94.
  * Erica et al. (2018) Erica, E., Handari, B. and Hertono, G. (2018), Agglomerative clustering and genetic algorithm in portfolio optimization, in ‘AIP Conference Proceedings’, Vol. 2023, AIP Publishing LLC, p. 020217.
  * Feir-Walsh and Toothaker (1974) Feir-Walsh, B. J. and Toothaker, L. E. (1974), ‘An empirical comparison of the anova f-test, normal scores test and kruskal-wallis test under violation of assumptions’, Educational and Psychological Measurement 34(4), 789–799.
  * Forgy (1965) Forgy, E. W. (1965), ‘Cluster analysis of multivariate data: efficiency versus interpretability of classifications’, biometrics 21, 768–769.
  * Forina and Tiscornia (1982) Forina, M. and Tiscornia, E. (1982), ‘Pattern-recognition methods in the prediction of italian olive oil origin by their fatty-acid content’, Annali di Chimica 72(3-4), 143–155.
  * Friedman et al. (2008) Friedman, J., Hastie, T. and Tibshirani, R. (2008), ‘Sparse inverse covariance estimation with the graphical lasso’, Biostatistics 9(3), 432–441.
  * Gallaugher et al. (2022) Gallaugher, M. P. B., Biernacki, C. and McNicholas, P. D. (2022), ‘Parameter-wise co-clustering for high-dimensional data’, Computational Statistics .
  * Gillespie et al. (2022) Gillespie, M., Jassal, B., Stephan, R., Milacic, M., Rothfels, K., Senff-Ribeiro, A., Griss, J., Sevilla, C., Matthews, L., Gong, C. et al. (2022), ‘The reactome pathway knowledgebase 2022’, Nucleic Acids Research 50(D1), D687–D692.
  * Golub et al. (1999) Golub, T. R., Slonim, D. K., Tamayo, P., Huard, C., Gaasenbeek, M., Mesirov, J. P., Coller, H., Loh, M. L., Downing, J. R., Caligiuri, M. A., Bloomfield, C. D. and Lander, E. S. (1999), ‘Molecular classification of cancer: Class discovery and class prediction by gene expression monitoring’, Science 286(5439), 531–537.
  * Grootendorst (2022) Grootendorst, M. (2022), ‘Bertopic: Neural topic modeling with a class-based tf-idf procedure’.
  * Gu and Liu (2008) Gu, J. and Liu, J. S. (2008), ‘Bayesian biclustering of gene expression data’, BMC genomics 9(1), 1–10.
  * Hartigan (1972) Hartigan, J. A. (1972), ‘Direct clustering of a data matrix’, Journal of the american statistical association 67(337), 123–129.
  * Hemmecke et al. (2009) Hemmecke, R., Köppe, M., Lee, J. and Weismantel, R. (2009), Nonlinear integer programming, in ‘50 Years of Integer Programming 1958-2008’, Springer Berlin Heidelberg, pp. 561–618.
  * Hubert and Arabie (1985) Hubert, L. and Arabie, P. (1985), ‘Comparing partitions’, Journal of Classification 2, 193–218.
  * Jurasz et al. (2004) Jurasz, P., Alonso-Escolano, D. and Radomski, M. W. (2004), ‘Platelet–cancer interactions: mechanisms and pharmacology of tumour cell-induced platelet aggregation’, British journal of pharmacology 143(7), 819.
  * Khan et al. (2014) Khan, K., Rehman, S. U., Aziz, K., Fong, S. and Sarasvady, S. (2014), Dbscan: Past, present and future, in ‘The fifth international conference on the applications of digital information and web technologies (ICADIWT 2014)’, IEEE, pp. 232–238.
  * Kluger et al. (2003) Kluger, Y., Basri, R., Chang, J. and Gerstein, M. (2003), ‘Spectral biclustering of microarray data: Coclustering genes and conditions’, Genome research 13, 703–16.
  * Ledoit and Wolf (2004) Ledoit, O. and Wolf, M. (2004), ‘A well-conditioned estimator for large-dimensional covariance matrices’, Journal of Multivariate Analysis 88(2), 365–411.
  * 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 1(1), 24–45.
  * Madeira et al. (2008) Madeira, S. C., Teixeira, M. C., Sa-Correia, I. and Oliveira, A. L. (2008), ‘Identification of regulatory modules in time series gene expression data using a linear time biclustering algorithm’, IEEE/ACM Transactions on Computational Biology and Bioinformatics 7(1), 153–165.
  * Martella et al. (2008) Martella, F., Alfò, M. and Vichi, M. (2008), ‘Biclustering of gene expression data by an extension of mixtures of factor analyzers’, The International Journal of Biostatistics 4.
  * McInnes et al. (2017) McInnes, L., Healy, J. and Astels, S. (2017), ‘hdbscan: Hierarchical density based clustering.’, J. Open Source Softw. 2(11), 205.
  * McNicholas and Murphy (2008) McNicholas, P. D. and Murphy, T. B. (2008), ‘Parsimonious gaussian mixture models’, Statistics and Computing 18, 285–296.
  * Milligan and Cooper (1987) Milligan, G. W. and Cooper, M. C. (1987), ‘Methodology review: Clustering methods’, Applied psychological measurement 11(4), 329–354.
  * Murtagh and Legendre (2014) Murtagh, F. and Legendre, P. (2014), ‘Ward’s hierarchical agglomerative clustering method: which algorithms implement ward’s criterion?’, Journal of classification 31(3), 274–295.
  * Nakai and Kanehisa (1991) Nakai, K. and Kanehisa, M. (1991), ‘Uci machine learning repository’, URL: https://archive. ics. uci. edu/ml/datasets/ecoli .
  * Padilha and Campello (2017) Padilha, V. A. and Campello, R. J. (2017), ‘A systematic comparative evaluation of biclustering techniques’, BMC bioinformatics 18(1), 1–25.
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M. and Duchesnay, E. (2011), ‘Scikit-learn: Machine learning in Python’, Journal of Machine Learning Research 12, 2825–2830.
  * Rand (1971) Rand, W. M. (1971), ‘Objective criteria for the evaluation of clustering methods’, Journal of the American Statistical Association 66(336), 846–850.
  * Reynolds (2009) Reynolds, D. A. (2009), ‘Gaussian mixture models.’, Encyclopedia of biometrics 741(659-663).
  * Rousseeuw (1987) Rousseeuw, P. J. (1987), ‘Silhouettes: A graphical aid to the interpretation and validation of cluster analysis’, Journal of Computational and Applied Mathematics 20, 53–65.
  * Rugeles et al. (2017) Rugeles, D., Zhao, K., Gao, C., Dash, M. and Krishnaswamy, S. (2017), Biclustering: An application of dual topic models, in ‘Proceedings of the 2017 SIAM International Conference on Data Mining’, SIAM, pp. 453–461.
  * Schwarz (1978) Schwarz, G. (1978), ‘Estimating the dimension of a model’, The Annals of Statistics 6(2), 461–464.
  * Shapira et al. (2005) Shapira, M., Ben-Izhak, O., Linn, S., Futerman, B., Minkov, I. and Hershko, D. D. (2005), ‘The prognostic impact of the ubiquitin ligase subunits skp2 and cks1 in colorectal carcinoma’, Cancer 103(7), 1336–1346.
  * Sievert and Shirley (2014) Sievert, C. and Shirley, K. (2014), Ldavis: A method for visualizing and interpreting topics, in ‘Proceedings of the workshop on interactive language learning, visualization, and interfaces’, pp. 63–70.
  * Sokal and Michener (1958) Sokal, R. R. and Michener, C. D. (1958), ‘A statistical method for evaluating systematic relationships’, University of Kansas science bulletin 38, 1409–1438.
  * Stoica and Selen (2004) Stoica, P. and Selen, Y. (2004), ‘Model-order selection: a review of information criterion rules’, IEEE Signal Processing Magazine 21(4), 36–47.
  * Tu and Subedi (2022) Tu, W. and Subedi, S. (2022), ‘A family of mixture models for biclustering’, Statistical Analysis and Data Mining: The ASA Data Science Journal 15, 206 – 224.
  * Vijaya et al. (2019) Vijaya, Sharma, S. and Batra, N. (2019), Comparative study of single linkage, complete linkage, and ward method of agglomerative clustering, in ‘2019 International Conference on Machine Learning, Big Data, Cloud and Parallel Computing (COMITCon)’, pp. 568–573.
  * Wallach (2006) Wallach, H. M. (2006), Topic modeling: beyond bag-of-words, in ‘Proceedings of the 23rd international conference on Machine learning’, pp. 977–984.
  * Wang et al. (2013) Wang, Y. K., Print, C. G. and Crampin, E. J. (2013), ‘Biclustering reveals breast cancer tumour subgroups with common clinical features and improves prediction of disease recurrence’, BMC Genomics 14, 1–15.
  * Wong et al. (2017) Wong, M., Mutch, D. and McNicholas, P. (2017), ‘Two-way learning with one-way supervision for gene expression data’, BMC Bioinformatics 18.
  * Yan and Jurasz (2016) Yan, M. and Jurasz, P. (2016), ‘The role of platelets in the tumor microenvironment: from solid tumors to leukemia’, Biochimica et Biophysica Acta (BBA)-Molecular Cell Research 1863(3), 392–400.
  * Zhang et al. (2012) Zhang, S., Wong, H.-S. and Shen, Y. (2012), ‘Generalized adjusted rand indices for cluster ensembles’, Pattern Recognition 45(6), 2214–2226.

##  Appendix A Standard Deviations for recovering the block-diagonal
covariance matrix

Figure 12: Standard deviations over ten runs of recovering block-diagonal
covariance matrix with only positive (top) / negative (bottom) off-diagonal
entries with five methods including MLE (empirical), Ledoit-Wolf, Graphical
Lasso, Hierarchical (proposed) & UCUU method from Tu, Subebi 2021 Tu and
Subedi (2022)

##  Appendix B Accuracy for selecting the number of components

Figure 13: Accuracy of the number of groups when the covariance structure is
block diagonal (Left) and a random positive definite matrix (Right) for the
GMM and the proposed bi-clustering model when $G=3$ while varying the number
of observations within each component.