# Consistent Biclustering

Cheryl J<EMAIL_ADDRESS>[ Patrick O.
<EMAIL_ADDRESS>[ New York University  Leonard N. Stern
School of Business
New York University
44 West 4th Street
New York, NY 10012-1126

E-mail: e2

###### Abstract

Biclustering, the process of simultaneously clustering observations and
variables, is a popular and effective tool for finding structure in a high-
dimensional dataset. A variety of biclustering algorithms exist, and they have
been applied successfully to data sources ranging from gene expression arrays
to review-website data. Currently, while biclustering appears to work well in
practice, there have been no theoretical guarantees about its performance. We
address this shortcoming with a theorem providing sufficient conditions for
asymptotic consistency when both the number of observations and the number of
variables in the dataset tend to infinity. This theorem applies to a broad
range of data distributions, including Gaussian, Poisson, and Bernoulli. We
demonstrate our results through a simulation study and with examples drawn
from microarray analysis and collaborative filtering.

62H30,

62P10,

62P25,

Biclustering,

block model,

profile likelihood,

consistency,

microarray data,

collaborative filtering,

###### keywords:

[class=AMS]

###### keywords:

name

and

##  1 Introduction

Suppose we are given a data matrix $\boldsymbol{X}=[X_{ij}]$, and our goal is
to cluster the rows and columns of $\boldsymbol{X}$ into meaningful groups.
For example, $X_{ij}$ could be the log activation level of gene $j$ in patient
$i$; our goal is to seek groups of patients with similar genetic profiles,
while at the same time finding groups of genes with similar activation levels.
Alternatively, $X_{ij}$ can indicate whether or not user $i$ reviewed movie
$j$; our goal is to simultaneously cluster the users and the movies. Mirkin
(1996) termed the general clustering process “biclustering,” but it is also
known as direct clustering (Hartigan, 1972), block modeling (Arabie, Boorman
and Levitt, 1978), and co-clustering (Dhillon, 2001).

Empirical results from a broad range of disciplines indicate that biclustering
data is useful in practice. For example, Ungar and Foster (1998) and Hofmann
(1999) found that biclustering helps identify structure in latent class models
for collaborative filtering problems where the data is sparse and diverse
tastes make it difficult to cluster purely based on purchasing or review
habits.

Several biclustering applications exist in the biological sciences. Eisen et
al. (1998) was one of the first papers to note the benefits of clustering
genes and conditions in microarray data, finding that genes with similar
functions cluster together. Recently, Harpaz et al. (2010) applied
biclustering methods to a Food and Drug Adminstration report database,
identifying associations between certain active ingredients and adverse
medical reactions. Several other applications of biclustering exist; see Cheng
and Church (2000),Getz, Levine and Domany (2000), Lazzeroni and Owen (2002),
and Kluger et al. (2003) as well as Madeira and Oliveira (2004) for a
comprehensive survey.

Although their goals are the same, the references above use a variety of
different biclustering algorithms. Clearly, many of these algorithms work well
in practice, but they are often ad-hoc, and there are no rigorous guarantees
as to their performance. In particular, the lack of any notion of consistency
means that practitioners cannot be assured that their discoveries from
biclustering will generalize or be reproducible; collecting more data may lead
to completely different biclusters.

Our first objective in this report is to establish a probabilistic model for
the data matrix, so that biclustering can be formalized as an estimation
problem. Once we have done this, we study a class of biclustering algorithms
based on profile-likelihood, as proposed by Ungar and Foster (1998). We show
that these profile-based procedures are asymptotically consistent as the
dimensions of the matrix tend to infinity, under weak assumptions on the
elements of $\boldsymbol{X}$. Notably, our methods can handle both dense and
sparse data matrices. To our knowledge, this is the first general consistency
result for a biclustering algorithm.

The present treatment was inspired by recent developments in clustering
methods for symmetric binary networks. In that context, $\boldsymbol{X}$ is an
$n$-by-$n$ symmetric binary matrix, and the clusters for the rows of
$\boldsymbol{X}$ are the same as the clusters for the columns of
$\boldsymbol{X}$. Bickel and Chen (2009) used methods similar to those used
for proving consistency of M-estimators to derive results for network
clustering when $n$ tends to infinity. This work was later extended by Choi,
Wolfe and Airoldi (2012), who allow the number of clusters to increase with
$n$; Zhao, Levina and Zhu (2011), who allow for nodes not belonging to any
cluster; and Zhao, Levina and Zhu (2012), who incorporate individual-specific
effects. In parallel to this work, Rohe, Chatterjee and Yu (2011) study the
performance of spectral clustering for symmetric binary networks; and Rohe and
Yu (2012) study spectral methods for unsymmetric binary networks.

The main contribution of our work is recognizing that methods originally
developed for an extremely specialized context (symmetric binary networks) can
be generalized to handle clustering for arbitrary data matrices. Using a
Lindeberg-type condition, we are able to generalize the Bickel and Chen (2009)
results beyond Bernoulli random variables. We have also managed to simplify
the presentation of the main proof. To our knowledge, this is the first time
methodologies for binary networks have been used to study general biclustering
methods.

The main text of the paper is organized as follows. First Section 2 describes
the theoretical setup and Section 3 presents the consistency theorem. Section
4 presents the main steps of the proof. Next, Section 5 corroborates the
theoretical findings through a simulation study, and Section 6 presents an
application to a microarray dataset. Finally, Section 7 presents some
concluding remarks. The appendices include additional technical and empirical
results.

##  2 Block models and profile likelihood

As above, let $\boldsymbol{X}=[X_{ij}]\in\mathbb{R}^{m\times n}$ be a data
matrix. One way to formalize the biclustering problem is to follow the network
clustering literature and posit existence of $K$ _row classes_ and $L$ _column
classes_ , such that the mean value of entry $X_{ij}$ is determined solely by
the classes of row $i$ and column $j$. That is, there is an unknown row class
membership vector $\boldsymbol{c}\in K^{m}$, an unknown column class
membership vector $\boldsymbol{d}\in L^{n}$, and an unknown mean matrix
$\boldsymbol{M}=[\mu_{kl}]\in\mathbb{R}^{K\times L}$ such that

$\operatorname{E}X_{ij}=\mu_{c_{i}d_{j}};$

we refer to this model as a _block model_ , after the related model for
undirected networks proposed by Holland, Laskey and Leinhardt (1983). With a
block model, our goal is to estimate $\boldsymbol{c}$ and $\boldsymbol{d}$. We
do this by assigning labels to the rows and columns of $\boldsymbol{X}$,
codified in vectors $\boldsymbol{g}\in K^{m}$ and $\boldsymbol{h}\in L^{n}$.
Ideally, $\boldsymbol{g}$ and $\boldsymbol{h}$ match $\boldsymbol{c}$ and
$\boldsymbol{d}$. Note that we are assuming that the true numbers of row and
column clusters, $K$ and $L$, are fixed and known.

Analogously to Bickel and Chen (2009), we can employ profile-likelihood
(Murphy and van der Vaart, 2000) for the biclustering task. However, we cannot
assume that the elements of $\boldsymbol{X}$ are Bernoulli random variables.
Instead, we proceed by initially proposing a simple distribution for
$\boldsymbol{X}$ and we deriving the maximum profile likelihood estimator for
$\boldsymbol{c}$ and $\boldsymbol{d}$ under this model. Later, we consider a
far more general distribution for $\boldsymbol{X}$, and we show that the
simple profile-likelihood based estimator is still consistent, even under
model misspecification.

Suppose that the elements of $\boldsymbol{X}$ are independent draws whose
distribution is specified by a single-paremeter exponential family.
Conditional on $\boldsymbol{c}$ and $\boldsymbol{d}$, entry $X_{ij}$ has
density $g(x;\eta_{c_{i}d_{j}})$ with respect to some $\sigma$-finite measure
$\nu$, where

$g(x;\eta)=\exp\\{x\eta-\psi(\eta)\\};$

$\psi(\eta)$ is the cumulant function, and
$\eta_{kl}=(\psi^{\prime})^{-1}(\mu_{kl})$. With labels $\boldsymbol{g}$ and
$\boldsymbol{h}$, the complete data log-likelihood is

$\operatorname{l}(\boldsymbol{g},\boldsymbol{h},\boldsymbol{M})=mn\sum_{k,l}\hat{p}_{k}\,\hat{q}_{l}\,\\{\bar{X}_{kl}\,\eta_{kl}-\psi(\eta_{kl})\\},$

where

$\hat{p}_{k}=\frac{1}{m}\sum_{i}\operatorname{I}(g_{i}=k),\qquad\hat{q}_{l}=\frac{1}{n}\sum_{j}\operatorname{I}(h_{j}=l),$

and

$\bar{X}_{kl}=\frac{\sum_{i,j}X_{ij}\operatorname{I}(g_{i}=k,h_{j}=l)}{\sum_{i,j}\operatorname{I}(g_{i}=k,h_{j}=l)}.$

The profile log-likelihood is

$\operatorname{pl}(\boldsymbol{g},\boldsymbol{h})=\sup_{\boldsymbol{M}}\operatorname{l}(\boldsymbol{g},\boldsymbol{h},\boldsymbol{M})=mn\sum_{k,l}\hat{p}_{k}\,\hat{q}_{l}\,\psi^{\ast}(\bar{X}_{kl}),$

where $\psi^{\ast}(x)=\sup_{\eta}\\{x\eta-\psi(\eta)\\}$ is the convex
conjugate of $\psi$, known as the _rate function_ in the large deviations
literature (Dembo and Zeitouni, 1998). Following the above derivation, it is
natural to estimate $\boldsymbol{c}$ and $\boldsymbol{d}$ by the label vectors
which maximize $\operatorname{pl}(\boldsymbol{g},\boldsymbol{h})$.

In the sequel, we consider a far more general setting. We consider criterion
functions of the form

$F(\boldsymbol{g},\boldsymbol{h})=mn\sum_{kl}\hat{p}_{k}\,\hat{q}_{l}\,f(\bar{X}_{kl}/\rho),$
(2.1)

where $f$ is any smooth convex function and $\rho$ is a scale parameter.
Borrowing terminology from the large deviations literature, we refer to $f$ as
the rate function of the criterion; though, since we allow $f$ to take
negative values, we do not require that $f$ be a rate function in the
strictest sense. We permit the elements of $\boldsymbol{X}$ to have different
distributions, allowing for heteroscedasticity and model misspecification. We
show that under mild technical conditions, the maximizer of $F$ is a
consistent estimator of the true row and column classes.

##  3 Consistency results

To prove consistency, we need to work with a sequence $\boldsymbol{X}_{n}$ of
data matrices. Suppose that $\boldsymbol{X}_{n}\in\mathbb{R}^{m\times n}$ and
$m=m(n)$ with $n/m\to\gamma$ for some finite constant $\gamma>0$. Suppose that
for each $n$ there exists a row class membership vector $\boldsymbol{c}_{n}\in
K^{m}$ and a column class membership vector $\boldsymbol{d}_{n}\in L^{n}$. We
assume that there exist vectors $\boldsymbol{p}\in\mathbb{R}^{K}$ and
$\boldsymbol{q}\in\mathbb{R}^{L}$ such that $\hat{p}_{k}(\boldsymbol{c})\to
p_{k}$ and $\hat{q}_{l}(\boldsymbol{d})\to q_{l}$ as $n\to\infty$ almost
surely for all k and l. This assumption is satisfied, for example, if the
class labels are independently drawn from a multinomial distribution. When
there is no ambiguity, we omit the subscript $n$.

Assume that the mean of element $X_{ij}$ depends only on the row and column
memberships $c_{i}$ and $d_{j}$, so that

$\operatorname{E}(X_{ij}\mid\boldsymbol{c},\boldsymbol{d})=\mu_{c_{i}d_{j}}$

for some matrix $\boldsymbol{M}=[\mu_{kl}]\in\mathbb{R}^{K\times L}$, possibly
varying with $n$. To model sparsity in $\boldsymbol{X}$, we allow
$\boldsymbol{M}$ to tend to $\mathbf{0}$. To avoid degeneracy, we suppose that
there exists a sequence $\rho$ and a fixed matrix
$\boldsymbol{M}_{0}\in\mathbb{R}^{K\times L}$ such that
$\rho^{-1}\boldsymbol{M}\to\boldsymbol{M}_{0}$.

Define the normalized confusion matrices
$\boldsymbol{C}(\boldsymbol{g})\in\mathbb{R}^{K\times K}$ and
$\boldsymbol{D}(\boldsymbol{h})\in\mathbb{R}^{L\times L}$ by

$C_{ak}(\boldsymbol{g})=\frac{1}{m}\sum_{i}I(c_{i}=a,g_{i}=k),\qquad
D_{bl}(\boldsymbol{h})=\frac{1}{n}\sum_{j}I(d_{j}=b,h_{j}=l).$

Entry $C_{ak}$ is the proportion of nodes with class $a$ and label $k$; entry
$D_{bl}$ is defined similarly. These matrices are normalized so that
$\boldsymbol{C}^{T}\boldsymbol{1}=\hat{\boldsymbol{p}}$ and
$\boldsymbol{D}^{T}\boldsymbol{1}=\hat{\boldsymbol{q}}$.

We only consider nontrivial partitions; to this end, for $\varepsilon>0$,
define

$\mathcal{J}_{\varepsilon}=\\{\boldsymbol{g},\boldsymbol{h}:\hat{p}_{k}(\boldsymbol{g})>\varepsilon,\hat{q}_{l}(\boldsymbol{h})>\varepsilon\\}.$

For fixed convex rate function $f$, we let $F$ be a criterion function as in
(2.1).

###  3.1 Assumptions

Denote by $\mathcal{M}_{0}\in\mathbb{R}$ the convex hull of the entries of
$\boldsymbol{M}_{0}$. Let $\mathcal{M}$ be a neighborhood of
$\mathcal{M}_{0}$. We require the following assumptions:

  1. (A1)

The biclusters are identifiable. No two rows of $\boldsymbol{M}_{0}$ are
equal, and no two columns of $\boldsymbol{M}_{0}$ are equal.

  2. (A2)

The rate function is locally strictly convex. That is,
$f^{\prime\prime}(\mu)>0$ when $\mu\in\mathcal{M}$.

  3. (A3)

The third derivative of the rate function is locally bounded. That is, $\lvert
f^{\prime\prime\prime}(\mu)\rvert$ is bounded when $\mu\in\mathcal{M}$.

  4. (A4)

The average variance of the elements is of the same order as $\rho$.

$\limsup_{n\to\infty}\frac{1}{\rho
mn}\sum_{i,j}\operatorname{E}[(X_{ij}-\mu_{c_{i}d_{j}})^{2}\mid\boldsymbol{c},\boldsymbol{d}]<\infty.$

  5. (A5)

The mean matrix does not converge to zero too quickly. That is,
$\limsup_{n\to\infty}\rho n=\infty$.

  6. (A6)

The elements satisfy a Lindeberg condition. For all $\varepsilon>0$,

$\lim_{n\to\infty}\frac{1}{\rho^{2}mn}\sum_{i,j}\operatorname{E}[(X_{ij}-\mu_{c_{i}d_{j}})^{2}\operatorname{I}(\lvert
X_{ij}-\mu_{c_{i}d_{j}}\rvert>\varepsilon\sqrt{mn}\rho)\mid\boldsymbol{c},\boldsymbol{d}]=0.$

Assumption (A4) is trivially satisfied for dense data and is satisfied for
Binomial and Poisson data so long as
$\rho^{-1}\boldsymbol{M}\to\boldsymbol{M}_{0}$. However, this assumption
cannot handle arbitrary sparsity. For example, if the elements of
$\boldsymbol{X}$ are distributed as Negative Binomial random variables, then
assumption (A4) requires that the product of the mean and the dispersion
parameter does not tend to infinity. In other words, for sparse count data,
the counts cannot be too heterogeneous.

Assumption (A5) places a sparsity restriction on the mean matrix. Zhao, Levina
and Zhu (2012) used the same assumption to establish weak consistency for
network clustering. A variant Lyaponuv’s condition (Varadhan, 2001) implies
(A6). That is, if

$\lim_{n\to\infty}\frac{1}{(\rho\sqrt{mn})^{2+\delta}}\sum_{i,j}\operatorname{E}|X_{ij}-\mu_{c_{i}d_{j}}|^{2+\delta}=0$

for some $\delta>0$, then (A6) follows. In particular, for dense data ($\rho$
bounded away from zero), uniformly bounded $(2+\delta)$ absolute central
moments for any $\delta>0$ is sufficient. For certain types of sparse data
(Bernoulli or Poisson data with $\rho$ converging to zero), (A5) is a
sufficient condition for (A6).

###  3.2 Main results

Given the above setup and assumptions, it follows that the maximizer of
$F(\boldsymbol{g},\boldsymbol{h})$ is a consistent estimator of the true row
and column labels.

######  Theorem 3.1.

Fix any $\varepsilon>0$ with $\varepsilon<\min_{a}\\{p_{a}\\}$ and
$\varepsilon<\min_{b}\\{q_{b}\\}$. Let
$(\boldsymbol{\hat{g}},\boldsymbol{\hat{h}})$ satisfy
$F(\boldsymbol{\hat{g}},\boldsymbol{\hat{h}})=\max_{\mathcal{J}_{\varepsilon}}F(\boldsymbol{g},\boldsymbol{h})$.
If assumptions (A1)–(A6) hold, then all limit points of
$\boldsymbol{C}(\boldsymbol{\hat{g}})$ and
$\boldsymbol{D}(\boldsymbol{\hat{h}})$ are permutations of diagonal matrices,
i.e. the proportions of mislabeled rows and columns converge to zero in
probability.

The statement of the theorem is somewhat awkward, because we can permute the
row or column labels and get the same value of the criterion function, $F$.
Thus, $F$ has multiple maximizers, and it is only possible to recover the true
classes up to a permutation of labels.

From the discussion in Section 2, it follows that maximizing the profile log-
likelihood associated with a single-parameter exponential family may be a
reasonable biclustering procedure. Furthermore, the proof of this result does
not require the distribution of the data to be correctly specified and allows
for possible model misspecification so long as (A1)-(A6) are satisfied. This
implies that this result can be applied to binary matrices, count data, and
continuous data which could be reasonably modeled by Binomial, Poisson, and
Gaussian distributed data, respectively, making this a suitable result for our
motivating examples.

Under stronger distributional assumptions, we can use the methods of the proof
to establish finite-sample results. For example, if we assume that the
elements of $\boldsymbol{X}$ are Gaussian, then the following result holds.
The proof of the finite-sample result follows the same outline as the
asymptotic result; see Appendix B.

######  Theorem 3.2.

Fix any $\varepsilon>0$. Let $(\hat{\boldsymbol{g}},\hat{\boldsymbol{h}})$
satisfy
$F(\hat{\boldsymbol{g}},\hat{\boldsymbol{h}})=\max_{\mathcal{J}_{\varepsilon}}F(\boldsymbol{g},\boldsymbol{h})$.
If the elements of $\boldsymbol{X}$ are independent Gaussian random variables
with constant variance $\sigma^{2}$ and assumptions (A1)-(A3) hold, then for
any
$0<\delta<\min\Big{\\{}1,\frac{8c\sigma\min\\{K^{2},L^{2}\\}}{\tau\varepsilon^{2}}\Big{\\}}$,

$\Pr\Big{(}\big{(}\boldsymbol{C}(\boldsymbol{\hat{g}}),\boldsymbol{D}(\boldsymbol{\hat{h}})\big{)}\notin\mathcal{P}_{\delta}\cap\mathcal{Q}_{\delta}\Big{)}\leq
2K^{m+1}L^{n+1}\exp\Big{\\{}-\frac{T_{n}\tau^{2}\varepsilon^{4}\delta^{2}}{256c^{2}\sigma^{2}\min\\{K^{4},L^{4}\\}}\Big{\\}},$

where $c=\sup|f^{\prime}(\mu)|$ for $\mu$ in $\mathcal{M}$,

$T_{n}=\inf_{k,l}\Big{\\{}\sum_{i}\sum_{j}\operatorname{I}(g_{i}=k,h_{j}=l)|(\boldsymbol{g},\boldsymbol{h})\in\mathcal{J}_{\varepsilon}\Big{\\}}.$

and $\mathcal{P}_{\delta}$ and $\mathcal{Q}_{\delta}$ are neighborhoods of the
optimal confusion matrices defined in the sequel in equation (4.1).

To prove the main theorem, we first establish that in the limit, $F$ is close
to a nonrandom “population version,” $G$. Then, we establish that $G$ is
maximized at the true class labels. Finally, we show that outside any
neighborhood around the true class labels, $G$ is smaller than at the true
values.

As alluded to in the Introduction, the main thrust of the proof is similar to
that used in the literature on clustering for symmetric binary networks
initiated by Bickel and Chen (2009) and extended by Choi, Wolfe and Airoldi
(2012), Zhao, Levina and Zhu (2011) and Zhao, Levina and Zhu (2012). There are
three main points of departure from this previous work. First, the literature
on clustering for symmetric binary networks works only with Bernoulli random
variables, whereas we consider arbitrary distributions that satisfy a
Lindeberg-type condition. Second, these papers work with uniformly Lipschitz
criterion functions, while we consider potentially unbounded criterion
functions. Lastly, the results for symmetric binary networks require an
additional assumption: that the criterion function is maximized at the
population parameters. We show that this condition is satisfied as long as the
parameters are identifiable and the rate function is strictly convex.

##  4 Proof of consistency theorem

In this section, we continue with the setup and notation of Section 3. To
prove Theorem 3.1, we need to introduce some additional notation.

Define expectation matrix
$\boldsymbol{E}(\boldsymbol{g},\boldsymbol{h})\in\mathbb{R}^{K\times L}$ with

$E_{kl}(\boldsymbol{g},\boldsymbol{h})=\operatorname{E}(\bar{X}_{kl}(\boldsymbol{g},\boldsymbol{h})\mid\boldsymbol{c},\boldsymbol{d})=\frac{[\boldsymbol{C}^{T}\,\boldsymbol{M}\,\boldsymbol{D}]_{kl}}{[\boldsymbol{C}^{T}\boldsymbol{1}]_{k}[\boldsymbol{D}^{T}\boldsymbol{1}]_{l}},$

where $\boldsymbol{C}=\boldsymbol{C}(\boldsymbol{g})$ and
$\boldsymbol{D}=\boldsymbol{D}(\boldsymbol{h})$. Also, define normalized
residual matrix
$\boldsymbol{R}(\boldsymbol{g},\boldsymbol{h})\in\mathbb{R}^{K\times L}$ by

$\boldsymbol{R}(\boldsymbol{g},\boldsymbol{h})=\rho^{-1}\\{\bar{\boldsymbol{X}}(\boldsymbol{g},\boldsymbol{h})-\boldsymbol{E}(\boldsymbol{g},\boldsymbol{h})\\}.$

The weak law of large numbers establishes that for fixed $\boldsymbol{g}$ and
$\boldsymbol{h}$, the convergence
$R_{kl}(\boldsymbol{g},\boldsymbol{h})\overset{\mathit{P}}{\to}0$ holds. We
can prove a stronger result, that this convergence is uniform over all
$\boldsymbol{g}$ and $\boldsymbol{h}$.

######  Lemma 4.1.

Under assumptions (A1)–(A6), for all $\varepsilon>0$,

$\sup_{\mathcal{J}_{\varepsilon}}\|\boldsymbol{R}(\boldsymbol{g},\boldsymbol{h})\|_{\infty}\overset{\mathit{P}}{\to}0,$

where $\|\mathbf{A}\|_{\infty}=\max_{k,l}|A_{kl}|$ for any matrix
$\mathbf{A}$.

With Lemma 4.1 (proved in Appendix A), we can establish that in the limit,
$F(\boldsymbol{g},\boldsymbol{h})$ is close to its “population version,” which
depends only on $\boldsymbol{C}$ and $\boldsymbol{D}$. To define this
population version, first, for each $M$, define $\mathcal{S}_{M}$ to be the
set of $M\times M$ matrices with nonnegative entries summing to one:
$\mathcal{S}_{M}=\\{\boldsymbol{X}\in\mathbb{R}_{+}^{M\times
M}:\mathbf{1}^{T}\boldsymbol{X}\boldsymbol{1}=1\\}.$ Next, define function
$G_{\boldsymbol{M}_{0}}:\mathcal{S}_{K}\times\mathcal{S}_{L}\to\mathbb{R}$ to
be the population version of $F$:

$\displaystyle G_{\boldsymbol{M}_{0}}(\boldsymbol{C},\boldsymbol{D})$
$\displaystyle=\sum_{k,l}[\boldsymbol{C}^{T}\boldsymbol{1}]_{k}\,[\boldsymbol{D}^{T}\boldsymbol{1}]_{l}\,f\Big{(}\frac{[\boldsymbol{C}^{T}\boldsymbol{M}_{0}\boldsymbol{D}]_{kl}}{[\boldsymbol{C}^{T}\boldsymbol{1}]_{k}\,[\boldsymbol{D}^{T}\boldsymbol{1}]_{l}}\Big{)}.$

######  Lemma 4.2.

$F$ is close to its population version in the sense that, for all
$\varepsilon>0$,

$\sup_{\mathcal{J}_{\varepsilon}}|F(\boldsymbol{g},\boldsymbol{h})-G_{\boldsymbol{M}_{0}}(\boldsymbol{C},\boldsymbol{D})|\overset{\mathit{P}}{\to}0.$

This is a direct consequence of Lemma 4.1 and assumption (A3), which implies
that $f$ is locally Lipschitz; see Appendix A for details.

Once we have established that $F$ is close to its population version, our next
task is to show that the population version is maximized at the true class
labels. To do this, for $\delta>0$, first define

$\displaystyle\mathcal{P}_{\delta}$
$\displaystyle=\\{\boldsymbol{C}\in\mathcal{S}_{K}:\max_{a\neq
a^{\prime}}C_{ak}C_{a^{\prime}k}<\delta\\}$ (4.1a) and
$\displaystyle\mathcal{Q}_{\delta}$
$\displaystyle=\\{\boldsymbol{D}\in\mathcal{S}_{L}:\max_{b\neq
b^{\prime}}D_{bl}D_{b^{\prime}l}<\delta\\}.$ (4.1b)

A permutation of a diagonal matrix has only one non-zero entry in each column,
so taking $\delta$ close to zero restricts the confusion matrices to be close
to permutations of diagonal matrices.

######  Lemma 4.3.

If $\min_{a}\\{[\boldsymbol{C}\boldsymbol{1}]_{a}\\}>\eta,$
$\min_{b}\\{[\boldsymbol{D}\boldsymbol{1}]_{b}\\}>\eta,$ and
$(\boldsymbol{C},\boldsymbol{D})\notin\mathcal{P}_{\delta}\times\mathcal{Q}_{\delta}$,
then $G_{\boldsymbol{M}_{0}}(\boldsymbol{C},\boldsymbol{D})$ is small, in the
sense that

$G_{\boldsymbol{M}_{0}}(\boldsymbol{C},\boldsymbol{D})-\sum_{a,b}[\boldsymbol{C}\boldsymbol{1}]_{a}[\boldsymbol{D}\boldsymbol{1}]_{b}f([\boldsymbol{M}_{0}]_{ab})\leq-\kappa\eta^{2}\delta,$

where $\kappa$ is a constant independent of $\delta$ and $\eta$.

We prove Lemma 4.3 in Appendix A.

Next, we prove the consistency theorem.

###### Proof of Theorem 3.1.

Fix $\delta>0$ and define $\mathcal{P}_{\delta}$ and $\mathcal{Q}_{\delta}$ as
in (4.1). We will show that if
$(\boldsymbol{g},\boldsymbol{h})\in\mathcal{J}_{\varepsilon}$ and if
$(\boldsymbol{C}(\boldsymbol{g}),\boldsymbol{D}(\boldsymbol{h}))\notin(\mathcal{P}_{\delta},\mathcal{Q}_{\delta})$,
then $F(\boldsymbol{g},\boldsymbol{h})<F(\boldsymbol{c},\boldsymbol{d})$ with
probability approaching one. Moreover, this inequality holds uniformly over
all such choices of $(\boldsymbol{g},\boldsymbol{h})$. Since $\delta$ is
arbitrary, this implies that $\boldsymbol{C}(\boldsymbol{\hat{g}})$ and
$\boldsymbol{D}(\boldsymbol{\hat{h}})$ converge to permutations of diagonal
matrices, i.e. the proportions of misclassified rows and columns converge to
zero.

Set
$r_{n}=\sup_{\mathcal{J}_{\varepsilon}}|F(\boldsymbol{g},\boldsymbol{h})-G_{\boldsymbol{M}_{0}}(\boldsymbol{C}(\boldsymbol{g}),\boldsymbol{D}(\boldsymbol{h})|$.
Suppose $(\boldsymbol{g},\boldsymbol{h})\in\mathcal{J}_{\varepsilon}$. In this
case,

$\displaystyle
F(\boldsymbol{g},\boldsymbol{h})-F(\boldsymbol{c},\boldsymbol{d})$
$\displaystyle\leq
2r_{n}+\\{G_{\boldsymbol{M}_{0}}(\boldsymbol{C}(\boldsymbol{g}),\boldsymbol{D}(\boldsymbol{h}))-G_{\boldsymbol{M}_{0}}(\boldsymbol{C}(\boldsymbol{c}),\boldsymbol{D}(\boldsymbol{d}))\\}$
$\displaystyle=2r_{n}+\big{\\{}G_{\boldsymbol{M}_{0}}(\boldsymbol{C}(\boldsymbol{g}),\boldsymbol{D}(\boldsymbol{h}))-\sum_{a,b}[\boldsymbol{C}\boldsymbol{1}]_{a}[\boldsymbol{D}\boldsymbol{1}]_{b}f([\boldsymbol{M}_{0}]_{ab})\big{\\}}.$

Pick $\eta>0$ smaller than $\min_{a}\\{p_{a}\\}$ and $\min_{b}\\{q_{b}\\}$. By
assumption, the true row and column class proportions converge to
$\boldsymbol{p}$ and $\boldsymbol{q}$. Thus, for all $\boldsymbol{g}\in K^{m}$
and $\boldsymbol{h}\in L^{n}$, for $n$ large enough,
$[\boldsymbol{C}(\boldsymbol{g})]_{a}\geq\eta$ and
$[\boldsymbol{D}(\boldsymbol{h})]_{b}\geq\eta$; this holds uniformly over all
choices of $(\boldsymbol{g},\boldsymbol{h})$.

Applying Lemma 4.3, to the second term in the inequality, we get that with
probability approaching one,

$F(\boldsymbol{g},\boldsymbol{h})-F(\boldsymbol{c},\boldsymbol{d})\leq
2r_{n}-\kappa\eta^{2}\delta$

for all $(\boldsymbol{g},\boldsymbol{h})\in\mathcal{J}_{\varepsilon}$ such
that
$(\boldsymbol{C}(\boldsymbol{g}),\boldsymbol{D}(\boldsymbol{h}))\notin\mathcal{P}_{\delta}\times\mathcal{Q}_{\delta}$.
By Lemma 4.2, $r_{n}\overset{\mathit{P}}{\to}0$. Thus, with probability
approaching one,
$(\boldsymbol{C}(\boldsymbol{\hat{g}}),\boldsymbol{D}(\boldsymbol{\hat{h}}))\in\mathcal{P}_{\delta}\times\mathcal{Q}_{\delta}$.
Since this result holds for all $\delta$, all limit points of
$\boldsymbol{C}(\boldsymbol{\hat{g}})$ and
$\boldsymbol{D}(\boldsymbol{\hat{h}})$ must be permutations of diagonal
matrices. ∎

##  5 Empirical evaluation

We study the performance of profile log-likelihood for biclustering Bernoulli,
Poisson, Gaussian, and Non-standardized Student’s t data. We consider the
following three rate functions:

$\displaystyle f_{\text{Bernoulli}}(\mu)$
$\displaystyle=\mu\log\mu+(1-\mu)\log(1-\mu),$ (5.1a) $\displaystyle
f_{\text{Poisson}}(\mu)$ $\displaystyle=\mu\log\mu-\mu,$ (5.1b) $\displaystyle
f_{\text{Gaussian}}(\mu)$ $\displaystyle=\mu^{2}/2.$ (5.1c)

We use the appropriate rate function for the Bernoulli and Gaussian
simulations. To study performance under model misspecification, we compare the
performance of the Poisson and Gaussian rate functions for the Poisson
simulations and use the Gaussian rate function for the t simulations.

For the Poisson and Gaussian rate functions (5.1b) and (5.1c), the maximizer
of the criterion function (2.1) does not depend on the scale factor $\rho$.
This is immediately obvious in the Gaussian case. For the Poisson case, the
function
$f_{Poisson}(\mu/\rho)=\frac{1}{\rho}\mu\log(\mu)-\frac{1}{\rho}\mu(1+\log(\rho))$.
When summed over all biclusters, the second term in this sum is equal to a
constant so the value of $\mu$ which maximizes $f_{Poisson}(\mu/\rho)$ does
not depend on the value of $\rho$.

This is not the case for the Binomial rate function (5.1a), but the parameter
$\rho$ is not identifiable in practice so it does not make sense to try to
estimate it. For our simulations we use which maximizes $\rho=1$ in the
fitting procedure, regardless of the true scale factor for the mean matrix
$\boldsymbol{M}$. Even though in the simulations assumption (A1) doesn’t hold
for this choice of $\rho$, we still get consistency, because the maximizer of
the criterion with $f_{\text{Bernoulli}}(\mu)$ is close to the maximizer with
$f_{\text{Poisson}}(\mu)$. See Perry and Wolfe (2012) for discussion of a
related phenomenon.

To maximize the profile log-likelihood, we compute initial partitions for the
rows and columns. Then, we iteratively update the cluster assignments in a
greedy manner, based on the Kernighan-Lin heuristic (Kernighan and Lin, 1970)
employed by Newman (2006):

  1. 1.

For each row and column, we compute the optimal label assignment while keeping
the labels of all other rows and columns fixed; we also record the improvement
made by making this assignment.

  2. 2.

In order of the local improvements recorded in step 1, we perform the label
reassignments determined in step 1. Note that these assignments are no longer
locally optimal since the labels of many of the rows and columns change during
this step.

  3. 3.

Out of those labels considered in step 2, we choose the one which has the
highest profile likelihood.

We iteratively perform steps 1–3 until the profile likelihood converges. The
algorithm is not guaranteed to converge to a global optimum, but it seems to
perform well in practice. If implemented efficiently, each update to the log-
likelihood can be done in a constant number of operations. Therefore, the most
complex part of this algorithm is the ordering of the local improvements. On
average, using a conventional sorting algorithm this requires $O(n\log(n))$
operations. For comparison, the algorithmic costs for spectral-based
biclustering algorithms are at least $O(n^{2})$, the cost of computing the top
singular vectors of the data matrix using an indirect method (Golub and Loan,
1996).

In our simulations, we report the proportion of misclassified rows and columns
by the profile-likelihood-based method (PL), which Theorem 3.1 guarantees to
be consistent. For the Poisson simulations we report the profile-liklihood-
based method with the Poisson rate function (PL-Pois) and with the Gaussian
rate function (PL-Gaus). We also report misclassification errors for $k$-means
applied separately to the rows and columns (KM) and for the DI-SIM
biclustering algorithm (DS) of Rohe and Yu (2012). For the profile-likelihood-
based method, we compute initial partitions for the rows and columns by
applying $k$-means separately to the rows and the columns. In principle any
initialization can be used. The algorithm is run once for each simulation.

For the Poisson simulation, we simulate from a block model with $K=2$ row
clusters and $L=3$ column clusters. We vary the number of columns, $n$,
between 200 to 1500 and we take the number of rows as $m=\gamma n$ where
$\gamma\in\\{0.5,1,2\\}$.

We set the row and column class memberships as independent multinomials with
probabilities $\boldsymbol{p}=(0.3,0.7)$ and $\boldsymbol{q}=(0.2,0.3,0.5)$.
We choose the matrix of block parameters to be

$\boldsymbol{M}=[\mu_{ab}]=\frac{b}{\sqrt{n}}\begin{pmatrix}0.92&0.77&1.66\\\
0.17&1.41&1.45\end{pmatrix},$

where $b$ is chosen between 5 and 20; the entries of the matrix were chosen
randomly, uniformly on the interval $[0,2]$. Note that
$\boldsymbol{M}\to\boldsymbol{0}$, so the data matrix is sparse. We generate
the data conditional on the row and column classes as
$X_{ij}\mid\boldsymbol{c},\boldsymbol{d}\sim\mathrm{Poisson}(\mu_{c_{i}d_{j}}).$

Figure 1 presents the average bicluster misclassification rates for each
sample size and Tables 1 and 2 report the standard deviations. In all of the
scenarios considered, the biclustering based on the profile log-likelihood
criterion performs at least as well as the other methods and shows signs of
convergence even if the rate function is misspecified. Although the $k$-means
algorithm and the DI-SIM algorithm perform well in some settings, in other
settings the performance is not well-behaved. We also see that the standard
deviation of the misclassification rate tends to zero for the PL bicustering,
but not for the other two algorithms. The results in McSherry (2001) and Rohe
and Yu (2012) suggest that the DI-SIM algorithm may be a consistent
biclustering algorithm, but require strong assumptions that do not hold in our
simulations. Specifically, the former requires that the variance of the
entries in $\boldsymbol{M}$ to be $\gg\log^{6}(n)/n$ and the later requires
that the minimum expected degree eventually exceeds
$\sqrt{2}n/\sqrt{\log(n)}$. Therefore it is not clear if the convergence of
the DI-SIM algorithm is slow or if the results do not apply in this setting.

Appendix C describes in detail the simulations for Bernoulli, Gaussian, and t
data. For the Bernoulli simulation, the profile-likelihood-based method
outperforms the other procedures. For the Gaussian and t-distributed data, the
DI-SIM algorithm is much more competitive, but our algorithm still beats it in
almost all cases.

Overall, the simulation results corroborate the conclusions of Theorem 3.1 and
support the use of biclustering based on the profile log-likelihood criterion.

Figure 1: Average misclassification rates for Poisson example over 100
simulations. Table 1: Standard deviations for Poisson example over 100
simulations for the profile likelihood methods based on Poisson and Gaussian
criteria. $m=0.5n$
---
| PL-Pois | PL-Norm
n | b=5 | b=10 | b=20 | b=5 | b=10 | b=20
200 | 0.0683 | 0.0177 | 0.0554 | 0.1376 | 0.0254 | 0.0566
500 | 0.0112 | 0.0046 | 0.0392 | 0.0161 | 0.0069 | 0.0391
1000 | 0.0050 | 0.0013 | 0.0002 | 0.0072 | 0.0022 | 0.0004
1500 | 0.0034 | 0.0007 | 0.0001 | 0.0046 | 0.0012 | 0.0002
$m=n$
n | b=5 | b=10 | b=20 | b=5 | b=10 | b=20
200 | 0.0158 | 0.0041 | 0.0777 | 0.0481 | 0.0061 | 0.0771
500 | 0.0039 | 0.0008 | 0.0000 | 0.0064 | 0.0014 | 0.0000
1000 | 0.0015 | 0.0002 | 0.0000 | 0.0019 | 0.0004 | 0.0000
1500 | 0.0007 | 0.0000 | 0.0000 | 0.0012 | 0.0001 | 0.0000
$m=2n$
n | b=5 | b=10 | b=20 | b=5 | b=10 | b=20
200 | 0.0063 | 0.0009 | 0.1259 | 0.0010 | 0.0002 | 0.0124
500 | 0.0012 | 0.0000 | 0.0552 | 0.0002 | 0.0000 | 0.0054
1000 | 0.0002 | 0.0001 | 0.0000 | 0.0001 | 0.0000 | 0.0000
1500 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000
Table 2: Standard deviations for Poisson example over 100 simulations for the
K-means and DI-SIM algorithms. $m=0.5n$
---
| KM | DS
n | b=5 | b=10 | b=20 | b=5 | b=10 | b=20
200 | 0.1325 | 0.0635 | 0.0487 | 0.0488 | 0.0355 | 0.0385
500 | 0.0666 | 0.0094 | 0.0357 | 0.0209 | 0.0090 | 0.0384
1000 | 0.0134 | 0.0029 | 0.0005 | 0.0090 | 0.0394 | 0.1069
1500 | 0.0061 | 0.0015 | 0.0445 | 0.0065 | 0.0387 | 0.1113
$m=n$
n | b=5 | b=10 | b=20 | b=5 | b=10 | b=20
200 | 0.1205 | 0.0102 | 0.0811 | 0.0298 | 0.0521 | 0.0798
500 | 0.0124 | 0.0021 | 0.0000 | 0.0121 | 0.0534 | 0.0926
1000 | 0.0039 | 0.0401 | 0.0405 | 0.0043 | 0.1162 | 0.1303
1500 | 0.0017 | 0.0002 | 0.0561 | 0.0018 | 0.0761 | 0.1050
$m=2n$
n | b=5 | b=10 | b=20 | b=5 | b=10 | b=20
200 | 0.0104 | 0.0009 | 0.1100 | 0.0056 | 0.0007 | 0.0858
500 | 0.0017 | 0.0000 | 0.0668 | 0.0017 | 0.0002 | 0.0844
1000 | 0.0003 | 0.0000 | 0.1020 | 0.0002 | 0.0000 | 0.1174
1500 | 0.0001 | 0.0000 | 0.0934 | 0.0001 | 0.0000 | 0.1077

##  6 Application

In this section we use profile-likelihood-based biclustering to reveal
structure in two high-dimensional datasets. For each example, we maximize the
profile log-likelihood using the algorithm described in Section 5.

###  6.1 AGEMAP

Biclustering is commonly used for microarray data to visualize the activation
patterns of thousands of genes simultaneously. It is used to identify patterns
and discover distinguishing properties between genes and individuals. We use
the AGEMAP dataset (Zahn et al., 2007) to illustrate this process.

AGEMAP is a large microarray data set containing the log expression levels for
40 mice across 8932 genes measured on 16 different tissue types. For this
analysis, we restrict attention to two tissue types: cerebellum and cerebrum.
The 40 mice in the dataset belong to four age groups, with five males and five
females in each group. One of the mice is missing data for the cerebrum tissue
so it has been removed from the dataset.

To study the relationship between mice and gene expression levels, we
bicluster the 39 $\times$ 17864 residual matrix computed from the least
squares fit to the multiple linear regression model

$Y_{ij}=\beta_{0j}+\beta_{1j}A_{i}+\beta_{2j}S_{i}+\varepsilon_{ij},$

where $Y_{ij}$ is the log-activation of gene $j$ in mouse $i$, $A_{i}$ is the
age of mouse $i$, $S_{i}$ indicates if mouse $i$ is male, $\varepsilon_{ij}$
is a random error, and $(\beta_{0j},\beta_{1j},\beta_{2j})$ is a gene-specific
coefficient vector.

The entries of the residual matrix are not independent (for example, the sum
of each column is zero). Also, the responses of many genes are likely
correlated with each other. Thus, the block model required by Theorem 3.1 is
not in force, so its conclusion will not hold unless the dependence between
the residuals is negligible. In light of this caveat, the example should be
considered as exploratory data analysis.

We perform biclusterring using profile log-likelihood based on the Gaussian
criterion (5.1c). Based on the preliminary analysis in Perry and Owen (2010),
it appears that there are three mice groups. To determine an appropriate
number of gene groups we experiment with values between two and five. Using
four gene groups appears to give a reasonable representation of the data. The
heatmap presented in Figure 2 shows the results.

Figure 2: Heatmap generated from AGEMAP residual data reflecting the varying
expression patterns in the different biclusters. The colors for the matrix
entries correspond encode to the first quartile, the middle two quartiles, and
the upper quartile.

Although the expression levels for the fourth gene group appear to be fairly
neutral across the three mouse groups, the first three gene groups each have a
visually apparent pattern. It appears that a mouse can have high expression
levels at most two of the first three gene groups. Mouse group 1 has high
expression for gene groups 2 and 3; mouse group 2 has high expression for gene
group 1; and mouse group 3 has high expression for gene groups 1 and 3.

The three clusters of mice agree with those found by Perry and Owen (2010).
That analysis identified the mouse clusters, but could not attribute meaning
to them. The bicluster-based analysis has deepend our understanding of the
three mouse clusters while suggesting some interesting interactions between
the genes.

###  6.2 MovieLens

Collaborative filtering uses the behavior of similar consumers to provide
better recommendations of products to new individuals. Since consumer habits
likely vary depending on products, biclustering can help identify structure in
the data and identify groups of consumers and groups of products with similar
patterns. As an application of this we apply biclustering to the MovieLens
data set (GroupLens Research Project, 1998). Appendix D describes the results
from this application which splits individuals who rate movies from all time
periods from individuals who only rate recent movies.

##  7 Discussion

In this report we developed a statistical setting for studying the performance
of biclustering algorithms. Under the assumption that the data follows a
stochastic block model, we derived sufficient conditions for an algorithm
based on maximizing the profile-likelihood to be consistent. This is the first
theoretical guarantee for any biclustering algorithm which can be applied to a
broad range of data distributions and can handle both sparse and dense data
matrices. These results can easily be extended to symmetric matrices. Our
empirical comparisons demonstrated that the method performs well in a variety
of situations and can outperform existing procedures.

It is important to note that the theoretical results operate under the
assumption that the true number of row and column classes are known. In
practice, this is a simplifying assumption and the optimal number of
biclusters needs to be inferred from the data. While this is an open problem
in the biclustering literature, the formalization of biclustering as an
estimation process may help identify parallels that can be drawn from
classical model selection procedures to this setting and is an interesting
topic for future research.

Applying the profile-likelihood based biclustering algorithm to real data
revealed several interesting findings. Biclustering the genes and mice in the
AGEMAP data exposed an interesting pattern in the expression of certain genes
and we found that at most two gene groups can have high expression levels for
any one mouse. Our results from the MovieLens dataset identified a distinct
difference in movie review behavior and split individuals who only rate recent
movies from individuals who rate movies from all time periods. The consistency
theorem proved in this report gives conditions under which we can have
confidence in the robustness of these findings.

##  Appendix A Additional technical results

######  Lemma A.1.

For each $n$, let $X_{n,m}$, $1\leq m\leq n$, be independent random variables
with $\operatorname{E}X_{n,m}=0$. Let $\rho_{n}$ be a sequence of positive
numbers. Let $\mathcal{I}_{n}$ be a subset of the powerset $2^{[n]}$, with
$\inf\\{|I|:I\in\mathcal{I}_{n}\\}\geq L_{n}.$ Suppose

  1. (i)

$\frac{1}{n\rho_{n}}\sum_{m=1}^{n}\operatorname{E}\lvert X_{n,m}\rvert^{2}$ is
uniformly bounded in $n$;

  2. (ii)

For all $\varepsilon>0$,
$\frac{1}{n\rho_{n}^{2}}\sum_{m=1}^{n}\operatorname{E}(\lvert
X_{n,m}\rvert^{2};\lvert X_{n,m}\rvert>\varepsilon\sqrt{n}\rho_{n})\to 0;$

  3. (iii)

$\varlimsup_{n\to\infty}\frac{n}{L_{n}}<\infty;$

  4. (iv)

$\varlimsup_{n\to\infty}\frac{\log\lvert\mathcal{I}_{n}\rvert}{\sqrt{n}}<\infty.$

  5. (v)

$\varlimsup_{n\to\infty}\rho_{n}\sqrt{n}=\infty.$

Then

$\sup_{I\in\mathcal{I}_{n}}\Big{|}\frac{1}{\rho_{n}\lvert I\rvert}\sum_{m\in
I}X_{n,m}\Big{|}\overset{\mathit{P}}{\to}0.$

###### Proof.

Let $\varepsilon>0$ be arbitrary. Define
$Y_{n,m}=\rho_{n}^{-1}X_{n,m}\operatorname{I}(\lvert
X_{n,m}\rvert\leq\varepsilon\sqrt{n}\rho_{n}),$ and note that

$\Pr(Y_{n,m}\neq\rho_{n}^{-1}X_{n,m}\text{ for some $1\leq m\leq
n$})\leq\sum_{m=1}^{n}\Pr(\lvert X_{n,m}\rvert>\varepsilon\sqrt{n}\rho_{n})\\\
\leq\frac{1}{\varepsilon^{2}n\rho_{n}^{2}}\sum_{m=1}^{n}\operatorname{E}(\lvert
X_{n,m}\rvert^{2};\lvert X_{n,m}\rvert>\varepsilon\sqrt{n}\rho_{n}).$

Fix an arbitrary $t>0$. Set $\mu_{n,m}=\operatorname{E}Y_{n,m}$ and for
$I\in\mathcal{I}_{n}$ define

$\mu_{n}(I)=\frac{1}{\lvert I\rvert}\sum_{m\in I}\mu_{n,m}.$

For $I\in\mathcal{I}_{n}$, write

$\Pr\Big{(}\sum_{m\in I}Y_{n,m}>t\,\lvert I\rvert\Big{)}=\Pr\Big{(}\sum_{m\in
I}(Y_{n,m}-\mu_{n,m})>\lvert I\rvert\big{(}t-\mu_{n}(I)\big{)}\Big{)}.$

Note that since $\operatorname{E}X_{n,m}=0$, it follows that

$\lvert\mu_{n,m}\rvert=\lvert-\operatorname{E}(\rho_{n}^{-1}X_{n,m};\lvert
X_{n,m}\rvert>\varepsilon\sqrt{n}\rho_{n})\rvert\leq\frac{1}{\varepsilon\sqrt{n}\rho^{2}_{n}}\operatorname{E}(\lvert
X_{n,m}\rvert^{2};\lvert X_{n,m}\rvert>\varepsilon\sqrt{n}\rho_{n}).$

Thus, by (ii) and (iii) we have that
$\sup_{I\in\mathcal{I}_{n}}\\{\lvert\mu_{n}(I)\rvert\\}\to 0;$ in particular,
for $n$ large enough,
$\sup_{I\in\mathcal{I}_{n}}\\{\lvert\mu_{n}(I)\rvert\\}<\frac{t}{2}.$
Consequently, for $n$ large enough, uniformly for all $I$,

$\Pr\Big{(}\sum_{m\in I}Y_{n,m}>t\,\lvert
I\rvert\Big{)}\leq\Pr\Big{(}\sum_{m\in I}(Y_{n,m}-\mu_{n,m})>t\,\lvert
I\rvert/2\Big{)}.$

Similarly,

$\Pr\Big{(}\sum_{m\in I}Y_{n,m}<-t\,\lvert
I\rvert\Big{)}\leq\Pr\Big{(}\sum_{m\in I}(Y_{n,m}-\mu_{n,m})<-t\,\lvert
I\rvert/2\Big{)}.$

We apply Bernstein’s inequality to the bound. Define
$\sigma_{n,m}^{2}=\operatorname{E}(Y_{n,m}-\mu_{n,m})^{2}$ and
$v_{n}(I)=\sum_{m\in I}\sigma_{n,m}^{2}$. Note that $\lvert
Y_{n,m}-\mu_{n,m}\rvert\leq 2\varepsilon\sqrt{n}$. By Bernstein’s inequality,

$\Pr\Big{(}\Big{|}\sum_{m\in I}(Y_{n,m}-\mu_{n,m})\Big{|}>t\,\lvert
I\rvert/2\Big{)}\leq 2\exp\Big{\\{}-\frac{t^{2}\lvert
I\rvert^{2}/8}{v_{n}(I)+\varepsilon t\lvert I\rvert\sqrt{n}/3}\Big{\\}}.$

By (i), (iv), and (v), it follows that for $n$ large enough,
$v_{n}(I)<\varepsilon t\lvert I\rvert\sqrt{n}/3,$ so

$\Pr\Big{(}\Big{|}\sum_{m\in I_{n}}(Y_{n,m}-\mu_{n,m})\Big{|}>t\lvert
I\rvert/2\Big{)}\leq 2\exp\Big{\\{}-\frac{t\lvert
I\rvert}{8\varepsilon\sqrt{n}}\Big{\\}}.$

We use this bound for each $I$ to get the union bound:

$\Pr\Big{(}\sup_{I\in\mathcal{I}_{n}}\Big{|}\frac{1}{\lvert I\rvert}\sum_{m\in
I}Y_{n,m}\Big{|}>t\Big{)}\leq
2\lvert\mathcal{I}_{n}\rvert\exp\Big{\\{}-\frac{tL_{n}}{8\varepsilon\sqrt{n}}\Big{\\}}=2\exp\Big{\\{}\log\lvert\mathcal{I}_{n}\rvert-\frac{tL_{n}}{8\varepsilon\sqrt{n}}\Big{\\}}.$

By (iii) and (iv), it is possible to choose $\varepsilon$ such that the right
hand side goes to zero. It follows then that

$\Pr\Big{(}\sup_{I\in\mathcal{I}_{n}}\Big{|}\frac{1}{\rho_{n}\lvert
I\rvert}\sum_{m\in
I}X_{n,m}\Big{|}>t\Big{)}\leq\Pr(Y_{n,m}\neq\rho_{n}^{-1}X_{n,m}\text{ for
some $1\leq m\leq n$})\\\
+\Pr\Big{(}\sup_{I\in\mathcal{I}_{n}}\Big{|}\frac{1}{\lvert I\rvert}\sum_{m\in
I}Y_{n,m}\Big{|}>t\Big{)}\to 0.$

∎

###### Proof of Lemma 4.1.

For all $t>0$,

$\Pr\Big{(}\sup_{\mathcal{J}_{\varepsilon}}\lVert\boldsymbol{R}(\boldsymbol{g},\boldsymbol{h})\rVert_{\infty}>t\Big{)}\\\
\leq KL\Pr\Big{(}\sup_{I\in\mathcal{I}_{n}}\rho^{-1}\Big{|}\sum_{\\{i,j\\}\in
I}(X_{ij}-\mu_{c_{i}d_{j}})\Big{|}>t\lvert I\rvert\Big{)},$

where $\mathcal{I}_{n}\subset 2^{[n]}\times 2^{[m]}$ is the set of all
biclusters such that $\hat{p}_{k}>\varepsilon$ for all $k$ and
$\hat{q}_{l}>\varepsilon$ for all $l$. Since $\mathcal{I}_{n}$ is a subset of
the powerset $2^{[nm]}$, by Lemma A.1, it follows that

$\Pr\Big{(}\sup_{\mathcal{J}_{\varepsilon}}\lVert\boldsymbol{R}(\boldsymbol{g},\boldsymbol{h})\rVert_{\infty}>t\Big{)}\rightarrow
0.$

∎

###### Proof of Lemma 4.2.

The technical assumptions of $f$ imply that its first derivative is bounded.
Therefore, $f$ is locally Lipschitz continuous with Lipschitz constant
$c=\sup|f^{\prime}(\mu)|$ for $\mu$ in a neighborhood of $\mathcal{M}$ and

$\left|F(\rho^{-1}\bar{\boldsymbol{X}}(\boldsymbol{g},\boldsymbol{h}),\hat{\boldsymbol{p}}(\boldsymbol{g}),\hat{\boldsymbol{q}}(\boldsymbol{h}))-G_{\boldsymbol{M}_{0}}(\boldsymbol{C},\boldsymbol{D})\right|\leq
c\|\boldsymbol{R}(\boldsymbol{g},\boldsymbol{h})\|_{\infty}.$

From Lemma 4.1, the righthand side converges to zero almost surely and the
result follows. ∎

######  Lemma A.2.

Let $\boldsymbol{C}=[C_{ak}]\in\mathbb{R}^{A\times K}$ and
$\boldsymbol{D}=[D_{bl}]\in\mathbb{R}^{B\times L}$ be nonnegative matrices
whose entries sum to one. Let $\boldsymbol{M}=[\mu_{ab}]\in\mathbb{R}^{A\times
B}$ be a matrix with no two identical columns. Define $\mathcal{M}$ to be the
convex hull of the entries of $\boldsymbol{M}$ and let
$f:\mathbb{R}\to\mathbb{R}$ be a twice differentiable convex function with
$f^{\prime\prime}$ bounded away from zero in $\mathcal{M}$. Suppose that
$[\boldsymbol{C}\boldsymbol{1}]_{a}\geq\eta$ for all $a$ and that
$D_{bl}D_{b^{\prime}l}>\delta$ for some $l$ and $b\neq b^{\prime}$. There
exists a positive constant $C$ depending on $\boldsymbol{M}$ and $f$ such that

$\sum_{k=1}^{K}\sum_{l=1}^{L}[\boldsymbol{C}^{T}\boldsymbol{1}]_{k}[\boldsymbol{D}^{T}\boldsymbol{1}]_{l}f\bigg{(}\frac{[\boldsymbol{C}^{T}\boldsymbol{M}\boldsymbol{D}]_{kl}}{[\boldsymbol{C}^{T}\boldsymbol{1}]_{k}[\boldsymbol{D}^{T}\boldsymbol{1}]_{l}}\bigg{)}\\\
\leq\sum_{a=1}^{A}\sum_{b=1}^{B}[\boldsymbol{C}\boldsymbol{1}]_{a}[\boldsymbol{D}\boldsymbol{1}]_{b}f(\mu_{ab})-\frac{C\eta^{2}\delta}{K^{2}}.$

###### Proof.

Let $l$, and $b\neq b^{\prime}$ be such that $D_{bl}D_{b^{\prime}l}>\delta$.
Since no two columns of $\boldsymbol{M}$ are identical, there exists an $a$
such that $\mu_{ab}\neq\mu_{ab^{\prime}}$. Let $k$ be the index of the largest
element in row $a$ of matrix $\boldsymbol{C}$; this element must be at least
as large as the mean, i.e.

$C_{ak}\geq\frac{[\boldsymbol{C}\boldsymbol{1}]_{a}}{K}\geq\frac{\eta}{K}.$

Let
$W=[\boldsymbol{C}^{T}\boldsymbol{1}]_{k}[\boldsymbol{D}^{T}\boldsymbol{1}]_{l}$;
this is nonzero. Now, there exists $\mu_{\ast}\in\mathcal{M}$ such that

$[\boldsymbol{C}^{T}\boldsymbol{M}\boldsymbol{D}]_{kl}=C_{ak}D_{bl}\mu_{ab}+C_{ak}D_{b^{\prime}l}\mu_{ab^{\prime}}+(W-C_{ak}D_{bl}-C_{ak}D_{b^{\prime}l})\mu_{\ast}.$

Let $z=\frac{[\boldsymbol{C}^{T}\boldsymbol{M}\boldsymbol{D}]_{kl}}{W}$. Set
$\kappa=\inf_{\mu\in\mathcal{M}}f^{\prime\prime}(\mu)$ and define
$\boldsymbol{N}=[\nu_{ab}]\in\mathbb{R}^{A\times B}$ with
$\nu_{ab}=f(\mu_{ab})$. By a second-order Taylor expansion,

$\displaystyle\frac{[\boldsymbol{C}^{T}\boldsymbol{N}\boldsymbol{D}]_{kl}}{W}-f(z)$
$\displaystyle\geq\frac{1}{2}\kappa\Big{(}\frac{C_{ak}D_{bl}}{W}(\mu_{ab}-z)^{2}+\frac{C_{ak}D_{b^{\prime}l}}{W}(\mu_{ab^{\prime}}-z)^{2}\Big{)}$
$\displaystyle\geq\kappa\frac{C^{2}_{ak}D_{bl}D_{b^{\prime}l}}{W^{2}}\Big{(}\frac{1}{2}(\mu_{ab}-z)^{2}+\frac{1}{2}(z-\mu_{ab^{\prime}})^{2}\Big{)}$
$\displaystyle\geq\kappa\frac{C^{2}_{ak}D_{bl}D_{b^{\prime}l}}{W^{2}}\Big{(}\frac{1}{2}(\mu_{ab}-z)+\frac{1}{2}(z-\mu_{ab^{\prime}})\Big{)}^{2}$
$\displaystyle=\kappa\frac{C^{2}_{ak}D_{bl}D_{b^{\prime}l}}{4W^{2}}(\mu_{ab}-\mu_{ab^{\prime}})^{2}.$

It follows that

$\displaystyle[\boldsymbol{C}^{T}\boldsymbol{1}]_{k}[\boldsymbol{D}^{T}\boldsymbol{1}]_{l}f\bigg{(}\frac{[\boldsymbol{C}^{T}\boldsymbol{M}\boldsymbol{D}]_{kl}}{[\boldsymbol{C}^{T}\boldsymbol{1}]_{k}[\boldsymbol{D}^{T}\boldsymbol{1}]_{l}}\bigg{)}-[\boldsymbol{C}^{T}\boldsymbol{N}\boldsymbol{D}]_{kl}$
$\displaystyle\leq-\frac{\kappa}{2}(\mu_{ab}-\mu_{ab^{\prime}})^{2}\frac{C_{ak}^{2}D_{bl}D_{b^{\prime}l}}{W}$
$\displaystyle\leq-\frac{\kappa\eta^{2}\delta}{4K^{2}}(\mu_{ab}-\mu_{ab^{\prime}})^{2}.$

This inequality only holds for one particular choice of $k$ and $l$; for other
choices, the left hand side is nonpositive by Jensen’s inequality. The result
of the theorem follows, with $C$ defined by

$C=\frac{\kappa}{4}\min_{a,b\neq b^{\prime}}(\mu_{ab}-\mu_{ab^{\prime}})^{2}.$

∎

###### Proof of Lemma 4.3.

If $\boldsymbol{D}\notin\mathcal{Q}_{\delta}$, then for some $l$ and some
$b\neq b^{\prime}$, $D_{bl}D_{b^{\prime}l}\geq\delta$. By Lemma A.2,

$\sum_{k,l}[\boldsymbol{C}^{T}\boldsymbol{1}]_{k}\,[\boldsymbol{D}^{T}\boldsymbol{1}]_{l}\,f\Big{(}\frac{[\boldsymbol{C}^{T}\boldsymbol{M}_{0}\boldsymbol{D}]_{kl}}{[\boldsymbol{C}^{T}\boldsymbol{1}]_{k}\,[\boldsymbol{D}^{T}\boldsymbol{1}]_{l}}\Big{)}-\sum_{a,b}[\boldsymbol{C}\boldsymbol{1}]_{a}[\boldsymbol{D}\boldsymbol{1}]_{b}f([\boldsymbol{M}_{0}]_{a,b})\leq-\frac{\kappa_{0}\,\eta^{2}\delta}{K^{2}},$

where $\kappa_{0}$ depends only on $\boldsymbol{M}_{0}$ and $f$. Similarly, if
$\boldsymbol{C}\notin\mathcal{P}_{\delta}$, then the right hand side is
bounded by $-\kappa_{0}\,\eta^{2}\delta/L^{2}.$ The result of the lemma
follows with $\kappa=\kappa_{0}/\min\\{K^{2},L^{2}\\}$. ∎

##  Appendix B Finite Sample Results

In this appendix we derive a finite sample tail bound for the probability that
the class assignments that maximize the profile likelihood are close to the
true class labels. To proceed in this setting, we make stronger distributional
assumptions than in the asymptotic case. Specifically, we assume here that the
entries $X_{ij}|\boldsymbol{c},\boldsymbol{d}$ follow a Gaussian distribution
with mean $\mu_{c_{i}d_{j}}$ and finite variance $\sigma^{2}$. We proceed with
the notation from the main text.

######  Lemma B.1.

For all $\varepsilon>0$, if $t<\sigma$ then

$\Pr\Big{(}\sup_{\mathcal{J}_{\varepsilon}}\|\boldsymbol{R}(\boldsymbol{g},\boldsymbol{h})\|_{\infty}>t\Big{)}\leq
2K^{m+1}L^{n+1}\exp\Big{(}-\frac{L_{n}t^{2}}{4\sigma^{2}}\Big{)},$

and if $t\geq\sigma$ then

$\Pr\Big{(}\sup_{\mathcal{J}_{\varepsilon}}\|\boldsymbol{R}(\boldsymbol{g},\boldsymbol{h})\|_{\infty}>t\Big{)}\leq
2K^{m+1}L^{n+1}\exp\Big{(}-\frac{L_{n}t}{4\sigma}\Big{)}$

where $\|\boldsymbol{A}\|_{\infty}=\max_{k,l}|A_{kl}|$ for any matrix
$\boldsymbol{A}$.

###### Proof.

If the entries $X_{ij}$ follow a Gaussian distribution with mean
$\mu_{c_{i}d_{j}}$ and variance $\sigma^{2}$ then

$\operatorname{E}\Big{(}|X_{ij}-\mu_{c_{i}d_{j}}|^{l}\Big{)}\leq\frac{\sigma^{2}}{2}\sigma^{l-2}l!$

so the conditions of Bernstein’s inequality hold. It follows that for any
bicluster $I$, for all $t>0$,

$\displaystyle\Pr\Big{(}\lvert\sum_{i,j\in
I}X_{ij}-\mu_{c_{i}d_{j}}\rvert>t|I|\Big{)}$ $\displaystyle\leq
2\exp\Big{\\{}-\frac{|I|^{2}t^{2}}{2(\sigma^{2}|I|+\sigma|I|t)}\Big{\\}}$
$\displaystyle\leq 2\exp\Big{\\{}-\frac{L_{n}t^{2}}{4\max\\{\sigma^{2},\sigma
t\\}}\Big{\\}}.$

Applying a union bound,

$\Pr\Big{(}\sup_{\mathcal{J}_{\varepsilon}}\|\boldsymbol{R}(\boldsymbol{g},\boldsymbol{h})\|_{\infty}>t\Big{)}\leq
2K^{m+1}L^{n+1}\exp\Big{\\{}-\frac{L_{n}t^{2}}{4\max\\{\sigma^{2},\sigma
t\\}}\Big{\\}}.$

∎

Lemma B.1 is used to establish a finite sample bound on the difference between
$F(\boldsymbol{g},\boldsymbol{h})$ and its population version.

######  Lemma B.2.

Under assumptions (A2) and (A3), for any $t>0$,

$\Pr\Big{(}\sup_{\mathcal{J}_{\varepsilon}}|F(\boldsymbol{g},\boldsymbol{h})-G_{M_{0}}(\boldsymbol{c},\boldsymbol{d})|>t\Big{)}\leq\Pr\Big{(}\sup_{\mathcal{J}_{\varepsilon}}\|\boldsymbol{R}(\boldsymbol{g},\boldsymbol{h})\|_{\infty}>\frac{t}{c}\Big{)}$

where $c=\sup|f^{\prime}(\mu)|$ for $\mu$ in $\mathcal{M}$.

Lemma B.2 is a direct consequence of the fact that $f$ is locally Lipschitz
continuous under assumptions (A2) and (A3). The details are similar to proof
of Lemma 4.2.

The next step is to show that population version is maximized at the true
class labels.

######  Lemma B.3.

Choose $\tau>0$ such that $\min_{a\neq
a^{\prime},b}(\mu_{ab}-\mu_{a^{\prime}b})^{2}\geq\tau$ and $\min_{a,b\neq
b^{\prime}}(\mu_{ab}-\mu_{ab^{\prime}})^{2}\geq\tau$. Then for all
$\varepsilon>0$, for
$(\boldsymbol{g},\boldsymbol{h})\in\mathcal{J}_{\varepsilon}$ and
$(\boldsymbol{C},\boldsymbol{D})\notin\mathcal{P}_{\delta}\cap\mathcal{Q}_{\delta}$,
$G_{M_{0}}(\boldsymbol{C},\boldsymbol{D})$ is small in the sense that

$G_{M_{0}}(\boldsymbol{C},\boldsymbol{D})-\sum_{a,b}[\boldsymbol{C}\boldsymbol{1}]_{a}[\boldsymbol{D}\boldsymbol{1}]_{b}f([\boldsymbol{M}_{0}]_{ab})\leq-\frac{\tau\varepsilon^{2}\delta}{4\min\\{K^{2},L^{2}\\}}.$

The proof of Lemma B.3 is similar to the proof of Lemma 4.3 except that the
bound on the difference uses $\varepsilon$ in place of the random value
$\eta$. Some details follow.

###### Proof of Lemma B.3.

First note that in Lemma A.2, we can let $a$ be the index of the largest
element in column $k$ of matrix $\boldsymbol{C}$; then, since we are
restricted to the set $\mathcal{J}_{\varepsilon}$, this element must be at
least as large as

$C_{ak}\geq\frac{[\boldsymbol{C}^{T}\boldsymbol{1}]_{k}}{K}\geq\frac{\varepsilon}{K}.$

Therefore the bound in Lemma A.2 can be rewritten as

$\sum_{k=1}^{K}\sum_{l=1}^{L}[\boldsymbol{C}^{T}\boldsymbol{1}]_{k}[\boldsymbol{D}^{T}\boldsymbol{1}]_{l}f\bigg{(}\frac{[\boldsymbol{C}^{T}\boldsymbol{M}\boldsymbol{D}]_{kl}}{[\boldsymbol{C}^{T}\boldsymbol{1}]_{k}[\boldsymbol{D}^{T}\boldsymbol{1}]_{l}}\bigg{)}\\\
\leq\sum_{a=1}^{A}\sum_{b=1}^{B}[\boldsymbol{C}\boldsymbol{1}]_{a}[\boldsymbol{D}\boldsymbol{1}]_{b}f(\mu_{ab})-\frac{C\varepsilon^{2}\delta}{K^{2}}.$

Noting that for the Gaussian rate function $f^{\prime\prime}(\mu)=1$ for all
$\mu\in\mathcal{M}$, the remainder of the proof is similar to the proof of
Lemma 4.3. ∎

We establish a finite sample bound by combining these results.

###### Proof of Theorem 3.2.

Fix $\delta>0$ and define $\mathcal{P}_{\delta}$ and $\mathcal{Q}_{\delta}$ as
in Lemma 4.3. Set
$r_{n}=\sup_{\mathcal{J}_{\varepsilon}}|F(\boldsymbol{g},\boldsymbol{h})-G_{\boldsymbol{M}_{0}}(\boldsymbol{C}(\boldsymbol{g}),\boldsymbol{D}(\boldsymbol{h})|$.
Suppose $(\boldsymbol{g},\boldsymbol{h})\in\mathcal{J}_{\varepsilon}$. In this
case,

$\displaystyle
F(\boldsymbol{g},\boldsymbol{h})-F(\boldsymbol{c},\boldsymbol{d})$
$\displaystyle\leq
2r_{n}+\\{G_{\boldsymbol{M}_{0}}(\boldsymbol{C}(\boldsymbol{g}),\boldsymbol{D}(\boldsymbol{h}))-G_{\boldsymbol{M}_{0}}(\boldsymbol{C}(\boldsymbol{c}),\boldsymbol{D}(\boldsymbol{d}))\\}$
$\displaystyle=2r_{n}+\big{\\{}G_{\boldsymbol{M}_{0}}(\boldsymbol{C}(\boldsymbol{g}),\boldsymbol{D}(\boldsymbol{h}))-\sum_{a,b}[\boldsymbol{C}\boldsymbol{1}]_{a}[\boldsymbol{D}\boldsymbol{1}]_{b}f([\boldsymbol{M}_{0}]_{ab})\big{\\}}.$

Applying Lemma B.3 to the second term in the inequality, we get that

$F(\boldsymbol{g},\boldsymbol{h})-F(\boldsymbol{c},\boldsymbol{d})\leq
2r_{n}-\frac{\tau\varepsilon^{2}\delta}{4\min\\{K^{2},L^{2}\\}}$

for all $(\boldsymbol{g},\boldsymbol{h})\in\mathcal{J}_{\varepsilon}$ such
that
$(\boldsymbol{C},\boldsymbol{D})\notin\mathcal{P}_{\delta}\cap\mathcal{Q}_{\delta}$.
The result follows by applying Lemma B.2. ∎

##  Appendix C Additional empirical results

This appendix reports additional empirical results for Bernoulli, Gaussian,
and non-standardized Student’s t distributed data. Figures 3-5 present the
average bicluster misclassification rates for each sample size and Tables 3-5
report the standard deviations for the Bernoulli, Gaussian, and t simulations,
respectively. Since the normalization for the DI-SIM algorithm is only
sbpecified for non-negative data, the algorithm is run on the un-normalized
matrix for the Gaussian and non-standardized Student’s t examples.

For the Bernoulli simulation, we simulate from a block model with $K=2$ row
clusters and $L=3$ column clusters. We vary the number of columns, $n$,
between 200 to 1500 and we take the number of rows as $m=\gamma n$ where
$\gamma\in\\{0.5,1,2\\}$.

We set the row and column class membership probabilities as
$\boldsymbol{p}=(0.3,0.7)$ and $\boldsymbol{q}=(0.2,0.3,0.5)$. We choose the
matrix of block parameters to be

$\boldsymbol{M}=\frac{b}{\sqrt{n}}\begin{pmatrix}0.43&0.06&0.13\\\
0.10&0.34&0.17\end{pmatrix}.$

where the entries were selected to be on the same scale as Bickel and Chen
(2009). We vary $b$ between 5 and 20. We generate the data conditional on the
row and column classes as
$X_{ij}\mid\boldsymbol{c},\boldsymbol{d}\sim\mathrm{Bernoulli}(\mu_{c_{i}d_{j}}).$

Figure 3: Average misclassification rates for Bernoulli example over 100
simulations. Table 3: Standard deviations for Bernoulli example over 100
simulations. $m=0.5n$
---
| PL | KM | DS
n | b=5 | b=10 | b=20 | b=5 | b=10 | b=20 | b=5 | b=10 | b=20
200 | 0.1474 | 0.0281 | 0.0344 | 0.1230 | 0.1490 | 0.0982 | 0.0614 | 0.0503 | 0.0408
500 | 0.0471 | 0.0096 | 0.0029 | 0.1217 | 0.1490 | 0.1399 | 0.0321 | 0.0205 | 0.0494
1000 | 0.0111 | 0.0052 | 0.0010 | 0.1573 | 0.1426 | 0.1490 | 0.0223 | 0.0123 | 0.1075
1500 | 0.0073 | 0.0030 | 0.0005 | 0.1529 | 0.1720 | 0.1688 | 0.0189 | 0.0087 | 0.1229
$m=n$
n | b=5 | b=10 | b=20 | b=5 | b=10 | b=20 | b=5 | b=10 | b=20
200 | 0.0889 | 0.0113 | 0.0285 | 0.1234 | 0.1357 | 0.1112 | 0.0526 | 0.0248 | 0.0657
500 | 0.0121 | 0.0376 | 0.0003 | 0.1466 | 0.1483 | 0.1544 | 0.0233 | 0.0391 | 0.0905
1000 | 0.0056 | 0.0013 | 0.0001 | 0.1659 | 0.1687 | 0.1666 | 0.0132 | 0.0956 | 0.1267
1500 | 0.0032 | 0.0403 | 0.0000 | 0.1675 | 0.1720 | 0.1641 | 0.0090 | 0.0637 | 0.0923
$m=2n$
n | b=5 | b=10 | b=20 | b=5 | b=10 | b=20 | b=5 | b=10 | b=20
200 | 0.0275 | 0.0321 | 0.0000 | 0.1336 | 0.1611 | 0.1397 | 0.0373 | 0.0756 | 0.0998
500 | 0.0048 | 0.0009 | 0.0000 | 0.1598 | 0.1668 | 0.1481 | 0.0120 | 0.1112 | 0.1166
1000 | 0.0015 | 0.0001 | 0.0000 | 0.1590 | 0.1652 | 0.1485 | 0.0543 | 0.1230 | 0.1232
1500 | 0.0009 | 0.0000 | 0.0000 | 0.1591 | 0.1679 | 0.1617 | 0.0746 | 0.0859 | 0.0854

For the Gaussian simulation, we simulate from a block model with $K=2$ row
clusters and $L=3$ column clusters. We vary the number of columns, $n$,
between 50 to 400 and we take the number of rows as $m=\gamma n$ where
$\gamma\in\\{0.5,1,2\\}$.

We set the row and column class membership probabilities as
$\boldsymbol{p}=(0.3,0.7)$ and $\boldsymbol{q}=(0.2,0.3,0.5)$. We choose the
matrix of block parameters to be

$\boldsymbol{M}=b\begin{pmatrix}\phantom{-}0.47&0.15&-0.60\\\
-0.26&0.82&\phantom{-}0.80\end{pmatrix}$

where the entries were simulated from a uniform distribution on $[-1,1]$. We
vary $b$ between 0.5 and 2. We generate the data conditional on the row and
column classes as
$X_{ij}\mid\boldsymbol{c},\boldsymbol{d}\sim\mathrm{Gaussian}(\mu_{c_{i}d_{j}},\sigma=1).$

Figure 4: Average misclassification rates for Gaussian example over 500
simulations. Table 4: Standard deviations for Gaussian example over 500
simulations. $m=0.5n$
---
| PL | KM | DS
n | b=5 | b=10 | b=20 | b=5 | b=10 | b=20 | b=5 | b=10 | b=20
50 | 0.1300 | 0.0874 | 0.1136 | 0.1019 | 0.1068 | 0.1247 | 0.0961 | 0.1890 | 0.1781
100 | 0.0681 | 0.0808 | 0.1094 | 0.0737 | 0.1074 | 0.1116 | 0.1264 | 0.1750 | 0.1975
200 | 0.0406 | 0.0958 | 0.1077 | 0.0700 | 0.0966 | 0.1122 | 0.1699 | 0.1846 | 0.2108
400 | 0.0459 | 0.1074 | 0.1064 | 0.0754 | 0.1071 | 0.1106 | 0.1274 | 0.2111 | 0.2151
$m=n$
n | b=5 | b=10 | b=20 | b=5 | b=10 | b=20 | b=5 | b=10 | b=20
50 | 0.1121 | 0.0830 | 0.1181 | 0.1010 | 0.1156 | 0.1234 | 0.1076 | 0.1940 | 0.1865
100 | 0.0576 | 0.0964 | 0.1047 | 0.0690 | 0.1002 | 0.1121 | 0.1576 | 0.1727 | 0.2060
200 | 0.0576 | 0.1120 | 0.1137 | 0.0947 | 0.1119 | 0.1204 | 0.1516 | 0.2033 | 0.2106
400 | 0.0843 | 0.1102 | 0.1124 | 0.0819 | 0.1096 | 0.1281 | 0.1712 | 0.2144 | 0.2147
$m=2n$
n | b=5 | b=10 | b=20 | b=5 | b=10 | b=20 | b=5 | b=10 | b=20
50 | 0.0829 | 0.1014 | 0.1171 | 0.0858 | 0.1290 | 0.1305 | 0.1194 | 0.1604 | 0.1893
100 | 0.0527 | 0.1215 | 0.1255 | 0.0851 | 0.1212 | 0.1397 | 0.1696 | 0.1956 | 0.2050
200 | 0.0786 | 0.1126 | 0.1126 | 0.0911 | 0.1146 | 0.1364 | 0.1578 | 0.2084 | 0.2106
400 | 0.0848 | 0.0966 | 0.0995 | 0.0837 | 0.0992 | 0.1410 | 0.2010 | 0.2118 | 0.2139

For the non-standardized Student’s t simulation, we use the same parameters as
in the Gaussian simulation and we generate the data conditional on the row and
column classes as
$X_{ij}\mid\boldsymbol{c},\boldsymbol{d}\sim\mathrm{t}(\mu_{c_{i}d_{j}},\sigma=1)$
with four degrees of freedom.

Figure 5: Average misclassification rates for t example over 500 simulations.
Table 5: Standard deviations for t example over 500 simulations. $m=0.5n$
---
| PL | KM | DS
n | b=5 | b=10 | b=20 | b=5 | b=10 | b=20 | b=5 | b=10 | b=20
50 | 0.0039 | 0.0059 | 0.0042 | 0.0033 | 0.0050 | 0.0059 | 0.0034 | 0.0061 | 0.0080
100 | 0.0056 | 0.0029 | 0.0042 | 0.0040 | 0.0041 | 0.0051 | 0.0043 | 0.0065 | 0.0078
200 | 0.0031 | 0.0029 | 0.0054 | 0.0028 | 0.0047 | 0.0055 | 0.0055 | 0.0060 | 0.0088
400 | 0.0022 | 0.0040 | 0.0049 | 0.0032 | 0.0046 | 0.0050 | 0.0055 | 0.0080 | 0.0096
$m=n$
n | b=5 | b=10 | b=20 | b=5 | b=10 | b=20 | b=5 | b=10 | b=20
50 | 0.0047 | 0.0045 | 0.0041 | 0.0036 | 0.0044 | 0.0059 | 0.0037 | 0.0064 | 0.0074
100 | 0.0044 | 0.0034 | 0.0049 | 0.0039 | 0.0049 | 0.0054 | 0.0052 | 0.0066 | 0.0084
200 | 0.0025 | 0.0044 | 0.0054 | 0.0027 | 0.0053 | 0.0052 | 0.0062 | 0.0068 | 0.0092
400 | 0.0026 | 0.0049 | 0.0050 | 0.0045 | 0.0055 | 0.0050 | 0.0048 | 0.0092 | 0.0096
$m=2n$
n | b=5 | b=10 | b=20 | b=5 | b=10 | b=20 | b=5 | b=10 | b=20
50 | 0.0057 | 0.0037 | 0.0047 | 0.0043 | 0.0072 | 0.0058 | 0.0045 | 0.0101 | 0.0081
100 | 0.0034 | 0.0032 | 0.0048 | 0.0032 | 0.0057 | 0.0052 | 0.0060 | 0.0068 | 0.0090
200 | 0.0034 | 0.0039 | 0.0055 | 0.0034 | 0.0047 | 0.0054 | 0.0063 | 0.0073 | 0.0093
400 | 0.0041 | 0.0039 | 0.0050 | 0.0050 | 0.0045 | 0.0052 | 0.0060 | 0.0087 | 0.0095

Similar to the Poisson simulation, biclustering based on the profile log-
likelihood criterion performs at least as well as the other methods and shows
signs of convergence in all three examples. These results provide further
verification of the theoretical findings and support the use of biclustering
based on the profile log-likelihood criterion.

##  Appendix D Additional Application - Movielens

The MovieLens data set consists of 100,000 movie reviews on 1682 movies by 943
users. Each user has rated at least 20 movies and each movie is rated on a
scale from one to five. In addition to the review rating, the release date and
genre of each movie is available as well as some demographic information about
each user including gender, age, occupation and zip code.

In order to retain customers, movie-renting services strive to recommend new
movies to individuals who are likely to view them. Since most users only
review movies that they have already seen, we can use the structure of the
user-movie review matrix to identify associations between users and viewing
habits of movies. Specifically, we consider the 943$\times$1682 binary matrix
$\boldsymbol{X}$ where $X_{ij}=1$ if user $i$ has rated movie $j$ and
$X_{ij}=0$ otherwise. To find structure in $\boldsymbol{X}$, we biclustered
the rows and columns of $\boldsymbol{X}$ using the profile log-likelihood
based on the Bernoulli criterion (5.1a).

To determine a reasonable selection for the number of biclusters we varied the
number of user groups from two to three and the number of movie groups from
two to seven. Qualitatively, the model with three user groups and six movie
groups seems to provide a parsimonious description of the data.

Figure 6 presents the heatmap of the data based on the resulting bicluster
assignments, with the ordering of the clusters determined by the total number
of a reviews in each cluster. Table 6 reports the top ten movies in each
group. The eclectic mix of genres within each movie group suggests that the
rating behavior of users is not explained by genre alone. Figure 7 presents a
boxplot comparing the distributions of the movie release years for each group.
We can see a clear ordering of the movie groups by median release date.

Figure 6: Heatmap generated from MovieLens data reflecting the varying review
patterns in the different biclusters. Blue identfies movies with no review and
white identifies rated movies.

The median ages within the user group were 32, 31, and 30.5, and the
percentages of male users within each group were 67.6%, 72.8%, and 77.8%.
These statistics suggest that there is some age and gender effect on the
reviewing habits of the users.

Table 6: The top ten movies in each cluster based on the total number of reviews. Group 1 | Group 2
---|---
Jade (1995) | She’s the One (1996)
When the Cats Away (1996) | Jack (1996)
Jaws 3-D (1983) | The Preacher’s Wife (1996)
Bastard Out of Carolina (1996) | Striptease (1996)
Exit to Eden (1994) | Mirror Has Two Faces, The (1996)
The Ruling Class (1972) | Hercules (1997)
The Air Up There (1994) | Kids in the Hall: Brain Candy (1996)
Bad Taste (1987) | Jean de Florette (1986)
Stuart Saves His Family (1995) | The Fan (1996)
Cabin Boy (1994) | Extreme Measures (1996)
Group 3 | Group 4
The Firm (1993) | Courage Under Fire (1996)
The Abyss (1989) | Volcano (1997)
Die Hard: With a Vengeance (1995) | Murder at 1600 (1997)
Remains of the Day, The (1993) | Mars Attacks! (1996)
Sneakers (1992) | The People vs. Larry Flynt (1996)
The Professional (1994) | Starship Troopers (1997)
Clerks (1994) | Eraser (1996)
Reservoir Dogs (1992) | Das Boot (1981)
Like Water For Chocolate (1992) | Good Will Hunting (1997)
Chinatown (1974) | The Fifth Element (1997)
Group 5 | Group 6
Raiders of the Lost Ark (1981) | Star Wars (1977)
Pulp Fiction (1994) | Contact (1997)
The Silence of the Lambs (1991) | Fargo (1996)
The Empire Strikes Back (1980) | Return of the Jedi (1983)
Back to the Future (1985) | Liar Liar (1997)
The Fugitive (1993) | The English Patient (1996)
Indiana Jones and the Last Crusade (1989) | Scream (1996)
The Princess Bride (1987) | Toy Story (1995)
Forrest Gump (1994) | Air Force One (1997)
Monty Python and the Holy Grail (1974) | Independence Day (1996)
Figure 7: Boxplot comparing the different clusters based on movie release
dates.

Roughly speaking, user group 3 is consistently active across all movie groups
with increasing activity as the popularity of the movie increases. The
reviewing habits of user group 2 follow a similar pattern but to a lesser
extent. In contrast, user group 1 is consistently inactive with the only
exceptions being movie groups 4 and 6. From Figure 7, it appears that the
users in group 1 only rate recent movies whereas users in groups 2 and 3 rate
movies from all time periods.

The biclusters discovered here suggest that a movie-renting service should
recommend under-reviewed movies to individuals in user group 3, and it should
recommend new releases to users in group 1. This information can be used in an
ensemble-based recommendation engine like that of Töscher, Jahrer and Bell
(2009).

## References

  * Arabie, Boorman and Levitt (1978) [author] Arabie, PhippsP., Boorman, Scott AS. A. and Levitt, Paul RP. R. (1978). Constructing blockmodels: How and why. J. Math. Pscyh. 17 21–63.
  * Bickel and Chen (2009) [author] Bickel, P. J.P. J. and Chen, A.A. (2009). A Nonparametric View of Network Models and Newman-Girvan and Other Modularities. Proc. Nat. Acad. Sci. USA 106 21068–21073.
  * Cheng and Church (2000) [author] Cheng, Y.Y. and Church, G. M.G. M. (2000). Biclustering of expression data. Proceedings International Conference on Intelligent Systems for Molecular Biology ; ISMB. International Conference on Intelligent Systems for Molecular Biology 8 93–103.
  * Choi, Wolfe and Airoldi (2012) [author] Choi, D.D., Wolfe, P. J.P. J. and Airoldi, E. M.E. M. (2012). Stochastic blockmodels with growing number of classes. Biometrika 99 273–284.
  * Dembo and Zeitouni (1998) [author] Dembo, A.A. and Zeitouni, O.O. (1998). Large Deviations Techniques and Applications, Second ed. Springer-Verlag.
  * Dhillon (2001) [author] Dhillon, I. S.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) 26–29.
  * Eisen et al. (1998) [author] Eisen, Michael B.M. B., Spellman, Paul T.P. T., Brown, Patrick O.P. O. and Botstein, DavidD. (1998). Cluster analysis and display of genome-wide expression patterns. Proc. Nat. Acad. Sci. USA 95 14863–14868.
  * Getz, Levine and Domany (2000) [author] Getz, G.G., Levine, E.E. and Domany, E.E. (2000). Coupled Two-Way Clustering Analysis of Gene Microarray Data. Proc. Nat. Acad. Sci. USA 97 12079–12084.
  * Golub and Loan (1996) [author] Golub, Gene H.G. H. and Loan, Charles F. VanC. F. V. (1996). Matrix Computation. Johns Hopkins University Press.
  * GroupLens Research Project (1998) [author] GroupLens Research Project, (1998). MovieLens Dataset. University of Minnesota http://www.grouplens.org/data/.
  * Harpaz et al. (2010) [author] Harpaz, R.R., Perez, H.H., Chase, H. S.H. S., Rabadan, R.R., Hripcsak, G.G. and Friedman, C.C. (2010). Biclustering of Adverse Drug Events in the FDA’s Spontaneous Reporting System. Clinical Pharmacology & Therapeutics 89 243–250.
  * Hartigan (1972) [author] Hartigan, J. A.J. A. (1972). Direct Clustering of a Data Matrix. J. Amer. Statist. Assoc. 67 123–129.
  * Hofmann (1999) [author] Hofmann, T.T. (1999). Latent Class Models for Collaborative Filtering. In In Proceedings of the Sixteenth International Joint Conference on Artificial Intelligence 688–693.
  * Holland, Laskey and Leinhardt (1983) [author] Holland, P. W.P. W., Laskey, K. B.K. B. and Leinhardt, S.S. (1983). Stochastic Blockmodels: First Steps. Social Networks 5 109–137.
  * Kernighan and Lin (1970) [author] Kernighan, B. W.B. W. and Lin, S.S. (1970). An Efficient Heuristic Procedure for Partitioning Graphs. The Bell system technical journal 49 291–307.
  * Kluger et al. (2003) [author] Kluger, Y.Y., Basri, R.R., Chang, J. T.J. T. and Gerstein, M.M. (2003). Spectral Biclustering of Microarray Data: Coclustering Genes and Conditions. Genome Research 13 703–716.
  * Lazzeroni and Owen (2002) [author] Lazzeroni, L.L. and Owen, A.A. (2002). Plaid models for gene expression data. Statist. Sinica 12 61–86.
  * Madeira and Oliveira (2004) [author] Madeira, S. C.S. C. and Oliveira, A. L.A. L. (2004). Biclustering Algorithms for Biological Data Analysis: A Survey. IEEE T. Comput. Bi. 1 24–45.
  * McSherry (2001) [author] McSherry, FrankF. (2001). Spectral partitioning of random graphs. In Proceedings 2001 IEEE International Conference on Cluster Computing 529–537. IEEE Comput. Soc.
  * Mirkin (1996) [author] Mirkin, BorisB. (1996). Mathematical classification and clustering. Kluwer Academic Press.
  * Murphy and van der Vaart (2000) [author] Murphy, S. A.S. A. and van der Vaart, A. W.A. W. (2000). On Profile Likelihood. J. Amer. Statist. Assoc. 95 449–465.
  * Newman (2006) [author] Newman, M. E. J.M. E. J. (2006). Modularity and community structure in networks. Proc. Nat. Acad. Sci. USA 103 8577–8582.
  * Perry and Owen (2010) [author] Perry, Patrick O.P. O. and Owen, Art B.A. B. (2010). A Rotationn Test to Verify Latent Structure. J. Mach. Learn. Res. 11 603–624.
  * Perry and Wolfe (2012) [author] Perry, P. O.P. O. and Wolfe, P. J.P. J. (2012). Null Models for Network Data. Preprint arXiv:1201.5871.
  * Rohe, Chatterjee and Yu (2011) [author] Rohe, K.K., Chatterjee, S.S. and Yu, B.B. (2011). Spectral Clustering and the High-Dimensional Stochastic Blockmodel. Ann. Statist. 39 1878–1915.
  * Rohe and Yu (2012) [author] Rohe, K.K. and Yu, B.B. (2012). Co-clustering for Directed Graphs; the Stochastic Co-Blockmodel and a Spectral Algorithm. Preprint arXiv:1204.2296.
  * Töscher, Jahrer and Bell (2009) [author] Töscher, A.A., Jahrer, M.M. and Bell, R. M.R. M. (2009). The BigChaos Solution to the Netflix Grand Prize. Technical Report.
  * Ungar and Foster (1998) [author] Ungar, L.L. and Foster, D. P.D. P. (1998). A Formal Statistical Approach to Collaborative Filtering. In In CONALD 98.
  * Varadhan (2001) [author] Varadhan, S. R. S.S. R. S. (2001). Probability Theory (Courant Lecture Notes). American Mathematical Soceity.
  * Zahn et al. (2007) [author] Zahn, Jacob M.J. M., Poosala, SureshS., Owen, Art B.A. B., Ingram, Donald K.D. K., Lustig, AnaA., Carter, ArnellA., Weeraratna, Ashani T.A. T., Taub, Dennis D.D. D., Gorospe, MyriamM., Mazan-Mamczarz, KrystynaK., Lakatta, Edward G.E. G., Boheler, Kenneth R.K. R., Xu, XiangruX., Mattson, Mark P.M. P., Falco, GeppinoG., Ko, Minoru S. H.M. S. H., Schlessinger, DavidD., Firman, JeffreyJ., Kummerfeld, Sarah K.S. K., Wood, William H.W. H., Zonderman, Alan B.A. B., Kim, Stuart K.S. K. and Becker, Kevin G.K. G. (2007). AGEMAP: A Gene Expression Database for Aging in Mice. PLOS Genetics.
  * Zhao, Levina and Zhu (2011) [author] Zhao, Y.Y., Levina, E.E. and Zhu, J.J. (2011). Community extraction for social networks. P. Natl. Acad. Sci. USA 108 7321–7326.
  * Zhao, Levina and Zhu (2012) [author] Zhao, YunpengY., Levina, ElizavetaE. and Zhu, JiJ. (2012). Consistency of Community Detection in Networks Under Degree-Corrected Stochastic Block Models. Preprint arXiv:1110.3854.