# Learning Sparsity and Block Diagonal Structure in Multi-View Mixture Models

Iain Carmichael111Department of Statistics, University of Washington, Seattle,
WA

###### Abstract

Scientific studies increasingly collect multiple modalities of data to
investigate a phenomenon from several perspectives. In integrative data
analysis it is important to understand how information is heterogeneously
spread across these different data sources. To this end, we consider a
parametric clustering model for the subjects in a multi-view data set (i.e.
multiple sources of data from the same set of subjects) where each view
marginally follows a mixture model. In the case of two views, the dependence
between them is captured by a cluster membership matrix parameter and we aim
to learn the structure of this matrix (e.g. the zero pattern). First, we
develop a penalized likelihood approach to estimate the sparsity pattern of
the cluster membership matrix. For the specific case of block diagonal
structures, we develop a constrained likelihood formulation where this matrix
is constrained to be block diagonal up to permutations of the rows and
columns. To enforce block diagonal constraints we propose a novel optimization
approach based on the symmetric graph Laplacian. We demonstrate the
performance of these methods through both simulations and applications to data
sets from cancer genetics and neuroscience. Both methods naturally extend to
multiple views.

Keywords: Multi-view data, integrative clustering, graph Laplacian, structured
sparsity, EM-algorithm, model-based clustering, TCGA, neuron cell type

##  1 Introduction

Scientific studies often investigate a phenomenon from several perspectives by
collecting multiple modalities of data. For example, modern cancer studies
collect data from several genomic platforms such as RNA expression, microRNA,
DNA methylation and copy number variations (Network et al.,, 2012; Hoadley et
al.,, 2018). Neuroscientists investigate neurons using transcriptomic,
electrophysiological and morphological measurements (Tasic et al.,, 2018;
Gouwens et al.,, 2019, 2020). These integrative studies require methods to
analyze multi-view data: a fixed set of observations with several disjoint
sets of variables (views).

Classical multi-view methods such as canonical correlation analysis for
dimensionality reduction estimate joint information shared by all views
(Hotelling,, 1936). Similarly, many multi-view clustering methods assume there
is one consensus clustering (see Figure 2(d) below) that is present in each
data-view (Shen et al.,, 2009; Kumar et al.,, 2011; Kirk et al.,, 2012; Lock
and Dunson,, 2013; Gabasova et al.,, 2017; Wang and Allen,, 2019). A singular
focus on joint signals ignores the possibility that information is
heterogeneously spread across the views. For example, environmental factors
might show up in a clinical data-view, but not in a genomic data-view.
Contemporary multi-view methods examine how information is shared (or not
shared) by different views. Recent work in dimensionality reduction looks for
partially shared latent signals (Lock et al.,, 2013; Klami et al.,, 2014; Zhao
et al.,, 2016; Gaynanova and Li,, 2017; Feng et al.,, 2018). Similarly, recent
multi-view clustering methods investigate how clustering information is spread
across multiple views (Hellton and Thoresen,, 2016; Gao et al., 2019a, ; Gao
et al., 2019b, ).

A motivating example comes from breast cancer pathology where investigators
study tumors using both genomic and histological222Meaning a doctor or
algorithm visually examines an image of a tumor biopsy. information
(Carmichael et al.,, 2019). Breast cancer tumor subtypes can be defined using
either genomic information (e.g. the PAM50 molecular subtypes Parker et al.,
2009) or histological information (e.g. high, medium or low grade Hoda et al.,
2020). Some cluster information may be jointly shared by both data views e.g.
if histological subtype 1 corresponds to exactly genomic subtype 1. Other
information may be contained in one view but not another view e.g. if
histological subtype 2 correspond to genomic subtypes 2 and 3. See Figure 1(a)
below.

We develop an approach to learn how information is spread across views in a
multi-view mixture model (MVMM) (Bickel and Scheffer,, 2004). This model,
detailed in Section 2, makes two assumptions for a $V\geq 2$ view data set:

  1. 1.

Marginally, each view follows a mixture model i.e. there are $V$ sets of view-
specific clusters.

  2. 2.

The views are independent given the marginal view cluster memberships.

We further assume there may be some kind of “interesting relationship” between
clusters in different views. For example, in a two-view data set every
observation has two (hidden) cluster labels
$(y^{(1)},y^{(2)})\in[K^{(1)}]\times[K^{(2)}]$ where $K^{(v)}$ is the number
of clusters in the $v$th view and $[K]:=\\{1,\dots,K\\}$. The joint
distribution of the cluster labels is described by the cluster membership
probability matrix $\pi\in\mathbb{R}^{K^{(1)}\times K^{(2)}}_{+}$ where

$\pi_{k^{(1)},k^{(2)}}=P(y^{(1)}=k^{(1)},y^{(2)}=k^{(2)}),\text{ for
}k^{(1)}\in[K^{(1)}]\text{ and }k^{(2)}\in[K^{(2)}].$

The structure of this matrix captures how information is shared between the
two views. Figure 1(a) shows a hypothetical $\pi$ matrix. Many of the entries
are zero, meaning, for instance, an observation cannot be simultaneously in
cluster 1 in the first view and cluster 2 in the second view. In this example,
cluster 1 in the first view is exactly the same as cluster 1 in the second
view; this information is shared jointly by both views. On the other hand,
cluster 3 in the second view breaks up into clusters 3, 4 and 5 in the first
view; here there is information in the first view that is not contained in the
second view. In general $\pi$ may be anywhere from rank 1 (i.e. the views are
independent thus share no information) to diagonal (the consensus clustering
case where the views contain the same information). The goal of this paper is
to learn the structure of $\pi$ while simultaneously learning the cluster
parameters (e.g. cluster means).

(a) An example block diagonal $\pi$ matrix.

(b) The bipartite graph whose node sets are the rows and columns of $\pi$ and
whose edges the entries of $\pi$.

Figure 1:  The matrix $\pi$ (Figure 1(a)) captures the between-view
dependence. $\pi$ can equivalently be thought of as a bipartite graph (Figure
1(b)). The connected components of this graph represent “clusters of clusters”
that are related to each other. Note $\pi$ can be rectangular in general.

Section 2 formalizes the multi-view mixture model outlined above. Section 3
presents a penalized likelihood approach making use of the concave
$\log(\cdot+\delta)$ penalty to estimate the zero-pattern of $\pi$. Section 4
considers the case when $\pi$ has block diagonal structure and formulates a
block diagonally constrained maximum likelihood version of the MVMM. This
section develops an alternating minimization approach for imposing block
diagonal matrix constraints in general via the symmetric Laplacian. A detailed
discussion of this alternating algorithm and convergence results are provided
in Section C. An extension of this approach to block diagonal multi-arrays is
sketched in Section A. Section 5 presents a simulation study of the methods
developed in this paper. Section 6 applies these methods to the TCGA breast
cancer data set and an excitatory mouse neuron data set. The main algorithmic
ideas are presented in the body of the paper and detailed discussions are
provided in the appendix. Proofs and additional simulations are also provided
in the appendix.

The methods developed in this paper are implemented in a publicly available
python package www.github.com/idc9/mvmm. Code to reproduce the simulations as
well as supplementary data and figures can be found at
www.github.com/idc9/mvmm_sim. The code makes use of the following python
packages: Hunter, (2007); McKinney et al., (2010); Walt et al., (2011);
Pedregosa et al., (2011); Diamond and Boyd, (2016); Waskom et al., (2017);
Davidson-Pilon et al., (2020); Virtanen et al., (2020).

###  1.1 Summary of contributions and related work

We develop two novel methods that explore how information is shared between
views in the parametric multi-view mixture model of Bickel and Scheffer,
(2004). Both methods impose interpretable structure — sparsity (Section 3) or
block diagonal constraints (Section 4.3) — on the cluster membership matrix.
They also lead to challenging optimization issues. Our approaches to address
these challenges are of interest in applications beyond this paper.

Many existing multi-view clustering methods focus on the consensus clustering
case (see reference above). While the consensus clustering case is a special
case of the MVMM when $\pi$ is diagonal, our method allows for more flexible
relations among the clusters in each view. The work of Gao et al., 2019a ; Gao
et al., 2019b  takes an important step beyond consensus clustering by
developing a test for independence between the views in a two-view MVMM.

The joint and individual clustering (JIC) method developed by Hellton and
Thoresen, (2016) is a multi-view clustering algorithm based on dimensionality
reduction using JIVE (Lock et al.,, 2013). JIC identifies information that is
either shared by all views (joint clusters) or is only contained in one view
(individual clusters). An immediate difference between JIC and our methods is
that we work with parametric mixture models while JIC is based on
dimensionality reduction. Moreover, our methods take a different perspective
on how information can be shared among views (see Footnote 3).

The $\log(\cdot+\delta)$ penalized likelihood approach adopted in Section 3
was developed in Huang et al., (2017) for (single view) mixture-model model
selection. To fit mixture models with this penalty Huang et al., (2017)
suggests an EM algorithm where the M-step is approximated with a soft-
thresholding operation. This soft-thresholding approximation — based on a
heuristic argument — is used by a number of other papers (Yao et al.,, 2018;
Yu and Wang,, 2019; Bugdary and Maymon,, 2019) and similar approximations
appear elsewhere (Gu and Xu,, 2019). We provide rigorous justification for
this soft-thresholding approximation and show the algorithm is insensitive to
the choice of $\delta$ for small values of $\delta$ (Theorem 3.1).

The task of learning model parameters with block diagonal structure arises in
a variety of contexts including: graphical models (Marlin and Murphy,, 2009;
Tan et al.,, 2015; Devijver and Gallopin,, 2018; Kumar et al.,, 2019), co-
clustering (Han et al.,, 2017; Nie et al.,, 2017), subspace clustering (Feng
et al.,, 2014; Lu et al.,, 2018), principal components analysis (Asteris et
al.,, 2015), and community detection (Nie et al.,, 2016). Learning parameter
values and block diagonal structure simultaneously is a combinatorial problem
that is generally intractable except in certain special cases (Asteris et
al.,, 2015). Block diagonal constraints are often enforced with continuous
optimization approaches using the unnormalized graph Laplacian (Nie et al.,,
2016, 2017).

Sections 4 and C develop an approach to impose block diagonal constraints via
the symmetric graph Laplacian. This approach avoids the strong modeling
assumptions — that the row and column sums are known ahead of time — required
by the unnormalized Laplacian (Nie et al.,, 2016, 2017). By making use of an
extremal characterization of generalized eigenvalues we provide an alternating
algorithm for the penalized symmetric Laplacian Problem (14) that is no more
computationally burdensome than the analogous problem with the unnormalized
Laplacian (see Section D).

###  1.2 Notation

A multi-view random vector $x\in\mathbb{R}^{\sum_{v=1}^{V}d^{(v)}}$ is a
random vector where the variables have been partitioned into $V$ mutually
exclusive sets of sizes $d^{(1)},\dots,d^{(V)}$. We write
$x^{(v)}\in\mathbb{R}^{d^{(v)}}$ for the $v$th view i.e. $x$ is the
concatenation of the $x^{(1)},\dots,x^{(V)}$. We use superscript parenthesis,
e.g. $x^{(v)}$, to reference quantities related to a particular view.

For a matrix $V\in\mathbb{R}^{R\times C}$ let $V(r,:)\in\mathbb{R}^{C}$ denote
the $r$th row and let $V(:,c)=V_{c}\in\mathbb{R}^{R}$ denote the $c$th column.
For $v\in\mathbb{R}^{n}$, let $\text{diag}(v)\in\mathbb{R}^{n\times n}$ be the
diagonal matrix whose diagonal elements are given by $v$. Let
$\mathbf{1}_{n}\in\mathbb{R}^{n}$ be the vector of ones. For a set
$A\subseteq[n]$ let $\mathbf{1}_{A}\in\\{0,1\\}^{n}$ denote the vector with 1s
in the entries corresponding to elements of $A$ and 0s elsewhere. The
indicator function, $I:\mathbb{R}^{n}\to\\{0\\}\cup\infty$, of a set
$\mathcal{C}\subseteq\mathbb{R}^{n}$ is defined by $I(x)=0$ if
$x\in\mathcal{C}$ and $I(x)=\infty$ if $x\not\in\mathcal{C}$. For vectors
$a,b$ let $a\odot b$ denote the Haadamard (element-wise) product.

For a symmetric matrix $A\in\mathbb{R}^{n\times n}$ we write
$\lambda_{1}(A)\geq\lambda_{2}(A)\geq\dots$ for the eigenvalues sorted in
decreasing order and $\lambda_{(1)}(A)\leq\lambda_{(2)}(A)\leq\dots$ for the
eigenvalues sorted in increasing order. For two symmetric matrices
$A,B\in\mathbb{R}^{n\times n}$ we write
$\lambda_{1}(A,B)\geq\lambda_{2}(A,B)\geq\dots$ for the generalized
eigenvalues (i.e. numbers $\lambda$ where there exists a $v\in\mathbb{R}^{n}$
such that $Av=\lambda Bv$ with the normalization $v^{T}Bv=1$).

##  2 Multi-view mixture model specification

This section describes a multi-view mixture model for $V\geq 2$ views (Bickel
and Scheffer,, 2004; Gao et al., 2019a, ). This model assumes that marginally,
each view follows a mixture model and that the views are conditionally
independent given the view cluster memberships.

In detail, let $x^{(v)}\in\mathbb{R}^{d^{(v)}}$ denote the random vector for
the $v$th view. In the $v$th view there are $K^{(v)}$ view-specific clusters
and let $y^{(v)}\in[K^{(v)}]$ denote latent, view specific membership
assignment for the $v$th view. Let

$f(x^{(v)}|y^{(v)}=k)=\phi^{(v)}(x|\Theta^{(v)}_{k})\text{ for
}k\in[K^{(v)}],v\in[V],$

be the conditional distribution of the $k$th cluster in the $v$th view where
$\phi^{(v)}(\cdot|\theta)$ is a density function with parameter $\theta$ (e.g.
cluster means). Also let

$P\left(y=(k^{(1)},\dots,k^{(V)})\right)=\pi_{k^{(1)},\dots,k^{(V)}}\text{ for
}k^{(v)}\in[K^{(v)}],v\in[V],$

be the joint distribution of the view specific labels where
$y=(y^{(1)},\dots,y^{(V)})\in\mathbb{Z}^{V}_{+}$ is the latent cluster
membership vector and $\pi\in\mathbb{R}^{K^{(1)}\times\dots\times K^{(V)}}$ is
the cluster membership probability multi-array (non-negative entries summing
to 1). Then the probability density function of the joint distribution is

$f(x,y=(k^{(1)},\dots,k^{(V)})|\Theta,\pi)=\pi_{k^{(1)},\dots,k^{(V)}}\prod_{v=1}^{V}\phi^{(v)}(x^{(v)}|\Theta^{(v)}_{k^{(v)}}),$
(1)

where $\Theta:=\\{\\{\Theta^{(v)}_{k}\\}_{k=1}^{K^{(v)}}\\}_{v=1}^{V}$ is the
collection of view specific cluster parameters. We further assume that the
marginal view-specific cluster probabilities are strictly positive, i.e.

$0<\pi^{(v)}_{\textbf{k}}:=P(y^{(v)}=\textbf{k})=\sum_{j^{(1)}=1}^{K^{(1)}}\dots\sum_{j^{(v-1)}=1}^{K^{(v-1)}}\sum_{j^{(v+1)}=1}^{K^{(v+1)}}\dots\sum_{j^{(V)}=1}^{K^{(V)}}\pi_{j^{(1)},\dots,j^{(v-1)},\textbf{k},j^{(v+1)}\dots
j^{(V)}}$ (2)

for each $k\in[K^{(v)}],\text{ and }v\in[V]$.

The marginal distribution of the $v$th view, $x^{(v)}$, is a mixture model
with $K^{(v)}$ view-specific clusters (Figure 2(a)). The joint distribution,
$x$, is a mixture model with
$|\text{supp}(\pi)|\in[\min_{v\in[V]}(K^{(v)}),\prod_{v=1}^{V}K^{(v)}]$
overall clusters (Figures 2(b)-2(d)). In other words, looking at the joint
distribution there is one set333 In the JIC model the view joint distribution
has $V+1$ sets of clusters for a $V$-view data set — one set of joint clusters
and $V$ sets of view-individual clusters. For details see (Hellton and
Thoresen,, 2016). of $|\text{supp}(\pi)|$ clusters, but the clusters share
parameters.

######  Remark 2.1.

This model promotes parameter sharing; if $\pi$ is dense, the number of
overall clusters scales multiplicatively (e.g. like $O(K^{V})$) in the number
of view marginal clusters while the number of cluster parameters (e.g. cluster
means) scaled additively (e.g. like $O(VK)$).

Figure 2 shows three scenarios for a $V=2$ view data set. Both views are one
dimensional and marginally follow a Gaussian mixture model (GMM) with
$K^{(1)}=K^{(2)}=10$ clusters (Figure 2(a)). In the first scenario (Figure
2(b)) there is no information shared between the two views; $\pi$ is a rank 1
matrix. In the third scenario (Figure 2(d)) the two views capture the same
information i.e. the clusters in the first view are the same clusters as the
clusters in the second view. Here $\pi$ is a diagonal matrix (after
appropriately permuting the cluster labels). In the second scenario (Figure
2(c)) the two views have partially overlapping information. In this scenario
$\pi$ is the block diagonal matrix shown in Figure 1(a) above.

(a) Marginally, each view follows a mixture model with $10$ clusters.

(b) Scenario 1: Independent views; $100$ overall clusters.

(c) Scenario 2: Partial dependence between views; $16$ overall clusters.

(d) Scenario 3: Strong dependence between views; $10$ overall clusters. This
is the consensus clustering case.

Figure 2:  Three scenarios for the joint distribution of two views. The
marginal distribution of each view is a one dimensional mixture model (Figure
2(a)). In Figure 2(b) all possible combinations of view 1 and view 2 clusters
occur with equal probability. In Figure 2(c) only some combinations of view 1
and view 2 clusters are possible. In Figure 2(d) each cluster in the first
view is matched with exactly one cluster in the second view.

Suppose we are given $n$ samples $\\{x_{i}\\}_{i=1}^{n}$ with
$x_{i}\in\mathbb{R}^{\sum_{v=1}^{V}d^{(v)}}$ from a $V$-view data set and have
specified the number of view-specific clusters $K^{(1)},\dots,K^{(V)}$. If no
additional assumptions are placed on $\pi$, we fit the model by maximizing the
log likelihood of the observed data

$\ell(\\{x_{i}\\}_{i=1}^{n}|\Theta,\pi):=\sum_{i=1}^{n}\log
f(x_{i}|\Theta,\pi)$ (3)

using an EM algorithm (Dempster et al.,, 1977) that is detailed in Section
E.1, where

$f(x|\Theta,\pi):=\sum_{k^{(1)}=1}^{K^{(1)}}\dots\sum_{k^{(V)}=1}^{K^{(V)}}\pi_{k^{(1)},\dots,k^{(V)}}\prod_{v=1}^{V}\phi^{(v)}(x^{(v)}|\Theta^{(v)}_{k^{(v)}})$
(4)

is the probability density function of the observed data. The remainder of
this paper focuses on simultaneously estimating the model parameters,
$\Theta,\pi$, as well as the sparsity structure of $\pi$.

##  3 Sparsity inducing log penalty

This section develops a penalized likelihood approach to estimate the sparsity
structure of $\pi$ that avoids the exponential search space of naive
enumeration. We assume the number of view specific clusters,
$K^{(1)},\dots,K^{(V)}$, have been specified.

Consider fitting a standard, single-view mixture model with a sparsity
inducing penalty $p(\cdot)$ (e.g. Lasso or SCAD) on the entries of the cluster
membership probability vector, $\pi\in\mathbb{R}^{K}_{+}$. This raises several
issues. First, the Lasso penalty is constant since $\pi$ lives on the unit
simplex. Second, exact zeros in $\pi$ give a negative infinity the complete
data log likelihood (1), though this issue does not arise in the observed data
log-likelihood (3). If we use an EM algorithm to maximize the observed data
log-likelihood the M-step involves the following optimization problem

$\displaystyle\underset{\pi\in\mathbb{R}^{K}}{\text{minimize}}$
$\displaystyle-\sum_{k=1}^{K}a_{k}\log(\pi_{k})+\lambda\sum_{k=1}^{K}p(\pi_{k})$
(5) subject to $\displaystyle\pi\geq 0\text{ and }\pi^{T}\mathbf{1}_{K}=1,$

where $a\in\mathbb{R}^{K}_{+}$ is the output of the E-step (i.e. the expected
cluster assignments). The log in the first term of the objective function acts
as a barrier function that prevents the solution from having zeros.

Huang et al., (2017) provides theoretical justification for using the penalty
$p(\cdot)=\log(\delta+\cdot)$ for some small $\delta>0$. The following theorem
further justifies the use of this penalty for small $\delta$ by showing that
we can approximate the solution with a quantity that has exact zeros. This
theorem also leads to a computationally efficient approximation for the M-step
and suggests that the penalty is insensitive to the choice of $\delta$ for
small values of $\delta$.

######  Theorem 3.1.

Let $a_{1},\dots,a_{K}\geq 0$, $\sum_{k=1}^{K}a_{k}=1$, and
$0<\lambda<\frac{1}{K}$. Let $z^{\delta}\in\mathbb{R}_{+}^{K}$ be a solution
of the following problem for fixed $\delta>0$,

$\displaystyle\underset{z}{\textup{minimize}}$
$\displaystyle-\sum_{k=1}^{K}a_{k}\log(z_{k})+\lambda\sum_{k=1}^{K}\log(\delta+z_{k})$
(6) subject to $\displaystyle z\geq 0\text{ and }z^{T}\mathbf{1}_{K}=1.$

Then $\lim_{\delta\to 0}z^{\delta}=z^{0}\in\mathbb{R}^{K}$ where

$z_{k}^{0}:=\frac{(a_{k}-\lambda)_{+}}{\sum_{j=1}^{K}(a_{j}-\lambda)_{+}}\text{
for each }k\in[K].$ (7)

This theorem says that for small $\delta$ the global minimizer of (6) is close
to the normalized soft-thresholding operation (7). The condition
$\lambda<\frac{1}{K}$ guarantees the denominator of (7) is non-zero. The soft-
thresholding approximation presented in this theorem is proposed by Huang et
al., (2017) and used as a heuristic (Yao et al.,, 2018; Yu and Wang,, 2019;
Bugdary and Maymon,, 2019); we prove Theorem 3.1 in Section G.1.

Returning to the MVMM, we consider the following penalized likelihood problem

$\displaystyle\underset{\Theta,\pi}{\text{maximize}}$
$\displaystyle\ell(\\{x_{i}\\}_{i=1}^{n}|\Theta,\pi)-\lambda\left(\sum_{k^{(1)}=1}^{K^{(1)}}\dots\sum_{k^{(V)}=1}^{K^{(V)}}\log(\delta+\pi_{k^{(1)},\dots,k^{(V)}})\right),$
(8)

where $\ell$ is the observed data log likelihood (3) and $\delta>0$ is a small
value. This problem can be solved with an EM algorithm similar to the one
derived for the unpenalized model. The M-step of this EM algorithm solves a
problem in the form of (6). Based on Theorem 3.1 we approximate the M-step
using the normalized soft-thresholding operation. Details for this algorithm
can be found in Section E.2.

##  4 Enforcing block diagonal constraints

This section presents a constrained maximum likelihood approach to estimate
$\pi$ under the restriction that $\pi$ has a block diagonal structure.
Sections 4.1 and 4.2 discuss optimization with block diagonal constraints in a
general setting. Section 4.3 presents the particular case of the multi-view
mixture model.

For a fixed matrix we have to be careful about what “block diagonal” means
i.e. one could argue that the matrix $\text{diag}([1,1,0])$ has either 1, 2,
or 3 blocks. We take the convention that blocks must have at least one non-
zero entry and anything that can be a block is a block; thus
$\text{diag}([1,1,0])$ has 2 blocks. For a matrix $X$ whose rows/columns are
allowed to be permuted we say “$X$ is block diagonal with $NB(X)$ blocks up to
permutations” where

$\text{NB}(X):=\max\\{B|\text{the rows/columns of }X\text{ can be permuted to
create a }B\text{ block, block diagonal matrix}\\}$ (9)

Any permutation of the rows/columns of $X$ which achieves the above maximum is
called a maximally block diagonal permutation.

###  4.1 Spectral characterization of block diagonal matrices

This section gives a spectral characterization of block diagonal matrices up
to permutations. Let $A\in\mathbb{R}^{n\times n}_{+}$ be the adjacency matrix
of a weighted, undirected graph with no self loops. The unnormalized Laplacian
is

$L_{\text{un}}(A):=\text{diag}(\text{deg}(A))-A$ (10)

where $\text{deg}(A):=A\mathbf{1}_{n}\in\mathbb{R}^{n}_{+}$ is the vector of
the vertex degrees (Von Luxburg,, 2007). The symmetric, normalized Laplacian
is

$L_{\text{sym}}(A):=I-\text{diag}(\text{deg}(A))^{-1/2}A\text{
diag}(\text{deg}(A))^{-1/2}.$ (11)

When $\text{deg}(A)$ has zeros, the inverse is taken to be the Moore-Penrose
psueo-inverse thus the diagonal elements of $L_{\text{sym}}(\cdot)$ are always
equal to 1 even when there are degree zero (isolated) vertices.444This
convention is not always followed (Von Luxburg,, 2007), as discussed in
Section G.3. The eigenvalues of the symmetric Laplacian are equal555We have to
be careful when $\text{diag}(\text{deg}(A))$ is non-invertible; this issue is
addressed in Section B. to the generalized eigenvalues of
$(L_{\text{un}}(A),\text{diag}(\text{deg}(A))$.

For $X\in\mathbb{R}^{R\times C}_{+}$ let

$A_{\text{bp}}(X):=\begin{bmatrix}0&X\\\
X^{T}&0\end{bmatrix}\in\mathbb{R}^{(R+C)\times(R+C)}$

be the adjacency matrix of the weighted, bipartite graph $G(X)$ whose edge
weights are given by the entries of $X$ and whose vertex sets are the rows and
columns of $X$ (see Figure 1(b)). Note a row or column of zeros in $X$
corresponds to an isolated vertex in the graph.

Proposition 4.1 shows the connected components of this bipartite graph with at
least two vertices capture the block diagonal structure of $X$ up to
permutations; these connected components are in turn captured by the spectrum
of the symmetric, normalized Laplacian.

######  Proposition 4.1.

The following are equivalent for $1\leq
B+Z_{\text{row}}+Z_{\text{col}}\leq\min(R,C)$

  1. 1.

$X$ is block diagonal up to permutations with $B$ blocks and has
$Z_{\text{row}}$ rows and $Z_{\text{col}}$ columns of zeros.

  2. 2.

$G(X)$ has $B$ connected components with at least two vertices and
$Z_{\text{row}}+Z_{\text{col}}$ isolated vertices.

  3. 3.

$L_{\text{sym}}(A_{\text{bp}}(X))$ has exactly $B$ eigenvalues equal to 0.

  4. 4.

$L_{\text{un}}(A_{\text{bp}}(X))$ has exactly
$B+Z_{\text{row}}+Z_{\text{col}}$ eigenvalues equal to 0.

Additionally, the number of eigenvalues equal to 1 of the symmetric Laplacian
is at least $2\cdot(Z_{\text{row}}+Z_{\text{col}})$.

Section A generalizes this proposition to block diagonal multi-arrays.

Proposition 4.1 shows that the symmetric Laplacian gives more precise control
over the block diagonal structure of a matrix than the unnormalized Laplacian
does. The number of 0 eigenvalues of the symmetric Laplacian is exactly the
number of blocks while the number of zero eigenvalues of the unnormalized
Laplacian only upper bounds the number of blocks (see Figure 3(a)).

(a) $X\in\\{0,1\\}_{+}^{5\times 4}$ is block diagonal with two blocks and one
row of zeros.

(b) Spectrum of $L_{\text{sym}}(A_{\text{bp}}(X))$; two eigenvalues are equal
to zero.

(c) Spectrum of $L_{\text{un}}(A_{\text{bp}}(X))$; three eigenvalues are equal
to zero.

Figure 3:  The symmetric Laplacian’s spectrum counts the blocks of a matrix up
to permutations; the spectrum of the unnormalized Laplacian counts both blocks
and zero rows/columns.

###  4.2 Optimization with block diagonal constraints

This section considers the following block diagonally constrained optimization
problem

$\displaystyle\underset{X\in\mathbb{R}^{R\times C}}{\text{minimize}}$
$\displaystyle f(X)$ (12) subject to $\displaystyle X\geq 0\text{ and }X\text{
is block diagonal with at least }B\text{ blocks up to permutations},$

where $f:\mathbb{R}^{R\times C}\to\mathbb{R}$. The naive approach to solving
this problem involves iterating over all possible sparsity patterns with at
least $B$ blocks up to permutations and is likely computationally infeasible.
Based on Proposition 4.1, we see Problem (12) is equivalent to

$\displaystyle\underset{X\in\mathbb{R}^{R\times C}}{\text{minimize}}$
$\displaystyle f(X)$ (13) subject to $\displaystyle X\geq 0\text{ and
}L_{\text{sym}}(A_{\text{bp}}(X))\text{ has at least }B\text{ eigenvalues
equal to }0.$

To impose the rank constraint, we consider the following related problem

$\displaystyle\underset{X\in\mathbb{R}^{R\times C}}{\text{minimize}}$
$\displaystyle
f(X)+\alpha\sum_{j=1}^{B}\lambda_{(j)}\left(L_{\text{sym}}(A_{\text{bp}}(X))\right)$
(14) subject to $\displaystyle X\geq 0,$

for a sufficiently large value of $\alpha$. The non-linearity in
$L_{\text{sym}}(\cdot)$ makes this problem computationally challenging. We can
replace this nonlinearity with linear terms using a variational
characterization of generalized eigenvalues (Proposition B.2 and Corollary
B.1) to obtain,

$\displaystyle\underset{X\in\mathbb{R}^{R\times C},U\in\mathbb{R}^{(R+C)\times
B}}{\text{minimize}}$ $\displaystyle
f(X)+\alpha\text{Tr}\left(U^{T}L_{\text{un}}(A_{\text{bp}}(X))U\right)$ (15)
subject to $\displaystyle X\geq 0$ $\displaystyle
U^{T}\text{diag}(\text{deg}(A_{\text{bp}}(X)))U=I_{B},$

which typically has the same minimizers as (14).

######  Proposition 4.2.

Problems (12) and (13) are equivalent. If $(X,U)$ is a global minimizer of
(15) such that
$\sum_{j=1}^{B}\lambda_{(j)}\left(L_{\text{sym}}(A_{\text{bp}}(X))\right)=0$,
then $X$ is a global minimizer of (12), (13) and (14).

Proposition C.1 gives a similar statement for local solutions.

######  Remark 4.1.

Problem (12) is not guaranteed to have a solution. For example, let
$f(X)=||X-A||_{F}^{2}$ for some matrix $A\in\mathbb{R}^{B\times B}$. If $A=0$
then (12) has no solution (e.g. consider $X_{n}=\frac{1}{n}I_{B}$ then
$f(X_{n})\to 0$ and each $X_{n}$ satisfies the constraints, but $X_{n}\to 0$
which does not satisfy the constraints). On the other hand, if $NB(A)\geq B$
then $X=A$ is the global solution.

Informally, if the objective function $f$ does not encourage too many
rows/columns to be 0, Problem (12) will have a solution. When $f$ does not
have a solution, it may indicate that block diagonal constraints are not a
good modeling choice. For example, it does not make sense to ask for the
nearest block diagonal matrix to the zero matrix.

Problem (15) is amenable to an alternating minimization algorithm that
alternates between updating $U$ and updating $X$. When $X$ is fixed, a global
solution for $U$ is given by an eigen-decomposition. When $U$ is fixed, the
second term in the objective and the second term in the constraints of Problem
(15) are linear in $X$. This alternating algorithm is detailed in Section C.1
and includes the case where $f$ replaced by a surrogate function at each step.
While this algorithm is similar to the BSUM algorithm (Razaviyayn et al.,,
2013; Kumar et al.,, 2019), its convergence properties are more challenging to
study due the non-convexity of $f$ and the non-linearly coupled constraints.
Section C.2 studies the convergence behavior of this alternating algorithm
using Zangwill’s convergence theory.

Section D contrasts our approach based on the symmetric Laplacian with similar
approached based on the unnormalized Laplacian (Nie et al.,, 2016, 2017; Lu et
al.,, 2018; Kumar et al.,, 2019). Section A shows the approach discussed in
this section for matrices naturally extends to enforcing block diagonal
constraints on multi-arrays.

###  4.3 MVMM with block diagonal constraints

This section presents a constrained maximum likelihood problem that imposes a
block diagonal structure on $\pi$ for the MVMM for $V=2$ views. We decompose
$\pi=\epsilon\mathbf{1}\mathbf{1}^{T}+D$ where $\epsilon>0$ is a small
constant and $D$ is a block diagonal matrix. The $\epsilon$ term lets the
model have “outliers” e.g. observations that do not fall cleanly in the block
diagonal structure. It is also useful for computational reasons to avoid
issues with exact zeros similar to those discussed in Section 3. In
particular, we consider

$\displaystyle\underset{\Theta,D}{\text{minimize}}$
$\displaystyle-\ell(\\{x_{i}\\}_{i=1}^{n}|\Theta,\epsilon\mathbf{1}_{K^{(1)}}\mathbf{1}_{K^{(2)}}^{T}+D)$
(16) subject to $\displaystyle D\geq 0,\langle
D,\mathbf{1}_{K^{(1)}}\mathbf{1}_{K^{(2)}}^{T}\rangle=1-K^{(1)}K^{(2)}\epsilon$
$\displaystyle D\text{ has at least }B\text{ blocks up to permutations},$

where $\ell$ is the observed data log-likelihood (3) and
$0<\epsilon<\frac{1}{K^{(1)}K^{(2)}}$. Following Section 4.2, we replace the
block diagonal constraint with a penalty on the smallest generalized
eigenvalues of $L_{\text{sym}}(A_{\text{bp}}(D))$. An alternating EM algorithm
for the resulting problem is presented in Section E.3. Each step of this
algorithm requires an eigen-decomposition and solving a convex problem. Based
on the discussion in Section A, it is straightforward to extend block diagonal
constraints to the case of $V\geq 2$ multi-view mixture models.

##  5 Simulations

We examine the clustering performance of the log penalized MVMM (log-MVMM) and
the block diagonally constrained MVMM (bd-MVMM) on a synthetic data example.
The data in this section are sampled from a $V=2$ view Gaussian mixture model
where $\pi\in\mathbb{R}^{10\times 10}$ has five $2\times 2$ blocks (Figure
5(a) below) with $d^{(1)}=d^{(2)}=10$ features. Each view cluster has an
identity covariance matrix. The cluster means are sampled from isotropic
Gaussians with standard deviations $\sigma_{\text{mean}}^{(1)}$ and
$\sigma_{\text{mean}}^{(2)}$ for the first and second views respectively.
These parameters control the difficulty of the clustering problem e.g. if are
both large then the cluster means tend to be far apart. In this section we set
$\sigma_{\text{mean}}^{(1)}=1$ and $\sigma_{\text{mean}}^{(2)}=.5$ meaning the
clusters in the first view are better separated than those in the second view.
The simulations below are repeated 20 times with different seeds and the
cluster means are sampled once for each Monte-Carlo repetition.

The log-MVMM model is fit for a range of $\lambda$ values and we assume the
number of view clusters $K^{(1)}=K^{(2)}=10$ are known. The bd-MVMM is also
fit for a range of number of blocks and we set
$\epsilon=0.01\cdot\frac{1}{K^{(1)}\cdot K^{(2)}}$. Both the log-MVMM and bd-
MVMM are initialized by fitting the basic MVMM discussed in Section 2 for 10
EM iterations. All models fit in this section assume diagonal covariance
matrices for each cluster and use a small amount of covariance regularization
to prevent clusters from collapsing on a single observation. As baselines for
comparison we also fit the basic MVMM (MVMM) as well as a mixture model on the
concatenated data (cat-MM).

(a) ARI comparing the predicted vs. true cluster labels.

(b) ARI comparing the predicted vs. true block level labels.

Figure 4:  Clustering performance at the true hyper-parameter values; 20
components for log-MVMM, 20 components for cat-MM and 5 blocks for bd-MVMM.
The lines show the Monte-Carlo means; the shaded areas show
$\pm\frac{1}{\sqrt{20}}$ times the Monte-Carlo standard deviation.

We first compare each model when the true parameter values are known e.g.
total number of components for log-MVMM666If the true value does not show up
in the tuning sequence we pick the model with the closest value. and the true
number of blocks for the bd-MVMM. Figure 4 shows the results for a range of
training sample sizes ($n=$ 200, 500, 1000, 1500, 2000, 2500, 3000, 3500,
4000). Recall the adjusted Rand index (ARI) measures how well a vector of
predicted cluster labels corresponds to a vector of true cluster labels where
large values mean better correspondence (Rand,, 1971).

Figure 4(a) shows the ARI of each model’s predicted cluster labels compared to
the true cluster labels for an independent test set (note there are
$|\text{supp}(\pi)|=20$ true clusters). Here the bd-MVMM and log-MVMM perform
better than the two baselines (MVMM and cat-MVMM) for a range of sample sizes.
The prior information about the sparsity structure of $\pi$ helps these two
models estimate the cluster parameters. The performance gap is larger for
smaller sample sizes and narrows with enough data. Note the cat-MM catches up
slowly because it does not take the view structure into account.

Figure 4(b) evaluates the models’ ability to find the block structure of the
$\pi$ matrix. Here we group clusters together that are in the same block i.e.
there are 5 true block clusters. For the log-MVMM and bd-MVMM we predict block
cluster labels based on the block structure of the estimated $\widehat{\pi}$
and $\widehat{D}$ matrices respectively. As a baseline for comparison we apply
bipartite spectral clustering (Dhillon,, 2001) to the estimated
$\widehat{\pi}$ matrix from the MVMM. Here the bd-MVMM performs the best,
which is not surprising because it was designed to target this kind of
structure. The log-MVMM struggles because small mistakes on the
$\widehat{\pi}$ matrix can cause two blocks to be linked. Once the sample size
grows large enough the MVMM eventually comes close to the bd-MVMM.

Figures 5(b) and 5(c) show the estimated $\widehat{D}$ matrix from one Monte-
Carlo repetition. For smaller sample sizes the block diagonal structure is
almost correct (Figure 5(b)). With more samples the bd-MVMM finds the correct
block diagonal structure (Figure 5(c)).

(a) True $\pi\in\mathbb{R}^{10\times 10}$.

(b) Estimated $\widehat{D}$ with $n=2,000$ samples.

(c) Estimated $\widehat{D}$ with $n=2,500$ samples.

Figure 5:  True $\pi$ and estimated $\widehat{D}$ matrices. The cluster labels
have been permuted to reveal the block diagonal structure.

Next we evaluate the models after performing model selection using a modified
BIC criteria (Schwarz et al.,, 1978). After fitting log-MVMM for a range of
$\lambda$ values we select the best model using the following BIC criterion
suggested by (Huang et al.,, 2017)

$\text{BIC}=2\sum_{i=1}^{n}\ell(x_{i}|\widehat{\Theta},\widehat{\pi})-\left(\text{dof}(\widehat{\Theta})+|\text{supp}(\widehat{\pi})|-1\right)\log(n),$
(17)

where $\widehat{\Theta}$ and $\widehat{\pi}$ are the estimated cluster
parameters and $\pi$ matrices respectively and $\text{dof}(\cdot)$ is the
number of degrees of freedom of the cluster parameters. Huang et al., (2017)
provides results about the consistency of this model selection procedure for
single view Gaussian mixture models. For the bd-MVMM we use a similar formula
except the support of $\widehat{\pi}$ is replaced with
$|\text{supp}(\widehat{D})|$.

Figure 6(a) shows the BIC estimated number of components for log-MVMM and bd-
MVMM. The bd-MVMM does a good job with model selection (e.g. it usually picks
5 blocks), but log-MVMM tends to select too few clusters. Figure 6(b) and 6(c)
are similar to Figures 4(a) and 4(b), but the BIC selected parameter values
are used instead of the true values. Here bd-MVMM still outperforms the MVMM,
but by a smaller margin.

(a) Estimated total number of clusters.

(b) ARI of predicted vs. true cluster labels.

(c) ARI of predicted vs. true block level labels.

Figure 6:  Clustering performance of the BIC selected models for log-MVMM and
bd-MVMM.

This section focuses on the case when the signal to noise level is different
in each view. Additional simulations examining different $\pi$ matrices and
different noise levels are shown in Section F. These additional simulations
show that when the noise level is the same in each view
($\sigma_{\text{mean}}^{(1)}=\sigma_{\text{mean}}^{(2)}$) then log-MVMM and
bd-MVMM perform much closer to the MVMM.

##  6 Real data examples

This section applies the block diagonal MVMM to two different data sets. While
more detailed analysis is beyond the scope of this paper we provide additional
results and figures in the online supplementary material.

###  6.1 TCGA breast cancer

The TCGA breast cancer study (Network et al.,, 2012) collects data from 1,027
breast cancer patients on multiple genomic platforms including: RNA expression
(RNA), microRNA (miRNA), DNA methylation (DNA) and copy number (CP). We
closely follow the data processing guidelines from Hoadley et al., (2018),
leaving us with 3,217 RNA features, 383 miRNA features, 3,139 DNA features and
3,000 CP features.777We were unable to obtain the feature list for copy number
so we selected the top 3,000 features with the largest variance. Missing
values are filled in using 5 nearest neighbors imputation (Troyanskaya et
al.,, 2001). We first determine the number of clusters in each view by fitting
a Gaussian mixture model (with diagonal covariances) to each view marginally;
BIC selects 10 RNA clusters, 11 miRNA clusters, 25 DNA clusters and 32 CP
clusters.

Next we fit a $V=2$ view block diagonal MVMM to the following pairings: RNA
vs. miRNA, RNA vs. DNA and RNA vs. CP. BIC selects 1 block for RNA vs. miRNA,
1 block for RNA vs. DNA and 3 blocks for RNA vs. CP. Figure 7(a) shows the
estimated $\widehat{D}$ matrices for RNA vs. CP. The block diagonal structure
of these matrices suggests there is strong jointly defined subtypes in the RNA
and CP views. Note there is still joint information in RNA/miRNA and RNA/DNA
since the estimated $\widehat{D}$ matrices are not rank one.

(a) Estimate $\widehat{D}$ matrix for RNA vs. CP. The rows and columns are
permuted to reveal the 3 blocks.

(b) Contingency table for predicted block labels vs. known PAM50 subtype
labels.

(c) Kaplan-Meier curves comparing the RNA-CP blocks against PFI. A LogRank
test finds the block label are statistically significantly related to PFI.

Figure 7:  The first block picks out Basal like tumors as well as a few
Luminal A tumors. The third block picks out Luminal A tumors that tend to have
better survival based on PFI.

We next investigate the RNA/CP blocks using two additional clinical variables:
PAM50 subtype (Basal-like, Luminal A, Luminal B, Her2-enriched) and survival
measured by progression free interval (PFI) as recommended by (Liu et al.,,
2018). Figure 7(b) shows block 1 picks out the Basal-like subtype, which is
known to have a strong genomic signal in each platform (Network et al.,, 2012;
Hoadley et al.,, 2018). Figures 7(b) and 7(b) show block 3 picks out Luminal A
tumors that have better survival.

###  6.2 Neuron cell types

Integrative clustering has become increasingly important for neuron subtype
discovery (Gouwens et al.,, 2019, 2020). Neuroscientists are now able to
collect a variety of data modalities from individual mouse neurons including
transcriptomic, morphological and electrophysiological features. We apply the
bd-MVMM to an excitatory mouse neuron data set obtained from the Allen
Institute (Gouwens et al.,, 2020). This two-view data set has 44
electrophysiological (EPHYS) features and 69 transcriptomic features (RNA)
available for $n=4,269$ excitatory neurons. The EPHYS features were obtained
using sparse PCA on 12 raw electrophysiological time series recordings in an
awake mouse as discussed in (Gouwens et al.,, 2019). The RNA features are
obtained by first selecting the $4,019$ most differentially expressed genes as
in (Tasic et al.,, 2018), applying a log transform then extracting the top
$69$ PCA features. This PCA rank was selected using the singular value
thresholding method discussed in (Gavish and Donoho,, 2014). We first
determine the number of clusters in each view by fitting a Gaussian mixture
model (with diagonal covariances) to each view marginally; BIC selects 47
EPHYS clusters and 41 RNA clusters.

(a) Estimated $\widehat{D}$ matrix for EPHYS vs. RNA. The rows/columns are
permuted to reveal the block diagonal structure.

(b) Predicted block labels vs. known transcriptomic subtypes. Block 1 picks
out the “Sst Chodl” subtype; block 4 picks out “Lampp5 Lhx6” and “Lamp5 Lsp1”
subtypes.

Figure 8:  MVMM blocks for EPHYS vs. RNA.

(a) Visualization of EPHYS cluster 36. The grey line shows the overall mean
for each time series and the grey shaded area shows $\pm$ 1 standard
deviation. The red lines show the cluster mean of the raw EPHYS recordings as
described in the below caption.

(b) Visualization of the top 50 genes for RNA cluster 33. The values shown are
the standardized difference of the cluster mean minus the overall mean scaled
by the overall sample standard deviation.

Figure 9:  Block 4 identifies EPHYS cluster 36 (Figure 9(a)) with RNA cluster
33 (Figure 9(b)). Both figures show visualizations of the cluster means of the
raw variables. While the clustering algorithm was run on PCA features, we show
the means on the scale of the original features. To represent the mean on the
raw data scale we compute a weighted average of all observations, where the
weights are given by the cluster prediction probabilities (i.e. this is
essentially the M-step for the Gaussian mean parameter).

We next fit a block diagonal MVMM to this two-view data set and select 4
blocks using BIC. There is one large $38\times 43$ block (i.e. 38 RNA clusters
and 43 RNA clusters) and the other blocks are $1\times 1$, $1\times 1$ and
$2\times 1$. This block diagonal structure suggest there are a handful of
jointly well defined clusters while the rest of the information is mixed
between the two views. Previous research has identified 60 RNA clusters (Tasic
et al.,, 2018), which we compare to the predicted block labels found by the
MVMM (Figure 8(b)). This figure shows, for example, block 1 picks out the “Sst
Chodl” subtype and block 4 picks out “Lampp5 Lhx6” and “Lamp5 Lsp1” subtypes.

Figure 9 takes a closer look at block 4 that identifies RNA cluster 33 with
EPHYS cluster 36. Figure 9(a) shows a visualization of EPHYS cluster 36’s mean
for each of the raw EPHYS response variables. This cluster, for example, has a
higher than average “spiking width” and “spiking upstroke downstroke ratio”
responses. Figure 9(b) shows the RNA cluster 33’s mean on the scale of the
standardized residual from the overall mean (i.e. the value shown for each
variable is $\frac{\text{cluster mean}-\text{overall mean}}{\text{sample
standard deviation}}$).

##  7 Conclusion

We presented two methods to estimate the sparsity structure of the $\Pi$
matrix for the multi-view mixture model. The log-MVMM presented in Section 3
makes no assumption about the structure of the sparsity while the bd-MVMM
presented in Section 4 assumes there is a block diagonal structure. These
methods allow scientists to explore how cluster information is spread across
multi-view data sets.

The simulations in Sections 5 and F show the modified BIC criteria often works
well for the block diagonal MVMM, but tends to select too few clusters for the
log penalized MVMM. Future work may establish better model selection methods
e.g. based on Chen and Chen, (2008) or Fu and Perry, (2020).

The main computational bottleneck for the block diagonal MVMM is the convex
Problem (45) in the M-step. For simplicity we use an off the shelf second
order cone program solver (Domahidi et al.,, 2013; Diamond and Boyd,, 2016). A
better algorithm may significantly speed up this step.

## Acknowledgements

We thank Nathan Gowens, Katherine Hoadley, Jonathan Williams, and Daniela
Witten for their insightful discussion and guidance. This material is based
upon work supported by the National Science Foundation under Award No.
1902440.

##  Appendix A Block diagonal multi-arrays

The approach discussed in Section 4 for enforcing block diagonal constraints
on matrices extends to multi-arrays $X\in\mathbb{R}^{d_{1}\times\dots\times
d^{(V)}}$. Consider the following problem

$\displaystyle\underset{X\in\mathbb{R}^{d_{1}\times\dots\times
d^{(V)}}}{\text{minimize}}$ $\displaystyle f(X)$ (18) subject to
$\displaystyle X\geq 0\text{ and }X\text{ is block diagonal with at least
}B\text{ blocks up to permutations}.$

(a) The support of $X$ where the highlighted entries are non-zero. Note there
is a 2d slice of 0s.

(b) Same as Figure 10(a), but the axes have been permuted.

(c) The weight of edge $\\{(1,3),(2,3)\\}$ is
$A_{(1,3),(2,3)}=X_{3,3,2}+X_{3,3,3}$.

Figure 10:  A block diagonal multi-array $X\in\mathbb{R}^{6\times 5\times 5}$
with three blocks up to permutations.

First we extend definitions for matrices given in Section 4 to multi-arrays.
In gory detail, a block of a multi-array is a $V$-hypercube of coordinates,
$\mathcal{B}=[L_{1},U_{1}]\times\dots\times[L_{V},U_{V}]\subseteq\mathbb{Z}^{V}$,
such that $X_{i_{1},\dots,i_{V}}=0$ if there is a $k\in[V]$ such that
$i_{k}\notin[L_{k},U_{k}]$ but $i_{j}\in[L_{j},U_{j}]$ for any $j\neq k$.
Figure 10(a) shows a block diagonal multi-array $X\in\mathbb{R}^{6\times
5\times 5}$ with three blocks e.g. the upper right block is
$[4,5]\times[4,5]\times[4,5]$. For a fixed block diagonal multi-array, we take
the convention that blocks must have at least one non-zero entry and anything
that can be a block is a block. For a multi-array $X$ whose axes are allowed
to be permuted we say “$X$ is block diagonal with $NB(X)$ blocks up to
permutations” where

$\text{NB}(X):=\max\\{B|\text{the axes of }X\text{ can be permuted to create a
}B\text{ block, block diagonal multi-array}\\}$

Any permutation of the axes of $X$ which achieves the above maximum is called
a maximally block diagonal permutation. The multi-arrays in Figures 10(a) and
10(b) both have three blocks up to permutations.

Next we construct a graph that captures the permutation invariant block
diagonal structure of a multi-array. Let $G(X)$ be a weighted, $V$-partite
graph whose vertex set is

$\mathcal{V}:=\\{(v,k)|k\in[d^{(v)}]\text{ and }v\in[V]\\}$

i.e. the generalization of rows and columns to multi-arrays. There can only be
an edge between two vertices on different axes i.e. $(a,j^{(a)})$ and
$(b,j^{(b)})$ where $a\neq b$. The weight of such an edge is given by

$A_{(a,j^{(a)}),(b,j^{(b)})}(X):=\sum_{k^{(1)}=1}^{K^{(1)}}\dots\sum_{k^{(a-1)}=1}^{K^{(a-1)}}\sum_{k^{(a+1)}=1}^{K^{(a+1)}}\dots\sum_{k^{(b-1)}=1}^{K^{(b-1)}}\sum_{k^{(b+1)}=1}^{K^{(b+1)}}\dots\sum_{k^{(V)}=1}^{K^{(V)}}X_{k^{(1)},\dots,j^{(a)},\dots,j^{(b)},\dots
k^{(V)}}$

i.e. summing over all entries of $X$ where the $a$th axis is fixed at
$j^{(a)}$ and the $b$th axis is fixed at $j^{(b)}$. Here
$A(X)\in\mathbb{R}^{\sum_{v=1}^{V}d^{(V)}\times\sum_{v=1}^{V}d^{(V)}}$ is the
adjacency matrix of $G(X)$. This adjacency matrix $A(X)$ is equivalent to the
hypergraph adjacency matrix given in Zhou et al., (2007).

The edge $\\{(a,j_{a}),(b,j_{b})\\}$ is present in $G(X)$ if and only if there
there is a tuple $(k^{(1)},\dots,j^{(a)},\dots,j^{(b)},\dots,k^{(V)})$ such
that $X_{k^{(1)},\dots,j^{(a)},\dots,j^{(b)},\dots,k^{(V)}}\neq 0$. For
example in Figure 10(a) the edge $\\{(2,4),(3,4)\\}$ is present, but the edge
$\\{(2,4),(3,5)\\}$ is not. A vertex $(v,k)$ is isolated if
$X_{j^{(1)},\dots,j^{(V)}}=0$ for all $\\{j^{(1)},\dots,j^{(V)}|j^{(v)}=k\\}$
i.e. there is a $V-1$ dimensional slice of zeros (e.g. $(1,6)$ is the only
isolated vertex in Figure 10(a)).

The symmetric Laplacian of this graph captures the block diagonal structure of
$X$ as follows.

######  Proposition A.1.

The following are equivalent for $1\leq
B+\sum_{v=1}^{V}Z^{(v)}\leq\min(d^{(1)},\dots,d^{(V)})$

  1. 1.

$X$ is block diagonal up to permutations with $B$ blocks and has $Z^{(v)}$
$V-1$ dimensional slices of zeros on the $v$th axis.

  2. 2.

$G(A)$ has $B$ connected components with at least two vertices and
$\sum_{v=1}^{V}Z^{(v)}$ isolated vertices.

  3. 3.

$L_{\text{sym}}(A(X))$ has exactly $B$ eigenvalues equal to 0.

  4. 4.

$L_{\text{un}}(A_{\text{bp}}(X))$ has exactly $B+\sum_{v=1}^{V}Z^{(v)}$
eigenvalues equal to 0.

Additionally, the number of eigenvalues equal to 1 of the symmetric Laplacian
is at least $\sum_{v=1}^{V}Z^{(v)}$.

We now have that Problem (18) is equivalent to

$\displaystyle\underset{X\in\mathbb{R}^{d^{(1)}\times\dots\times
d^{(v)}}}{\text{minimize}}$ $\displaystyle\ell(X)$ (19) subject to
$\displaystyle X\geq 0\text{ and }L_{\text{sym}}(A(X))\text{ has }B\text{
eigenvalues equal to 0}.$

Following Section 4.2, solve the related problem

$\displaystyle\underset{X\in\mathbb{R}^{d^{(1)}\times\dots\times
d^{(V)}},U\in\mathbb{R}^{\sum_{v=1}^{V}d^{(v)}\times B}}{\text{minimize}}$
$\displaystyle f(X)+\alpha\text{Tr}\left(U^{T}L_{\text{un}}(A(X))U\right)$
(20) subject to $\displaystyle X\geq 0\text{ and
}U^{T}\text{diag}(\text{deg}(A(X)))U=I_{B},$

for a sufficiently large value of $\alpha$. Note that $A(\cdot)$ is a linear
function so the second term in the objective and the constraints are linear in
$A$. An alternating algorithm similar to the one discussed in Section C.1 can
be used to solve this problem.

##  Appendix B Extremal characterization of weighted sums of generalized
eigenvalues

For a pair of symmetric matrices $A,B\in\mathbb{R}^{n\times n}$ denote the
matrices whose columns are the largest generalized $K$ eigenvectors of $(A,B)$
by

$\mathcal{GE}_{K}(A,B):=\\{U|AU_{k}=\lambda_{k}BU_{k}\text{ for
}k\in[K],U^{T}BU=I_{K}\\}\subseteq\mathbb{R}^{n\times K}$

where $\lambda_{1}\geq\dots\geq\lambda_{K}$ are the largest generalized
eigenvalues of $(A,B)$. Similarly, let $\mathcal{GE}_{(K)}(A,B)$ be the
analogous set for the smallest $K$ generalized eigenvalues. The following
proposition shows that the generalized eigenvalues of $(A,B)$ can still be
well defined when $B$ has a non-trivial kernel and can be ordered (since they
are real).

######  Proposition B.1.

If $\text{ker}(B)\subseteq\ker(A)$ and $m=n-\text{dim}(\text{ker}(B))$ then
$(A,B)$ has $m$ real generalized eigenvalues. These generalized eigenvalues
are given by the eigenvalues of $B^{-1/2}AB^{-1/2}$ excluding the eigenvalues
whose eigenvectors live in the kernel of $B$ where the inverse is taken to be
the Moore-Penrose pseudo inverse.

We adapt a Proposition from Marshall et al., (1979) to obtain an extremal
characterization for weighted sums of the largest (smallest) generalized
eigenvalues. This is a generalization of the famous Fan’s theorem (Fan,,
1949).

######  Proposition B.2.

Let $A,B\in\mathbb{R}^{n\times n}$ be symmetric. Assume $B$ is positive semi-
definite, $\text{ker}(B)\subseteq\text{ker}(A)$ and $K\leq
n-\text{dim}(\text{ker}(B))$. If $w\in\mathbb{R}^{K}$, such that $w_{1}\geq
w_{2}\geq\dots\geq w_{K}$, then

$\displaystyle\sum_{j=1}^{K}w_{j}\lambda_{j}(A,B)=$
$\displaystyle\underset{U\in\mathbb{R}^{n\times K}}{\text{maximum}}$
$\displaystyle\text{Tr}\left(U^{T}AU\text{diag}(w)\right)$ (21) subject to
$\displaystyle U^{T}BU=I_{K},$

where the maximum is attained by any matrix in $U\in\mathbb{R}^{n\times
K}\in\mathcal{GE}_{K}(A,B)$. Similarly,

$\displaystyle\sum_{j=1}^{K}w_{j}\lambda_{(j)}(A,B)=$
$\displaystyle\underset{U\in\mathbb{R}^{n\times K}}{\text{minimum}}$
$\displaystyle\text{Tr}\left(U^{T}AU\text{diag}(w)\right)$ (22) subject to
$\displaystyle U^{T}BU=I_{K},$

where the minimum is attained by any matrix $U\in\mathbb{R}^{n\times
K}\in\mathcal{GE}_{(K)}(A,B)$.

This proposition allows $U$ to be low rank (as opposed to the full $n\times n$
matrix), permits weighted sums of generalized eigenvalues and says $B$ does
not have to be full rank. Note $w$ is allowed to have negative entries which
may be of interest for some applications (e.g. a penalty that encourages some
eigenvalues to be large).

Proposition B.1 shows how the eigenvalues of
$L_{\text{sym}}(A_{\text{bp}}(X))$ are related to the generalized eigenvalues
of
$(L_{\text{un}}(A_{\text{bp}}(X)),\text{diag}(\text{deg}(A_{\text{bp}}(X))))$;
the latter are the subset of the former excluding the one eigenvalues that
come from degree zero nodes. The following corollary shows how problems (14)
and (15) and are related; they are the same so long as the solution does not
have too many rows/columns of zeros.

######  Corollary B.1.

Let $X\in\mathbb{R}^{R\times C}+$ and $1\leq K\leq\min(R,C)$. Let
$\widetilde{R}$ denote the number of non-zero rows of $X$ (similarly for
$\widetilde{C}$). If $K\leq\widetilde{R}+\widetilde{C}$ then

$\lambda_{(k)}\left(L_{\text{sym}}(A_{\text{bp}}(X))\right)=\lambda_{(k)}\left(L_{\text{un}}(A_{\text{bp}}(X)),\text{diag}(\text{deg}(A_{\text{bp}}(X)))\right)$

for each $1\leq k\leq K$.

##  Appendix C Alternating algorithm for block diagonal constraints

This section considers the following weighted nuclear norm regularized problem
for some $1\leq K\leq\min(R,C)$,

$\displaystyle\underset{X\in\mathbb{R}^{R\times C}}{\text{minimize}}$
$\displaystyle
f(X)+\alpha\sum_{j=1}^{K}w_{j}\lambda_{(j)}\left(L_{\text{sym}}(A_{\text{bp}}(X))\right)$
(23) subject to $\displaystyle X\geq 0,$

where $w\in\mathbb{R}^{K}_{+}$ is a positive weight vector with
$w_{1}\geq\dots\geq w_{K}$ e.g. $w_{k}\propto\frac{1}{k}$. Problem (14) is of
course recovered by setting $w=\mathbf{1}_{K}$.

The mild generalization of (14) allows us to put more weight on smaller
eigenvalues, which can lead to better estimators (Chen et al.,, 2013; Gu et
al.,, 2014). In some applications one might also want to consider (23) where
$\alpha$ is the (continuous) hyper-parameter controlling the amount of block
diagonal regularization instead of (13), which has the (discrete) hyper-
parameter of the number of blocks. We implemented this idea in the context of
the block diagonal MVMM simulations in Section 5. Unfortunately, this
simulation leads to a null result; we found that continuous block diagonal
regularization (e.g. where $w$ is exponentially or polynomially decaying) was
not faster or more accurate than the block diagonally constrained method of
Section 4.3.

The following proposition makes the connection between the constrained (12)
and (14)/(15). While global or local solutions to (12) are difficult to find,
these points are typically contained in a larger set of points (global or
local solutions to (15)), which are easier to find.

######  Proposition C.1.

  1. 1.

Suppose $X$ is a global (local) solution of (14) such that
$\sum_{j=1}^{B}\lambda_{(j)}\left(L_{\text{sym}}(A_{\text{bp}}(X))\right)=0$.
Then $X$ is a global (local) solution of (12).

  2. 2.

Suppose $X$ is a global (local) solution of (14) such that the largest $B$ row
sums and the largest $B$ column sums of $X$ are strictly positive. Then there
exists a coordinate-wise minimizer $U$ such that $(X,U)$ is a global (local)
minimizer of (15).

The first claim shows that if we can find a solution to the penalized Problem
(14) that is block diagonal (i.e. $\alpha$ is large enough to induce the rank
constraint) then we have a solution to (12). The second claim shows the
solutions of (14) are typically888If $\alpha$ is large, the solutions to (14)
typically satisfy the column/row sum condition in the second claim (since zero
rows/columns give large eigenvalues of 1 by Proposition G.3). also solutions
to the extremal representation problem (15); the latter are easier for our
algorithm to locate.

Section C.1 presents an alternating algorithm for (23) and Section C.2
discusses convergence properties of this algorithm.

###  C.1 Alternating algorithm for (23)

Following Proposition 21, we reformulate the weighted nuclear norm problem
(23) as

$\displaystyle\underset{X\in\mathbb{R}^{R\times C},U\in\mathbb{R}^{(R+C)\times
K}}{\text{minimize}}$ $\displaystyle
f(X)+\alpha\text{Tr}\left(U^{T}L_{\text{un}}(A_{\text{bp}}(X))U\text{diag}(w)\right)$
(24) subject to $\displaystyle X\geq
0,U^{T}\text{diag}(\text{deg}(A_{\text{bp}}(X)))U=I_{K}.$

####  C.1.1 $U$ subproblem

For fixed $X$, the $U$ subproblem in (24) is a generalized eigen-problem.
Corollary C.1 shows a global solution of this problem can be obtained through
a low rank SVD of a smaller matrix. For $X\in\mathbb{R}^{R\times C}_{+}$ let

$T_{\text{sym}}(X):=\text{diag}(X\mathbf{1}_{C})^{-1/2}X\text{diag}(X^{T}\mathbf{1}_{R})^{-1/2}\in\mathbb{R}^{R+C}.$
(25)

When $X$ has 0 rows or columns the inverse is taken to be the Moore-Penrose
pseudo-inverse. Note this matrix is the upper diagonal elements of
$L_{\text{sym}}(A_{\text{bp}}(X))$.

The Lasso penalty in the $X$ update discussed below can lead to exact zeros
and, in principle, can introduce some rows/columns of zeros. The following
corollary shows that the $U$ update can handle the case when $X$ has some
rows/columns of 0s. Note that if the initial value of Algorithm 1 satisfies
the condition $K\leq\widetilde{R}+\widetilde{C}$ then the output of each
successive step will also satisfy this condition.

######  Corollary C.1.

For $X\in\mathbb{R}^{R\times C}_{+}$ consider the following problem,

$\displaystyle\underset{U\in\mathbb{R}^{(R+C)\times K}}{\text{minimize}}$
$\displaystyle\text{Tr}\left(U^{T}L_{\text{un}}(A_{\text{bp}}(X))U\text{diag}(w)\right)$
(26) subject to $\displaystyle
U^{T}\text{diag}(\text{deg}(A_{\text{bp}}(X)))U=I_{K}$

for some $K\leq R+C$.

Case 1: Suppose $X$ has no rows or columns of zeros. Let
$U_{\text{left}}\in\mathbb{R}^{R\times\min(R,C)}$ and
$U_{\text{right}}\in\mathbb{R}^{C\times\min(R,C)}$ be the matrix of the left
and right singular vectors of $T_{\text{sym}}(X)$. Let $U_{\text{left},j}$
denote the left singular vector corresponding to the $j$th largest singular
value and let $U_{\text{left},(j)}$ denote the left singular vector
corresponding to the $j$th smallest singular value. If $R\geq C$ let
$Q\in\mathbb{R}^{R\times(\max(R,C)-\min(R,C))}$ be a orthonormal basis matrix
of $\text{col-span}(U_{\text{left}})^{\perp}$. If $R<C$ let
$Q\in\mathbb{R}^{C\times(\max(R,C)-\min(R,C))}$ be a orthonormal basis matrix
of $\text{col-span}(U_{\text{right}})^{\perp}$.

Let the columns of $U^{*}\in\mathbb{R}^{(R+C)\times K}$ be given by
$U^{*}_{k}=\text{diag}(\text{deg}(A_{\text{bp}}(X)))^{-1/2}\Xi$ where

$\Xi=\begin{cases}\begin{bmatrix}U_{\text{left},k}\\\
U_{\text{right},k}\end{bmatrix}&1\leq k\leq\min(R,C)\\\\[1.0pt] \\\
\begin{bmatrix}Q_{j}\\\ \mathbf{0}_{C}\end{bmatrix}&k=\min(R,C)+j,\text{ for
}1\leq j\leq\max(R,C)-\min(R,C),\text{ and }R\geq C\\\\[1.0pt] \\\
\begin{bmatrix}\mathbf{0}_{R}\\\ Q_{j}\end{bmatrix}&k=\min(R,C)+j,\text{ for
}1\leq j\leq\max(R,C)-\min(R,C),\text{ and }R<C\\\\[1.0pt] \\\
\begin{bmatrix}U_{\text{left},(j)}\\\
-U_{\text{right},(j)}\end{bmatrix}&k=\max(R,C)+j,\text{ for }j\geq
1.\end{cases}$ (27)

Then $U^{*}$ is a global minimizer of (26).

Case 2: Suppose $X$ has $\widetilde{R}$ and $\widetilde{C}$ non-zero rows and
columns and $K\leq\widetilde{R}+\widetilde{C}$. Let
$\widetilde{X}\in\mathbb{R}^{\widetilde{R}\times\widetilde{C}}$ denote $X$
after removing the zero rows and columns and let $\widetilde{U}$ be the
solution obtained using (27) applied to $\widetilde{X}$. Then a global
solution of (26) can be obtained by adding appropriate zero rows to
$\widetilde{U}$.

####  C.1.2 $X$ subproblem

For fixed $U$, the constraints and second term in the objective of Problem
(24) are linear in $X$. Let matrix $M(U,w)\in\mathbb{R}^{R\times C}$ be the
matrix such that

$\text{Tr}\left(U^{T}L_{\text{un}}(A_{\text{bp}}(X))U\text{diag}(w)\right)=\langle
X,M(U,w)\rangle.$

Writing $U=\begin{bmatrix}U_{\text{rows}}\\\ U_{\text{col}}\end{bmatrix}$
where $U_{\text{rows}}\in\mathbb{R}^{R\times K}$ and
$U_{\text{cols}}\in\mathbb{R}^{C\times K}$, we see the $r,c$th element of
$M(U,w)$ is

$[M(U,w)]_{rc}=||\text{diag}(w)^{1/2}\left(U_{\text{rows}}(r,:)-U_{\text{cols}}(c,:)\right)||_{2}^{2}.$
(28)

Let

$c_{\text{diag}}(U):=U\odot U\in\mathbb{R}^{(R+C)\times K}$ (29)

be the matrix whose elements are the squares of $U$; this matrix gives the
diagonal elements of the linear equality constraint of (24). Also let

$c_{\text{utri}}(U)\in\mathbb{R}^{(R+C)\times{K\choose 2}}\text{ be the matrix
whose columns are given by }U_{\ell}\odot U_{j},1\leq\ell<j\leq K.$ (30)

This matrix gives the upper-triangular of the linear equality constraints of
(24); the lower triangular constraints are redundant. Note some of the
constraints of $c_{\text{utri}}(U)$ may be redundant999E.g. when
$A_{\text{bp}}(X)$ has multiple connected components Proposition G.1 gives one
source of redundancy. and can be removed to improve numerical performance.

For fixed $U$, the $X$ subproblem for Problem (24) is given by

$\displaystyle\underset{X}{\text{minimize}}$ $\displaystyle f(X)+\alpha\langle
X,M(U,w)\rangle$ (31) subject to $\displaystyle X\geq 0$
$\displaystyle\begin{bmatrix}c_{\text{diag}}(U)^{T}\\\
c_{\text{utri}}(U)^{T}\end{bmatrix}\text{diag}(\text{deg}(A_{\text{bp}}(X)))=\begin{bmatrix}\mathbf{1}_{K}\\\
\mathbf{0}_{{K\choose 2}}\end{bmatrix}.$

If $f$ is convex then (31) is a convex problem because the second term in the
objective and the constraints are linear. Because $X$ is constrained to be
positive, the second term in the objective puts a weighted lasso penalty on
the entries of $X$ whose weights are given by $M(U,w)$.

For complicated objective functions (e.g. the log-likelihood of a mixture
model) the full $X$ updates may be computationally intractable. We therefore
consider surrogate updates obtained by replacing $f$ with a surrogate function
that has the same first order behavior (Razaviyayn et al.,, 2013).

######  Definition C.1.

A surrogate function $Q(X|Y)$ satisfies $Q(X|X)=f(X)$, $Q(X|Y)\geq f(Y)$,
$Q(\cdot|Y)$ is continuous and assume
$\frac{d}{dX}Q(X|Y)\big{|}_{X=Y}=\frac{d}{dX}f(X)\big{|}_{X=Y}$ for all $X,Y$.

Given the current guess, $X_{\text{current}}$, we update $X$ by solving the
following problem,

$\displaystyle\underset{X}{\text{argmin}}$ $\displaystyle
Q(X|X_{\text{current}})+\alpha\langle X,M(U,w)\rangle$ (32) subject to
$\displaystyle X\geq 0$ $\displaystyle\begin{bmatrix}c_{\text{diag}}(U)^{T}\\\
c_{\text{utri}}(U)^{T}\end{bmatrix}\text{diag}(\text{deg}(A_{\text{bp}}(X)))=\begin{bmatrix}\mathbf{1}_{K}\\\
\mathbf{0}_{{K\choose 2}}\end{bmatrix}.$

####  C.1.3 Alternating algorithm

Let $\textsc{Update-X}(X_{\text{current}},U)$ be an algorithm that solves
either the full update (31) or the surrogate update (32).

1

Input: $\alpha\geq 0$, $K\leq\min(R,C)$, $w\in\mathbb{R}^{K}_{+}$

Output: $X$

2

3Initialize $X^{0}$.

4while  _Stopping criteria not satisfied_ do

$U^{s+1}\leftarrow$ smallest $K$ generalized eigenvectors of

// Computed as in Corollary C.1

5
$\left(L_{\text{sym}}(A_{\text{bp}}(X^{s})),\text{diag}(\text{deg}(A_{\text{bp}}(X^{s})))\right)$
(33)

6 $X^{s+1}\leftarrow\textsc{Update-X}(X^{s},U^{s+1})$

7 $s\leftarrow s+1$

Algorithm 1 Alternating algorithm for the weighted nuclear norm Problem (24)

######  Remark C.1.

If Update-X solves either (31) or (32), each step of this algorithm decreases
the objective function of (24) (and (23)). In the former case, Algorithm 1 is
an alternating minimization algorithm while in the latter case it is a block
successive upper bound minimization algorithm with coupled constraints between
blocks (Razaviyayn et al.,, 2013).

####  C.1.4 Algorithm Intuition

The second term in (31) puts a weighted lasso penalty on the entries of $X$.
These weights, which come from (28), encourage $X$ to be more block diagonal.

Suppose $X_{\text{current}}\in\mathbb{R}^{R+C}_{+}$ is exactly block diagonal
up to permutations with $B$ blocks, $K=B$ and $w=\mathbf{1}_{B}$. Let
$\mathbf{1}_{A_{1}},\dots,\mathbf{1}_{A_{B}}\in\mathbb{R}^{R+C}$ denote the
indicator vectors of the blocks and let
$d_{b}=\mathbf{1}_{R+C}^{T}\mathbf{1}_{A_{b}}$ be the total degree of the
$b$th block for each $b\in[B]$. By Proposition G.1,

$U=\begin{bmatrix}\frac{1}{\sqrt{d_{1}}}\mathbf{1}_{A_{1}}&\dots&\frac{1}{\sqrt{d_{B}}}\mathbf{1}_{A_{B}}\end{bmatrix}$

is a $U$ global minimizer of (24). In this case

$[M(U,\mathbf{1}_{B})]_{r,c}=\begin{cases}0&\text{if row }r\text{ and column
}c\text{ are in the same block}\\\
\frac{1}{d_{\text{row}}(r)}+\frac{1}{d_{\text{col}}(c)}&\text{if row }r\text{
and column }c\text{ are in different blocks},\end{cases}$

where $d_{\text{row}}(r)=d_{b}$ where the $r$th row belongs to the $b$th block
(similarly for $d_{\text{col}}(r)$). The second term in (31) only penalizes
edges that go between blocks and does not penalize edges within a block.

If $X_{\text{current}}$ has a row or column of zeros the corresponding
eigenvalue of the symmetric Laplacian will be 1 (i.e. large) and the
corresponding eigenvector will not be included in the smallest $K$
eigenvectors that comprise $U$. Therefore, the algorithm does not want to
encourage rows/columns of zeros.

###  C.2 Convergence of Algorithm 1

We show Algorithm 1 converges to a coordinate-wise minimizer when Update-X
does a full update by solving (31). If Update-X does a surrogate update
solving (32) Algorithm 1 converges to a coordinate-wise stationary point
(defined below). The non-linear coupled constraints of Problem (32) make the
convergence analysis tricky e.g. the BSUM framework does not apply (Razaviyayn
et al.,, 2013).

Consider a constrained optimization problem with two blocks of variables; let
$f(x,y)$ be the objective function and let $g(x,y),h(x,y)$ be vector valued
functions corresponding to the equality and inequality constraints (all
functions are assumed to be continuous). Let $I(x,y)$ denote the indicator
function of the constraint set $\\{(x,y)|g(x,y)=0,h(x,y)\leq 0\\}$. Recall a
stationary point of an optimization problem is one that satisfies the KKT
conditions (Boyd et al.,, 2004); assuming appropriate constraint qualification
all local minimizers are stationary points.

######  Definition C.2.

Let $\mathcal{L}(y):=\\{x*|x^{*}\text{ is a local minimizer of
}\underset{x}{minimize}\;f(x,y)+I(x,y)\\}$ be the set of $x$ coordinate local
minimizers for fixed $y$. Let $\mathcal{S}(y):=\\{x^{*}|x^{*}$ is a stationary
point of $\underset{x}{minimize}\;f(x,y)+I(x,y)\\}$ be the set of $x$
coordinate stationary points for fixed $y$. Let
$\mathcal{G}(x):=\\{y*|y^{*}\in\underset{y}{argmin}\;f(x,y)+I(x,y)\\}$ be the
set of $y$ coordinate global minimizers for fixed $x$. Finally let,

$\mathcal{LG}:=\\{(x,y)|x\in\mathcal{L}(y),y\in\mathcal{G}(x)\\}$ (34)
$\mathcal{SG}:=\\{(x,y)|x\in\mathcal{S}(y),y\in\mathcal{G}(x)\\}$ (35)

denote the set of $x,y$ pairs where $x$ is a coordinate-wise local minimizer
(stationary point) and $y$ is a coordinate-wise minimizer.

######  Assumption C.1.

Assume the objective function $f:\mathbb{R}^{R\times C}_{+}\to\mathbb{R}$ is
continuous and the level set $S_{X^{0}}:=\\{X|f(X)\leq f(X^{0})+\alpha
w^{T}\mathbf{1}_{K}\\}$ is compact where $X^{0}$ is the point at which the
algorithm is initialized.

######  Assumption C.2.

Assume there exists an $\eta>0$ such that the iterates, $X^{s}$, of Algorithm
1 are contained in the set
$\mathcal{RC}_{\eta}:=\\{X|\text{deg}(A_{\text{bp}}(X))\geq\eta\mathbf{1}_{R+C}\\}$
for large enough $s$.

This technical assumption typically hold in practice since the algorithm does
not encourage rows/columns of zeros as discussed above. Alternatively, this
assumption can be enforced by adding the linear constraints
$\text{deg}(A_{\text{bp}}(X))\geq\eta\mathbf{1}_{R+C}$ to (23) and (24). The
updates for the algorithm still work even when some rows/columns of $X$ are
identically zero as long as there are at least $K$ total non-zero rows/columns
at each step101010We lose the convergence guarantees, however, because the
sequence is not guaranteed to be in a compact set..

######  Assumption C.3.

Assume one of the following,

  1. 1.

Update-X returns a global minimizer of the full update problem (31).

  2. 2.

There exists a surrogate function $Q$ and Update-X returns a global minimizer
of the surrogate update problem (32).

######  Proposition C.2.

Let $\\{X^{s},U^{s}\\}_{s=1}^{\infty}$ be any sequence generated by Algorithm
1 and suppose Assumptions C.1 and C.2 hold. Under Assumption C.3.1 all limit
points of $\\{X^{s},U^{s}\\}_{s=1}^{\infty}$ are elements of $\mathcal{LG}$.
Under Assumption C.3.2, all limit points of $\\{X^{s},U^{s}\\}_{s=1}^{\infty}$
are elements of $\mathcal{SG}$. In addition,
$\lim_{s\to\infty}f(X^{s})+\alpha\sum_{j=1}^{K}w_{j}\lambda_{(j)}\left(L_{\text{sym}}(A_{\text{bp}}(X^{s}))\right)\to
f(X^{*})+\alpha\sum_{j=1}^{K}w_{j}\lambda_{(j)}\left(L_{\text{sym}}(A_{\text{bp}}(X^{*}))\right)$
for all limit points $X^{*}$.

Proposition C.2 does not guarantee we find a minimizer of (23). The following
proposition shows these local minimizers are contained in the solution set we
actually are guaranteed to find.

######  Proposition C.3.

Suppose $X$ is a local minimizer of (23). Then $\text{ there exists a }U\text{
such that }(X,U)\in\mathcal{LG}\subseteq\mathcal{SG}$, where $\mathcal{LG}$
and $\mathcal{SG}$ correspond to Problem (24). If $X$ is a global minimizer of
(23), then there exists a $U$ such that $(X,U)$ is a global minimizer of (24).

##  Appendix D Choice of Laplacian

Many existing approaches to imposing block diagonal constraints use the
unnormalized Laplacian instead of the symmetric Laplacian (Feng et al.,, 2014;
Nie et al.,, 2016, 2017; Lu et al.,, 2018; Kumar et al.,, 2019). This section
shows that approaches based on the unnormalized Laplacian require stronger
modeling assumptions and do not have computational advantages over our
approach based on the symmetric Laplacian.

(a) $X_{\epsilon}\in\mathbb{R}^{4\times 4}$ is the matrix with two $2\times 2$
blocks of ones and whose off-diagonal elements are equal to $\epsilon$.

(b) Spectrum of $L_{\text{sym}}(A_{\text{bp}}(X_{\epsilon}))$ for a range of
values of $\epsilon\in[0,1]$.

(c) Spectrum of $L_{\text{un}}(A_{\text{bp}}(X_{\epsilon}))$ for a range of
values of $\epsilon\in[0,1]$.

Figure 11:  As $X_{\epsilon}$ approaches a 2 block, block diagonal matrix an
eigenvalue of both the symmetric and unnormalized Laplacian approaches 0.

(a) $X_{\epsilon}\in\mathbb{R}^{4\times 4}$ is the matrix whose first row is
equal to $\epsilon$ and whose remaining elements are equal to 1.

(b) Same as Figure 11(b). Here the spectrum is the same for every value of
$\epsilon$.

(c) Same as Figure 11(c).

Figure 12:  As a row approaches 0, an eigenvalue of the unnormalized Laplacian
approaches 0; the spectrum of the symmetric Laplacian is unaffected.

For any $X\in\mathbb{R}^{R\times C}_{+}$ Proposition 4.1 shows

$L_{\text{sym}}(A_{\text{bp}}(X))\text{ has exactly }B\text{ eigenvalues equal
to }0\iff X\text{ has exactly }B\text{ blocks up to permutations}$

while

$L_{\text{un}}(A_{\text{bp}}(X))\text{ has exactly }B\text{ eigenvalues equal
to }0\iff X\text{ has at most }B\text{ blocks up to permutations}.$

######  Remark D.1.

Consider replacing $L_{\text{sym}}(\cdot)$ with $L_{\text{un}}(\cdot)$ in
(13). We observed that in practice, using the unnormalized Laplacian for the
block diagonal MVMM often leads to unsatisfactory solutions with too many
rows/columns of 0s.

Figures 11 and 12 illustrate the difference between
$L_{\text{sym}}(A_{\text{bp}}(X))$ and $L_{\text{un}}(A_{\text{bp}}(X))$. When
the symmetric Laplacian has small eigenvalues, then $X$ is close to block
diagonal. When the unnormalized Laplacian has small eigenvalues, $X$ is either
close to block diagonal or has rows/columns of zeros.

It is easier to enforce the constraint “at least $B$ eigenvalues are 0” as
opposed to exactly $B$ eigenvalues are 0. For the symmetric Laplacian this
inequality constraint leads to

$L_{\text{sym}}(A_{\text{bp}}(X))\text{ has at least }B\text{ eigenvalues
equal to }0\iff X\text{ has at least }B\text{ blocks up to permutations}.$

If the inequality constraint is placed on the eigenvalues of the unnormalized
Laplacian we cannot make a corresponding statement about the block diagonal
structure of the matrix.

One approach to ensuring the exact correspondence between the 0 eigenvalues of
the unnormalized Laplacian and the block diagonal structure of $X$ is to
constrain the degrees to be a known, non-zero constant. Let
$c\in\mathbb{R}^{R\times C}_{+}$, with $c>0$ then

$\displaystyle L_{\text{un}}(A_{\text{bp}}(X))\text{ has exactly }B\text{
eigenvalues equal to }0\text{ and }\text{deg}(A_{\text{bp}}(X))=c$ (36)
$\displaystyle\iff$ $\displaystyle X\text{ has at exactly }B\text{ blocks up
to permutations and }\text{deg}(A_{\text{bp}}(X))=c.$

Assuming the degrees are known allows one to use the unnormalized Laplacian,
but requires stronger modeling assumptions.

Using the unnormalized Laplacian with the fixed degree constraint does not
provide computational advantages over our approach based on the symmetric
Laplacian. Each step of the alternating algorithm for the symmetric Laplacian
discussed in Section C.1 computes an eigen-decomposition then solves the
linearly perturbed subproblem (31). A similar algorithm for the unnormalized
Laplacian can also be developed (Nie et al.,, 2016). The eigen-decomposition
for the unnormalized Laplacian requires computing the smallest $K$
eigenvectors of an $\mathbb{R}^{(R+C)\times(R+C)}$ matrix. On the other hand,
the eigen-decomposition for the symmetric Laplacian can be obtained by
computing the largest $K$ singular vectors of a smaller $\mathbb{R}^{R+C}$
matrix (Corollary C.1). Additionally, when the fixed degree constraint is
applied for the unnormalized Laplacian the corresponding linearly perturbed
subproblem is in the same form as (31) (i.e. has linear constraints).

Note that $L_{\text{sym}}(A_{\text{bp}}(\cdot))$ is not a continuous function
near degree zero nodes due to the inverse so we have to be careful about how
we use it. In practice, we find this discontinuity is not a major issue and is
not even present in the extremal formulation of the problem (15). Minimizing
the eigenvalues of the symmetric Laplacian tends not to encourage rows or
columns to be zero, unlike the unnormalized Laplacian (see Figure 12).

##  Appendix E EM algorithms for the multi-view mixture model

This section provides EM algorithms to fit the various multi-view mixture
model problems described in the body of the paper. Many of the computations
(e.g. the E-step and the M-step for the cluster parameters $\Theta$) can be
done using standard single-view mixture model algorithms. This means we can
base implementations of the MVMM EM algorithms off of pre-existing mixture
modeling software such as sklearn (Pedregosa et al.,, 2011).

###  E.1 EM algorithm for the MVMM

We fit the MVMM described in Section 2 by minimizing the negative observed
data log likelihood

$\displaystyle\underset{\Theta,\pi}{\text{minimize}}$
$\displaystyle-\ell(\\{x_{i}\\}_{i=1}^{n}|\Theta,\pi)$ (37)

using an EM algorithm. The E-step constructs a surrogate function for the
original objective function at the current guess and the M-step minimizes this
surrogate function (Lange et al.,, 2000).

In detail, at the $s$th step, given current parameter estimates
$(\Theta^{s},\pi^{s})$ the E-step constructs

$\displaystyle Q^{s}\left(\Theta,\pi\right)$
$\displaystyle:=E\left[\sum_{i=1}^{n}\log
f\left(x_{i},y_{i}|\Theta^{s},\pi^{s}\right)\right]$ (38)
$\displaystyle=\sum_{i=1}^{n}\sum_{k^{(1)}=1}^{K^{(1)}}\dots\sum_{k^{(V)}=1}^{K^{(V)}}\gamma(k^{(1)},\dots,k^{(V)}|x_{i})\log\left(\pi_{k^{(1)},\dots,k^{(V)}}\prod_{v=1}^{V}\phi^{(v)}(x^{(v)}_{i}|\Theta^{(v)}_{k^{(v)}})\right),$

where $f$ is the complete data pdf (1) and the responsibilities are

$\displaystyle\gamma^{s}(k^{(1)},\dots,k^{(V)}|x_{i})$
$\displaystyle:=E\left[P(y=(k^{(1)},\dots,k^{(V)})|x_{i})\right]=\frac{\pi^{s}_{k^{(1)},\dots,k^{(V)}}\prod_{v=1}^{V}\phi^{(v)}(x_{i}^{(v)}|\Theta^{(v),s}_{k^{(v)}})}{f(x_{i}|\Theta^{s},\pi^{s})}\text{
for each }i\in[n].$ (39)

The parameters are then updated in the M-step by solving

$\Theta^{s+1},\pi^{s+1}=\underset{\Theta,\pi}{\text{argmax
}}Q^{s}\left(\Theta,\pi\right)$

From (38) we see that this optimization problem splits into $V+1$ separate
problems; one for $\pi$ and one for each set of view cluster parameters
$\Theta^{(v)}$, $v=1,\dots,V$. The $\pi$ update has an analytical solution
given by $\pi^{s+1}=a$ where $a\in R^{K^{(1)}\times\dots\times K^{(V)}}$ with

$a_{k^{(1)},\dots,k^{(V)}}:=\frac{1}{n}\sum_{i=1}^{n}\gamma^{s}(k^{(1)},\dots,k^{(V)}|x_{i}).$
(40)

The cluster parameters for the $v$th view are updated by solving the following
weighted maximum likelihood problem,

$\displaystyle\underset{\\{\Theta^{(v)}_{k}\\}_{k=1}^{K^{(v)}}}{\text{minimize}}$
$\displaystyle-\sum_{i=1}^{n}\sum_{k=1}^{K^{(v)}}\gamma(k^{(1)},\dots,k^{(V)}|x_{i})\log\left(\phi^{(v)}(x^{(v)}_{i}|\Theta^{(v)}_{k})\right),$
(41)

Note this problem is in exactly the same form as the M-step for a standard,
single-view mixture model making it straightforward to use pre-existing EM
implementations.

1

Input: $K^{(1)},\dots,K^{(V)}$

Data: $\\{x_{i}\\}_{i=1}^{n}$

Output: $\Theta,\pi$

2 Initialize $\Theta^{0}=\\{\Theta^{(v),0}\\}_{v=1}^{V},\pi^{0}$.

3while  _Stopping criteria not satisfied_ do

4

$Q^{s}(\cdot),a^{s}\leftarrow$
E-step$\left(\\{x_{i}\\}_{i=1}^{n},\Theta^{s},\pi^{s}\right)$

// From (38) and (39)

5

6 for  _for v=1, …, V_ do

$\Theta^{(v),s+1}\leftarrow\text{argmin}_{\Theta^{(v)}}Q^{s}(\Theta^{(v)})$

// Solves a problem in the form of (41)

7

8

9 $\pi^{s+1}\leftarrow a^{s}$

10 $s\leftarrow s+1$

Algorithm 2 EM algorithm for the MVMM

The view specific cluster parameters, $\Theta$, can be initialized using
standard mixture model initialization strategies. We initialize the $\pi$
matrix so that the entries all have the same value. We terminate the algorithm
when the objective function has stopped decreasing.

###  E.2 EM algorithm for the log penalized MVMM

This section discusses an EM algorithm for the log-penalized problem (8). This
EM algorithm is similar to the one described in Section E.1, but the M-step
solves a different problem. At each step we majorize the log-likelihood with
$Q^{s}(\Theta,\pi)$ from (38). The updates for $\Theta$ are the same as in
Section E.1. The update of $\pi$ leads to the following problem

$\displaystyle\underset{\pi\in\mathbb{R}^{K^{(1)}\times\dots\times
K^{(V)}}}{\text{minimize}}$
$\displaystyle-\sum_{k^{(1)}=1}^{K^{(1)}}\dots\sum_{k^{(V)}=1}^{K^{(V)}}a_{k^{(1)},\dots,k^{(V)}}\log(\pi_{k^{(1)},\dots,k^{(V)}})+\lambda\log(\delta+\pi_{k^{(1)},\dots,k^{(V)}})$
(42) subject to $\displaystyle\pi\geq 0\text{ and
}\sum_{k^{(1)}=1}^{K^{(1)}}\dots\sum_{k^{(V)}=1}^{K^{(V)}}\pi_{k^{(1)},\dots,k^{(V)}}=1,$

where $a$ is given by (40). Based on Theorem 3.1 we approximate the solution
to this problem with the normalized soft-thresholding operation

$\pi_{k^{(1)},\dots,k^{(V)}}=\frac{(a_{k^{(1)},\dots,k^{(V)}}-\lambda)_{+}}{\sum_{j^{(1)}=1}^{K^{(1)}}\dots\sum_{j^{(V)}=1}^{K^{(V)}}(a_{j^{(1)},\dots,j^{(V)}}-\lambda)_{+}}.$
(43)

1

Input: $K^{(1)},\dots,K^{(V)}$, $0<\lambda<\frac{1}{\prod_{v=1}^{V}K^{(v)}}$

Data: $\\{x_{i}\\}_{i=1}^{n}$

Output: $\Theta,\pi$

2 Initialize $\Theta^{0}=\\{\Theta^{(v),0}\\}_{v=1}^{V},\pi^{0}$.

3while  _Stopping criteria not satisfied_ do

4

$Q^{s}(\cdot),a^{s}\leftarrow$
E-step$\left(\\{x_{i}\\}_{i=1}^{n},\Theta^{s},\pi^{s}\right)$

// From (38) and (39)

5

6 for  _for v=1, …, V_ do

$\Theta^{(v),s+1}\leftarrow\text{argmin}_{\Theta^{(v)}}Q^{s}(\Theta^{(v)})$

// Solves a problem in the form of (41)

7

8

$\pi^{s+1}\leftarrow$ normalized soft-thresholding applied to $a^{s}$ as in
(43)

// Approximates Problem (42)

9

10 $s\leftarrow s+1$

Algorithm 3 EM algorithm for the log-penalized MVMM, Problem (8)

Algorithm 4 is initialized by running a few EM steps for the unconstrained
MVMM using the algorithm discussed in Section E.1. We terminate the algorithm
when the objective function of (8) has stopped decreasing. We specify a small
value of $\delta$ (e.g. $10^{-6}$) to monitor the convergence of the objective
function, but this value of $\delta$ plays no role in the EM updates due to
(43).

###  E.3 EM algorithm for the block diagonally constrained MVMM

Following Sections 4.2 and C we replace (16) with the following related
problem for a sufficiently large value of $\alpha$,

$\displaystyle\underset{\Theta,D,U}{\text{minimize}}$
$\displaystyle-\ell(\\{x_{i}\\}_{i=1}^{n}|\Theta,\epsilon\mathbf{1}_{K^{(1)}}\mathbf{1}_{K^{(2)}}^{T}+D)+\alpha\text{Tr}\left(U^{T}L_{\text{un}}(A_{\text{bp}}(D))U\right)$
(44) subject to $\displaystyle D\geq 0,\langle
D,\mathbf{1}_{K^{(1)}}\mathbf{1}_{K^{(2)}}^{T}\rangle=1-K^{(1)}K^{(2)}\epsilon$
$\displaystyle U^{T}\text{diag}(\text{deg}(A_{\text{bp}}(D)))U=I_{B}.$

We can solve this problem by alternating between updating $U$ and
$(\Theta,\pi)$. The $U$ variable is updated with an eigen-decomposition as in
Corollary C.1. To update $(\Theta,D)$ at the $s$th step we majorize the log-
likelihood with
$Q^{s}(\Theta,\epsilon\mathbf{1}_{K^{(1)}}\mathbf{1}_{K^{(2)}}^{T}+D)$ from
(38). The update for $\Theta$ is the same as in Section E.1. The M-step for
$D$ solves the following convex problem

$\displaystyle\underset{D}{\text{minimize}}$
$\displaystyle-\sum_{k^{(1)}=1}^{K^{(1)}}\sum_{k^{(2)}=1}^{K^{(2)}}a_{k^{(1)}k^{(2)}}\log(\epsilon+D_{k^{(1)}k^{(2)}})+\alpha\langle
D,M(V,\mathbf{1}_{B})\rangle$ (45) subject to $\displaystyle D\geq 0,\langle
D,\mathbf{1}_{K^{(1)}}\mathbf{1}_{K^{(2)}}^{T}\rangle=1-K^{(1)}K^{(2)}\epsilon$
$\displaystyle\begin{bmatrix}c_{\text{diag}}(V)^{T}\\\
c_{\text{utri}}(V)^{T}\end{bmatrix}\text{diag}(\text{deg}(A_{\text{bp}}(D)))=\begin{bmatrix}\mathbf{1}_{B}\\\
\mathbf{0}_{{B\choose 2}}\end{bmatrix},$

where $a$ is from (40) and
$M(U,\mathbf{1}_{B}),c_{\text{diag}}(U),c_{\text{utri}}(U)$ are from (28),
(29), (30). Let Update-D(U) be an algorithm that solves the convex Problem
(45).

1

Input: $K^{(1)},K^{(2)}$, $1\leq B\leq\min(K^{(1)},K^{(2)})$,
$0<\epsilon<\frac{1}{K^{(1)}K^{(2)}}$, $\alpha>0$

2

Data: $\\{x_{i}\\}_{i=1}^{n}$

Output: $\Theta,D$

3

4Initialize $\Theta^{0}=\\{\Theta^{(v),0}\\}_{v=1}^{2},D^{0}$.

Initialize $\alpha$

// E.g. from (47) below

5

6while  _Block diagonal stopping criteria is not satisfied_ do

7

8 while  _Optimization convergence stopping criteria not satisfied_ do

9

$U^{s+1}\leftarrow$ smallest $B$ generalized eigenvectors of

// Computed as in Corollary C.1

10
$\left(L_{\text{un}}(A_{\text{bp}}(D^{s})),\text{diag}(\text{deg}(A_{\text{bp}}(D^{s})))\right)$

$Q^{s}(\cdot)\leftarrow$
E-step$\left(\\{x_{i}\\}_{i=1}^{n},\Theta^{s},\epsilon\mathbf{1}_{K^{(1)}}\mathbf{1}_{K^{(2)}}^{T}+D^{s}\right)$

// From (38) and (39)

11

12 for  _for v=1, 2_ do

$\Theta^{(v),s+1}\leftarrow\text{argmin}_{\Theta^{(v)}}Q^{s}(\Theta^{(v)})$

// Solves a problem in the form of (41)

13

14

15 $D^{s+1}\leftarrow\textsc{Update-D}(U^{s+1})$

16 $s\leftarrow s+1$

17

Increase $\alpha$

// E.g. $\alpha\leftarrow 2*\alpha$

18

Algorithm 4 EM algorithm for the block diagonally constrained MVMM, Problem
(44)

Algorithm 4 is initialized by running a few EM steps for the unconstrained
MVMM using the algorithm discussed in Section E.1. Each step of the inner loop
of Algorithm 4 is guaranteed to decrease the objective function of (44)
therefore, we stop the inner loop when the objective function has stopped
decreasing. The convergence results discussed in Section C apply to the inner
loop of Algorithm 4.

For a given value of $\alpha$, the solution output by Algorithm 4 may have too
few 0 eigenvalues; in this case we increase $\alpha$ (e.g. multiplying it by
2) and re-run the inner loop. The following proposition motivates a heuristic
choice for the initial value of $\alpha$ as well as an initializer for
Update-D.

######  Proposition E.1.

Let $a,b,\epsilon>0$. The unique global minimizer,

$\displaystyle x^{*}=$ $\displaystyle\underset{x\in\mathbb{R}}{\text{argmin}}$
$\displaystyle-a\log(x+\epsilon)+bx$ (46) subject to $\displaystyle x\geq 0$

is given by

$x^{*}=\begin{cases}\frac{a}{b}-\epsilon&\text{ if }\frac{a}{b}-\epsilon>0\\\
0&\text{otherwise.}\end{cases}$

Let $(\Theta^{0},D^{0})$ be the initial guess in Algorithm 4. By ignoring
constraints the solution to (45) can be approximated by

$D_{k^{(1)},k^{(2)}}^{*}\approx\left(\frac{a_{k^{(1)},k^{(2)}}}{\alpha
M(U,\mathbf{1}_{B})_{k^{(1)},k^{(2)}}}-\epsilon\right)_{+},$

where $a$ and $U$ are obtained from $(\Theta^{0},D^{0})$ as above. This
suggests the following guess111111This median value gives a rough estimate for
the scale of $\alpha$ at which terms are set to 0. for $\alpha$

$\alpha=c\cdot\text{median}\left(\left\\{\frac{a_{k^{(1)},k^{(2)}}}{\epsilon
M(U,\mathbf{1}_{B})_{k^{(1)},k^{(2)}}}\right\\}_{k^{(1)}\in[K^{(1)}],k^{(2)}\in[K^{(2)}]}\right)$
(47)

for some small value of $c<0$ e.g. $c=0.01$.

Algorithm 4 can be sensitive to the initial choice of $\alpha$ and how fast
$\alpha$ is increased. Informally, if $\alpha$ is too large, the algorithm may
converge quickly to a bad local minimizer. If $\alpha$ is too small the
algorithm will take longer to converge.

##  Appendix F Additional simulations

This section expands on the simulations presented in Section 5. The setup here
is similar to the setup in Section 5. Here we look at three different $\Pi$
matrices (Figure 13) and at two different signal to noise settings. In the
first setting the views have uneven signal to noise ratio where
$\sigma_{\text{mean}}^{(1)}=1$ and $\sigma_{\text{mean}}^{(2)}=0.5$ (i.e. the
first view clusters are better separated than the second view clusters). In
the second setting the views have even signal to noise ratios where
$\sigma_{\text{mean}}^{(1)}=\sigma_{\text{mean}}^{(2)}=1$. The figures below
examine cluster level performance at the true parameter values (similar to
Figure 4(a)), block level performance at the true parameter values (similar to
Figure 4(b)) and the BIC estimated number of components (similar to Figure
6(a)). The details of these figures are explained in Section 5.

(a) $\Pi\in\mathbb{R}^{10\times 10}$ with five $2\times 2$ blocks. All entries
are equal.

(b) $\Pi\in\mathbb{R}^{10\times 10}$ with five $1\times 1$ blocks and one
$5\times 5$ block. All blocks have the same total weight.

(c) $\Pi\in\mathbb{R}^{10\times 10}$ with $18$ randomly selected, non-zero
entries with equal values.

Figure 13:  Three different $\Pi$ matrices examined by simulations in this
section.

The overall takeaway is that the bd-MVMM and log-MVMM usually outperform the
MVMM in the uneven setting. In the even setting the bd-MVMM and log-MVMM
sometimes still have an edge over the MVMM, but all three models are much
closer together. The log-MVMM sometimes struggles with the block level
performance because small errors in the support of the estimated $\Pi$ can
merge blocks together. The BIC criteria often, but not always works well for
the bd-MVMM. This BIC criteria is usually biased towards selecting too few
components for the log-MVMM.

Figure 14 shows the results for the $\Pi$ matrix shown in 13(a). The top row
shows the uneven setting and the bottom row shows the even setting. In the
uneven setting the log-MVMM and bd-MVMM out-perform the MVMM in both cluster
level and block level performance. In the even setting the log-MVMM and bd-
MVMM are only slightly better than the MVMM. In the uneven setting the log-
MVMM struggles with BIC based model selection, though it works well in the
even setting.

Figure 15 shows the results for the $\Pi$ matrix shown in 13(b). In the uneven
setting the bd-MVMM performs the best on both the cluster level and block
level labels, however, the log-MVMM struggles with the block labels. In this
uneven setting BIC does not work well for either model and selects too few
clusters in both cases. BIC may struggle with the bd-MVMM because the
individual clusters in the large $5\times 5$ block have smaller weights and
therefore breaking this block up into several blocks does not harm the model
fit as much as with the other $\Pi$. In the even setting all three models
perform similarly though the bd-MVMM has a slight edge at smaller sample
sizes. In this even setting BIC works well for bd-MVMM, but is still biased
down for log-MVMM.

Figure 16 shows the results for the sparse $\Pi$ matrix shown in 13(c). Again
in the uneven setting the log-MVMM out performs the MVMM, but in the even
setting the two models are much closer together. BIC is still biased towards
too few clusters for this setting.

Figure 14:  Results for the 5 block $\Pi$ matrix shown in Figure 13(a). In the
top row the view signal to noise ratios are uneven with
$\sigma_{\text{mean}}^{(1)}=1$ and $\sigma_{\text{mean}}^{(2)}=0.5$. In the
bottom row the view signal to noise ratios are even with
$\sigma_{\text{mean}}^{(1)}=\sigma_{\text{mean}}^{(2)}=1$. The first two
columns examine the cluster label and block label performance at the true
hyper-parameter values. The third column examines the BIC estimated number of
components.

Figure 15:  Results for the 6 block $\Pi$ matrix shown in Figure 13(b). The
top row shows the results for the uneven signal to noise ratio and the bottom
row shows the results for the even signal to noise ratio.

Figure 16:  Results for the sparse $\Pi$ shown in Figure 13(c). The top row
shows the results for the uneven signal to noise ratio and the bottom row
shows the results for the even signal to noise ratio. Only the results for
log-MVMM are shown.

##  Appendix G Proofs

###  G.1 Soft-thresholding with log penalty

Let $f(\pi)=a\log(\pi)-\lambda\log(\delta+\pi)$ (see Figure 17). The intuition
that Problem (6) sets terms to (approximately) 0 comes from the following.

  * •

The solution to the unpenalized Problem (6) when $\lambda=0$ is $\pi^{*}=a$.

  * •

If $a<\lambda$, $f(\pi)$ has a global maximum as
$\pi^{*}=\frac{a\delta}{\lambda-a}\propto\delta$.

As $\delta\to 0$, the terms where $a_{k}<\lambda$ go to zero. These terms
become negligible in the probability constraint so the terms
$\\{k|a_{k}>\lambda\\}$ solve the unpenalized problem (i.e. the problem if
$\lambda=0$) with coefficients $a_{k}-\lambda$.

(a) If $\lambda>a$ there is a global maximizer at $\frac{\delta
a}{\lambda-a}$.

(b) If $\lambda<a$ the function is strictly increasing and concave.

Figure 17:  Graph of $f(\pi)$ for two cases if $\lambda<a$ or $\lambda>a$.
When $\lambda>a$, there is a global maximizer, which is proportional to
$\delta$.

###### Proof.

of Theorem 3.1

First we check there exists at least one global minimizer. There exists a
$\xi>0$ such that $\pi\in[\xi,1-K\xi]^{K^{C}}\implies
f(\pi)<f(\frac{1}{K}\mathbf{1}_{K})$. Therefore, an optimal solution of (6) is
the same as an optimal solution for the restricted problem where
$\pi_{k}\in[\xi,1-K\xi]$ for $k=1,\dots,K$.. This restricted problem is a
continuous function over a compact set and must attain a minimum thus (6) has
at least one global minimizer.

Linear constraint qualification holds so the KKT conditions are first order
necessary. The Lagrangian of (6) is given by

$L(\pi,\eta)=\sum_{k=1}^{K}a_{k}\log(\pi_{k})-\lambda\sum_{k=1}^{K}\log(\delta+\pi_{k})-\eta\pi^{T}\mathbf{1}_{K}$

for $\eta\in\mathbb{R}$. We ignore the positivity constraint because the
$-\log(z_{k})$ terms ensure any stationary point is strictly positive. The
gradient of the Lagrangian is given by

$\frac{d\mathcal{L}}{d\pi_{k}}=\frac{a_{k}}{\pi_{k}}-\frac{\lambda}{\pi_{k}+\delta}-\eta\;\text{
for }k=1,\dots,K.$

Suppose $(\pi,\eta)$ is a stationary point of the Lagrangian. Setting
$\frac{d\mathcal{L}}{d\pi_{k}}=0$ leaves us with

$\eta\pi_{k}^{2}+(\lambda+\eta\delta-a_{k})\pi_{k}-a_{k}\delta=0$ (48)

Because $\sum_{k}a_{k}=1$ and $\lambda<\frac{1}{K}$, without loss of
generality $a_{1}\geq\frac{1}{k}>\lambda$ so $a_{1}>\lambda$. If $\eta=0$ then
by (48), $\pi_{k}=\frac{a_{k}\delta}{\lambda-a_{k}}$ for each $k$. In this
case, $\pi_{1}<0$ violates the positivity constraint so we conclude $\eta\neq
0$. Thus

$\pi_{k}\in\left\\{\frac{(a_{k}-\lambda-\eta\delta)\pm\sqrt{(a_{k}-\lambda-\eta\delta)^{2}+4a_{k}\eta\delta}}{2\eta}\right\\}$
(49)

Next we check that $\eta>0$ at any stationary point of the Lagrangian. Assume
for the sake of contradiction that $\eta<0$. Recall $\lambda<a_{1}$ so
$a_{1}-\lambda-\eta\delta>0$. Thus

$\frac{(a_{k}-\lambda-\eta\delta)+\sqrt{(a_{k}-\lambda-\eta\delta)^{2}+4a_{k}\eta\delta}}{2\eta}<0$

which violates the constraint $\pi_{1}\geq 0$. Furthermore,
$(a_{k}-\lambda-\eta\delta)-\sqrt{(a_{k}-\lambda-\eta\delta)^{2}+4a_{k}\eta\delta}>0$
so

$\frac{(a_{k}-\lambda-\eta\delta)-\sqrt{(a_{k}-\lambda-\eta\delta)^{2}+4a_{k}\eta\delta}}{2\eta}<0$

which again violates the constraint $\pi_{1}\geq 0$. Therefore we conclude
$\eta>0$ for $\pi$ to be a stationary point.

It can be checked that if a $(-)$ is chosen in (49) then $\pi_{1}<0$. Thus at
a stationary point

$\pi_{k}=\frac{(a_{k}-\lambda-\eta\delta)+\sqrt{(a_{k}-\lambda-\eta\delta)^{2}+4a_{k}\eta\delta}}{2\eta}$
(50)

Next we show $\eta\delta\to 0$ when $\delta\to 0$. Using the constraint
$1=\sum_{k=1}^{K}\pi_{k}$ we get

$2\eta=\sum_{k=1}^{K}\left(a_{k}-\lambda-\eta\delta+\sqrt{(a_{k}-\lambda-\eta\delta)^{2}+4a_{k}\eta\delta}\right)$
(51)

It can be checked that

$(a_{k}-\lambda-\eta\delta)^{2}+4a_{k}\eta\delta\leq(a_{k}+\lambda+\eta\delta)^{2},$

therefore

$\displaystyle 2\eta$
$\displaystyle\leq\sum_{k=1}^{K}|a_{k}-\lambda-\eta\delta|+\sqrt{(a_{k}-\lambda-\eta\delta)^{2}+4a_{k}\eta\delta}\leq\sum_{k=1}^{K}|a_{k}-\lambda-\eta\delta|+|a_{k}+\lambda+\eta\delta|$
$\displaystyle\leq 2K\eta\delta+\sum_{k=1}^{K}|a_{k}-\lambda|+a_{k}+\lambda.$

Thus

$\eta\leq\frac{1}{2}\frac{\sum_{k=1}^{K}|a_{k}-\lambda|+|a_{k}+\lambda|}{1-K\delta}$
(52)

which is upper bounded by a positive constant independent of $\delta$ because
$\delta<\frac{1}{K}$. We now conclude

$\lim_{\delta\to 0}\eta\delta=0$ (53)

Finally, using (51) and (53) we see

$\displaystyle\lim_{\delta\to 0}\eta$
$\displaystyle=\frac{1}{2}\lim_{\delta\to
0}\sum_{k=1}^{K}\left(a_{k}-\lambda-\eta\delta+\sqrt{(a_{k}-\lambda-\eta\delta)^{2}+4a_{k}\eta\delta}\right)$
$\displaystyle=\frac{1}{2}\sum_{k=1}^{K}a_{k}-\lambda+|a_{k}-\lambda|=\sum_{k=1}^{K}(a_{k}-\lambda)_{+},$

and the result follows from (50).

∎

###  G.2 Extremal characterization of generalized eigenvalues

###### Proof.

of Proposition B.1

Let $B=VDV^{T}$ be the eigen-decomposition of $B$ where
$V\in\mathbb{R}^{m\times m}$ and $D\in\mathbb{R}^{m\times m}$ is the diagonal
matrix of non-zero eigenvalues of $B$. Let $Z\in\mathbb{R}^{(n-m)\times n}$ be
such that $Q=[V;Z]$ is an orthonormal matrix (i.e. $Z$ is a basis for the
kernel).

Then $B^{-1/2}:=QD^{-1/2}Q^{T}$ is the Moore-Penrose inverse of the square
root of $B$ and

$B^{-1/2}AB^{-1/2}=\begin{bmatrix}\widetilde{A}&0\\\ 0&0\end{bmatrix},$

where $\widetilde{A}\in\mathbb{R}^{m\times m}$ since
$\text{ker}(B)\subseteq\text{ker}(A)$. Let $\lambda_{1},\dots,\lambda_{m}$ be
the eigenvalues of $\widetilde{A}$ with corresponding eigenvectors
$\widetilde{v}_{1},\dots\widetilde{v}_{m}\in\mathbb{R}^{m\times m}.$ Let
$v_{j}\in\mathbb{R}^{n}$ be the concatenation of $\widetilde{v}_{j}$ along
with $n-m$ zeros. Then the $v_{j}$ are eigenvectors of $B^{-1/2}AB^{-1/2}$
with eigenvalues $\lambda_{j}$. Setting $w_{j}=B^{1/2}v_{j}$ we see $w_{j}$
are generalized eigenvectors of $(A,B)$ with generalized eigenvalues
$\lambda_{j}$.

∎

###### Proof.

of Proposition B.2

First assume that $B$ is positive definite. Proposition 20.A.2.a from Marshall
et al., (1979) states

$\displaystyle\sum_{j=1}^{K}\lambda_{i}(A)\lambda_{i}(H)=$
$\displaystyle\underset{U\in\mathbb{R}^{n\times K}}{\text{maximum}}$
$\displaystyle\text{tr}\left(U^{T}AUH\right)$ (54) subject to $\displaystyle
U^{T}U=I_{K}.$

It is straightforward to check that this maximum value is attained by
$U=U_{A}(:,\;1:K)U_{H}$ where $U_{A}(:,\;1:K)\in\mathbb{R}^{n\times K}$ is an
orthonormal matrix of eigenvectors corresponding to the largest $K$
eigenvalues of $A$ and $U_{H}\in\mathbb{R}^{K\times K}$ is an orthonormal
matrix of eigenvectors of $H$.

If $B=I_{K}$ then (21) is a special case of (54) with $H=\text{diag}(w)$.
Recall that $(\lambda,u)$ is a generalized eigenvector of $(A,B)$ if and only
if $(\lambda,B^{1/2}u)$ is an eigenvector of $B^{-1/2}AB^{-1/2}$ and (21)
follows. Using this fact it is straightforward to extend the results for
general, positive definite $B$. Repeating this argument with
$H=-\text{diag}(w)$ we obtain (22).

Next we relax the positive definite assumptions; assume
$\text{ker}(B)\subseteq\text{ker}(A)$. Let $Q\in\mathbb{R}^{n\times n}$ be an
orthonormal matrix whose columns are $q_{1},\dots,q_{n}\in\mathbb{R}^{n}$.
Suppose $q_{1},\dots,q_{m}\in\mathbb{R}^{n}$ is an orthonormal basis of
$\text{ker}(B)^{\perp}$ and $q_{m+1},\dots,q_{n}\in\mathbb{R}^{n}$ is an
orthonormal basis of $\text{ker}(B)\subseteq\text{ker}(A)$ where
$m=n-\text{dim}(\text{ker}(B))\leq K$ by assumption. Then

$Q^{T}AQ=\begin{bmatrix}\widetilde{A}&0\\\ 0&0\end{bmatrix}\text{ and
}Q^{T}BQ=\begin{bmatrix}\widetilde{B}&0\\\ 0&0\end{bmatrix}$ (55)

where $\widetilde{A},\widetilde{B}\in\mathbb{R}^{m\times m}$ and
$\widetilde{B}$ is positive definite. Let $\widetilde{U}\in\mathbb{R}^{m\times
K}$ be generalized eigenvectors corresponding to the largest $K$ generalized
eigenvalues of $(\widetilde{A},\widetilde{B})$ and let

$U=Q^{T}\begin{bmatrix}\widetilde{U}\\\ 0\end{bmatrix}\in\mathbb{R}^{n\times
K}$

It is straightforward to see that the columns of $U$ are generalized
eigenvectors for $(A,B)$. We next check that $U$ is a global maximizer of
(21). First note

$U^{T}BU=\begin{bmatrix}\widetilde{U}\\\
0\end{bmatrix}QBQ^{T}\begin{bmatrix}\widetilde{U}\\\
0\end{bmatrix}=\begin{bmatrix}\widetilde{U}\\\
0\end{bmatrix}^{T}\begin{bmatrix}\widetilde{B}&0\\\
0&0\end{bmatrix}\begin{bmatrix}\widetilde{U}\\\
0\end{bmatrix}=\widetilde{U}^{T}\widetilde{B}\widetilde{U}=I_{K}$

so $U$ is feasible and its objective value is given by

$\displaystyle\text{Tr}\left(U^{T}AU\text{diag}(w)\right)$
$\displaystyle=\text{Tr}\left(\begin{bmatrix}\widetilde{U}\\\
0\end{bmatrix}^{T}Q^{T}AQ\begin{bmatrix}\widetilde{U}\\\
0\end{bmatrix}\text{diag}(w)\right)=\text{Tr}\left(\widetilde{U}^{T}\widetilde{A}\widetilde{U}\text{diag}(w)\right)$

Suppose that $W\in\mathbb{R}^{n\times K}$ is a global maximizer of (21). Let
$\widetilde{W}\in\mathbb{R}^{m\times K},\overline{W}\in\mathbb{R}^{n-m\times
K}$ such that

$Q^{T}W=\begin{bmatrix}\widetilde{W}\\\
\overline{W}\end{bmatrix}\in\mathbb{R}^{n\times K}$

Noting that $Q$ is orthonormal and using (55) we see

$I_{K}=W^{T}BW=W^{T}QQ^{T}BQ^{T}QW=\widetilde{W}^{T}\widetilde{B}\widetilde{W},$

and

$\displaystyle\text{Tr}\left(W^{T}AW\text{diag}(w)\right)$
$\displaystyle=\text{Tr}\left(W^{T}QQ^{T}AQ^{T}QW\text{diag}(w)\right)=\text{Tr}\left(W^{T}Q\begin{bmatrix}\widetilde{A}&0\\\
0&0\end{bmatrix}QW\text{diag}(w)\right)$
$\displaystyle=\text{Tr}\left(\widetilde{W}^{T}\widetilde{A}\widetilde{W}\text{diag}(w)\right).$

Assume for the sake of contradiction that
$\text{Tr}\left(W^{T}AW\text{diag}(w)\right)>\text{Tr}\left(U^{T}AU\text{diag}(w)\right)$.
Consider Problem (21) with $(\widetilde{A},\widetilde{B})$. By the first part
of the proof, $\widetilde{U}$ is a global maximizer (since $\widetilde{B}$ is
strictly positive definite). From the above discussion we have that
$\widetilde{W}$ is feasible for this problem with objective value
$\text{Tr}\left(W^{T}AW\text{diag}(w)\right)$, however, this contradicts the
fact that $\widetilde{U}$ is a global maximizer.

∎

Corollary B.1 is proved below in Section 4.1.

###  G.3 Spectrum of the symmetric Laplacian and block diagonal structure

We first give two propositions detailing the spectral properties of the
symmetric, normalized Laplacian. Recall the convention for isolated vertices
discussed in Section 4.1 that ensures the diagonal of $L_{\text{sym}}(\cdot)$
is always equal to 1.

######  Proposition G.1.

Let $A\in\mathbb{R}^{n\times n}_{+}$ be the adjacency matrix of a graph with
undirected, positively weighted edges and no self loops (i.e. $A$ is symmetric
and has 0s on the diagonal).

There is a one-to-one correspondence between the 0 eigenvalues of
$L_{\text{sym}}(A)$ and the connected components of $G$ with at least two
vertices. Let $A_{1},A_{2},\dots,[n]$ correspond to the connected components
with at least two vertices and let
$v_{i}=\text{diag}(\text{deg}(A))\mathbf{1}_{A_{i}}$ where
$\mathbf{1}_{A_{i}}$ is the vector with 1s in the entries corresponding to
indices in $A_{i}$ and 0s elsewhere. The eigenspace of 0 is spanned by
$v_{1},v_{2},\dots$.

The number of eigenvalues of $L_{\text{sym}}(A)$ equal to 1 is at least the
number of isolated vertices. The basis vector with a 1 in the entry
corresponding to an isolated vertex is an eigenvector with eigenvalue 1.

In general there is not a one-to-one correspondence between isolated vertices
and 1 eigenvalues e.g. consider

$L_{\text{sym}}(A_{\text{bp}}(\mathbf{1}_{m}\mathbf{1}_{m}^{T}))$

which has no isolated vertices but $2m$ eigenvalues equal to 1.

Note the difference between Proposition G.1 and Proposition 4 of Von Luxburg,
(2007). Some papers choose the convention that the symmetric Laplacian has a 0
on the diagonal for isolated vertices. When this alternative convention is
selected, the symmetric, normalized Laplacian and the normalized Laplacian
would treat isolated vertices the same.

######  Proposition G.2.

Let $X\in\mathbb{R}^{R\times C}_{+}$, then the eigenvalues of
$L_{\text{sym}}(A_{\text{bp}}(X))$ are

  * •

located in $[0,2]$,

  * •

symmetric around 1 meaning for every eigenvalue $\lambda=1-\eta$, $\eta\geq 0$
there is a corresponding eigenvalue at $1+\eta$,

  * •

given by $\\{1\pm\sigma_{i}(T_{\text{sym}}(X))\\}_{i=1}^{\min(R,C)}$ and
$R+C-2\min(R,C)$ 1s.

The singular values of $T_{\text{sym}}(X)$ are located in $[0,1]$.

Let $U\in\mathbb{R}^{R\times K}$, $V\in\mathbb{R}^{C\times K}$ be matrices of
the largest $K$ left and right singular vectors of $T_{\text{sym}}(X))$
respectively. Then $\begin{bmatrix}U\\\ V\end{bmatrix}$ is the matrix of the
smallest $K$ eigenvectors of $L_{\text{sym}}(A_{\text{bp}}(X))$ and
$\begin{bmatrix}U\\\ -V\end{bmatrix}$ is the matrix of the matrix of the
largest $K$ eigenvectors respectively.

###### Proof.

of Proposition G.1

Consider the subgraph corresponding to the $m$ vertices of this graph which
are contained in connected components with at least two vertices. Without loss
of generality assume that these are the first $m$ nodes of the graph and let
$\widetilde{A}\in\mathbb{R}^{m\times m}$ be the corresponding adjacency
matrix. Then

$\text{L}_{\text{sym}}(A)=\begin{bmatrix}\text{L}_{\text{sym}}(\widetilde{A})&0_{m\times(n-m)}\\\
0_{(n-m)\times m}&I_{n-m}\end{bmatrix}$

From this we see the unit vectors $e_{i}\in\mathbb{R}^{m}$ for $i=n+1,\dots,m$
are eigenvectors of $\text{L}_{\text{sym}}(A)$ with eigenvalue 1 and the claim
about isolated vertices follows.

If $(\lambda,\widetilde{v})$ with $\widetilde{v}\in\mathbb{R}^{m}$ is an
eigenvalue/eigenvector pair for $\text{L}_{\text{sym}}(\widetilde{A})$ then
$(\lambda,v)$ is an eigenvalue/eigenvector pair for $\text{L}_{\text{sym}}(A)$
with $v:=(\widetilde{v},0,\dots,0)\in\mathbb{R}^{n}$. Proposition 4 of Von
Luxburg, (2007) holds for $\text{L}_{\text{sym}}(\widetilde{A})$ which has no
isolated vertices. Therefore, an orthonormal set of $m$ eigenvectors of
$\text{L}_{\text{sym}}(\widetilde{A})$ give $m$ orthonormal eigenvectors of
$\text{L}_{\text{sym}}(A)$ with the same eigenvalues. From this we see the
claim about 0 eigenvalues of $\text{L}_{\text{sym}}(A)$ and the corresponding
eigenspace follows.

∎

###### Proof.

of Proposition G.2

By checking $\text{L}_{\text{sym}}(A_{\text{bp}}(X))$ is diagonally dominant
we see it is positive-semi definite. The diagonal elements of
$\text{L}_{\text{sym}}(A_{\text{bp}}(X))$ are equal to 1. Consider the first
row; if $\text{deg}(A_{\text{bp}}(X))_{1}=0$ then the first row of
$\text{L}_{\text{sym}}(A)$ is equal to the first standard basis vector. If
$\text{deg}(A_{\text{bp}}(X))_{1}>0$ then
$\text{L}_{\text{sym}}(A_{\text{bp}}(X))_{1j}=\frac{A_{\text{bp}}(X)_{1j}}{\sum_{i=1}^{n}A_{\text{bp}}(X)_{ij}}\leq
1$.

We see that

$\text{L}_{\text{sym}}(A_{\text{bp}}(X))=I-\begin{bmatrix}0&T_{\text{sym}}(X)\\\
T_{\text{sym}}(X)^{T}&0\end{bmatrix}.$

Note that the spectrum of the second matrix on the right hand side is
symmetric around 0. It is straightforward to check the remaining claims of the
proposition.

∎

###### Proof.

of Proposition 4.1

By inspecting the adjacency matrix $A_{\text{bp}}(X)$ it is clear there is a
one-to-one correspondence between the zero rows/columns of $X$ and the
isolated vertices in $G(X)$. Without loss of generality, assume $G(X)$ has no
isolated vertices and suppose $G(X)$ has $B$ connected components with at
least two vertices.

Let $\sigma_{r}$ be a permutation of the rows of $X$ such that the first rows
of $X$ belong to the first connected component, the next rows of $X$ belong to
the second connected component, etc. Let $\sigma_{c}$ be the analogous
permutation of the columns of $X$. Let $\widetilde{X}$ be the result of
applying these two permutations to $X$, then $\widetilde{X}$ is block diagonal
with $B$ blocks. We thus conclude that the number of connected components of
$G$ is a lower bound for $NB(X)$.

Now suppose there exists a permutation, $\sigma_{r}$ of the rows of $X$ and a
permutation $\sigma_{c}$ of the columns of $X$ such that the resulting matrix
has $C\geq B+1$ blocks. Let $\widetilde{X}$ be the result of applying these
two permutations to $X$, then $A_{\text{bp}}(\widetilde{X})$ is the adjacency
matrix of a bipartite graph with $C$ connected components. But this is a
contradiction, because shuffling the node labels of $G(X)$ induces a graph
isomorphism so the number of connected components must be conserved. Thus
claims 1 and 2 are equivalent.

Claims 2 and 3 are equivalent by Proposition G.1. Claims 1 and 4 are
equivalent because $G(X)$ has $B+Z_{\text{row}}+Z_{\text{col}}$ connected
components and Proposition 2 of Von Luxburg, (2007).

∎

###### Proof.

of Corollary B.1. Note that
$\text{dim}(\text{ker}(\text{diag}(\text{deg}(A_{\text{bp}}(X)))))=R+C-(\widetilde{R}+\widetilde{C})$.
Eigenvectors of $L_{\text{sym}}(A_{\text{bp}}(X))$ that live in the kernel of
$\text{diag}(\text{deg}(A_{\text{bp}}(X)))$ correspond to isolated nodes thus
give an eigenvalue of $1$. By Proposition G.1,
$\lambda_{(k)}(L_{\text{sym}}(A_{\text{bp}}(X)))\leq 1$ for $k\leq\min(R,C)$.
Therefore, for $k\leq\widetilde{R}+\widetilde{C}$, none of the
$\lambda_{(k)}(L_{\text{sym}}(A_{\text{bp}}(X)))$ correspond to eigenvectors
in the kernel of $\text{diag}(\text{deg}(A_{\text{bp}}(X)))$. The result
follows from Proposition B.1. ∎

###### Proof.

of Proposition A.1 The proof of Proposition 4.1 can be generalized to multi-
arrays. ∎

###### Proof.

of Corollary C.1 The result follows from Propositions B.2 and G.2.
$K\leq\widetilde{R}+\widetilde{C}$ ensures the assumption of Proposition B.2
are satisfied for case 3. ∎

###  G.4 Block diagonal optimization problem solution sets

###### Proof.

of Proposition 4.2

Problems (12) and (13) are equivalent by Proposition 4.1.

Suppose $(X,U)$ are such a global minimizer of (15). By Proposition 4.1,
$NB(X)\geq B$ so $X$ satisfies the constraints of (12). Assume for the sake of
contradiction that $Y$ is a better minimizer of (12) i.e. $f(Y)<f(X)$. Let
$U_{Y}$ be the smallest $B$ generalized eigenvectors of
$(L_{\text{un}}(A_{\text{bp}}(Y),\text{diag}(\text{deg}(A_{\text{bp}}(Y))))$.
Then

$\displaystyle
f(Y)+\alpha\text{Tr}\left(U_{Y}^{T}L_{\text{un}}(A_{\text{bp}}(Y))U_{Y}\right)=f(Y)+\alpha\sum_{j=1}^{B}\lambda_{(j)}\left(L_{\text{sym}}(A_{\text{bp}}(Y))\right)=f(Y)$
(56)
$\displaystyle<f(X)=f(X)+\alpha\sum_{j=1}^{B}\lambda_{(j)}\left(L_{\text{sym}}(A_{\text{bp}}(X))\right)=f(X)+\alpha\text{Tr}\left(U^{T}L_{\text{un}}(A_{\text{bp}}(X))U\right).$

Thus $(Y,U_{Y})$ is a better minimizer of (15) contradicting the fact that
$(X,U)$ is a global minimizer. ∎

###### Proof.

of Proposition C.1

Suppose $X$ satisfies the conditions of claim 1. Then $NB(X)\geq B$ by
Proposition 4.1 so $X$ satisfies the constraints of (12). Suppose $Y$ is a
better solution for (12), then
$\sum_{j=1}^{B}\lambda_{(j)}\left(L_{\text{sym}}(A_{\text{bp}}(Y))\right)=0$
and $f(Y)<f(X)$. Thus $Y$ is a better solution for (14).

Suppose $X$ conditions of claim 2. For fixed $X$, $U_{X}$ is a coordinate-wise
minimizer of (15) if and only if the columns of $U_{X}$ are the smallest $B$
generalized eigenvectors of
$(L_{\text{un}}(A_{\text{bp}}(X),\text{diag}(\text{deg}(A_{\text{bp}}(X))))$
and

$\text{Tr}\left(U_{X}^{T}L_{\text{un}}(A_{\text{bp}}(X))U_{X}\right)=\sum_{j=1}^{B}\lambda_{(j)}\left(L_{\text{sym}}(A_{\text{bp}}(X))\right).$

Note the row sum condition on $X$ ensures the kernel constraints of
Proposition B.2 hold. Suppose $(Y,\widetilde{U})$ are a better solution to
(15) than $(X,U_{X})$. Then $(Y,U_{Y})$ are also a better solution to (15)
than $(X,U_{X})$. But then $Y$ is a better solution to (14).

∎

###  G.5 Algorithm convergence

This section applies Zangwill’s global convergence theorem (Zangwill,, 1969)
to prove Proposition C.2. Following Sriperumbudur and Lanckriet, (2009), a
point-to-set mapping $\mathcal{A}:\mathcal{X}\to 2^{\mathcal{Y}}$ assigns a
subset $\mathcal{A}(x)\subseteq\mathcal{Y}$ to a point $x\in\mathcal{X}$. A
point-to-set mapping is closed if $x_{k}\to x^{*},y_{k}\to
y^{*},y_{k}\in\mathcal{A}(x_{k})$ together imply $y^{*}\in\mathcal{A}(x^{*})$;
this is a generalization of continuity for functions. A generalized fixed
point of $\mathcal{A}:\mathcal{X}\to 2^{\mathcal{X}}$ is a point $x$ such that
$x\in\mathcal{A}(x)$.

######  Lemma G.1.

Let $f:A\to B$ and $g:B\to C$ be closed, non-empty point-to-set maps. If $B$
is sequentially compact then $g\circ f:A\to C$ is closed.

###### Proof.

of Lemma G.1 Let $a_{k}\to a$, $c_{k}\to c$ and $c_{k}\in g(f(a_{k}))$. Let
$b_{k}\in f(a_{k})\cap g^{-1}(c_{k})$ where the inverse denotes the set of
pre-images. By assumption on $B$ there exists a convergent subsequence
$\\{b_{k_{i}}\\}_{i=1}^{\infty}$ such that $b_{k_{i}}\to b$ for some $b\in B$.
Since $f$ is closed, $b\in f(a)$; since $g$ is closed $c\in g(b)$, therefore
$c\in g(f(a))$. ∎

Let $\mathcal{A}(X,U)$ be the point-to-set map corresponding to Algorithm 1.
$\mathcal{A}:=\textsc{update-X}\circ\textsc{update-U}$ where
$\textsc{update}-U(X)$ solves Problem (24) for fixed $X$ and update
$\textsc{update}-X$ is either the full update (Assumption C.3.1) or surrogate
update (Assumption C.3.2.)

Let $\psi(X)$ be the objective function in (23) and let $\phi(X,U)$ be the
objective function in (24). Each step of Algorithm 1 decreases these objective
functions.

######  Proposition G.3.

Let $(X^{*},U^{*})\in\mathcal{A}(X,U)$ then $\phi(X^{*},U^{*})\leq\phi(X,U)$
and $\psi(X^{*})\leq\psi(X)$.

######  Lemma G.2.

Under Assumption C.3.1, if $(X,U)$ is a generalized fixed point of
$\mathcal{A}$ then $(X,U)\in\mathcal{LG}$.

Under Assumption C.3.2, if $(X,U)$ is a generalized fixed point of
$\mathcal{A}$ then $(X,U)\in\mathcal{SG}$.

###### Proof.

of Lemma G.2 By construction $U$ is a global minimizer of (24) for fixed $X$.

Under Assumption C.3.1, $X$ is a global minimizer of (24) for fixed $V$ thus
the first claim follows.

The constraints of (31)/(32) are affine in $X$ so the KKT conditions are first
order necessary. Under Assumption C.3.2 if $X^{*}$ is a generalized fixed
point of $\mathcal{A}$, $X^{*}$ is a minimizer of (32) and thus satisfies the
KKT conditions for (32). Because $Q$ and $f$ have the same first order
behavior by Assumption C.3.2, a KKT point of (32) is also a KKT point of (31)
and the second claim follows. ∎

###### Proof.

of Proposition C.2 We apply Zangwill’s global convergence theorem (Zangwill,,
1969) by checking the three conditions for Theorem 2 of Sriperumbudur and
Lanckriet, (2009). Let $\phi(X,U)$ be the objective function of (24) and let
$\Gamma$ be the set of generalized fixed points of $\mathcal{A}$.

Let $S_{X^{0}}$ denote the compact set in Assumption C.1. By Propositions B.2
and G.2 the second term in the objective function of (24) is upper bounded by
$\alpha w^{T}\mathbf{1}_{K}$. Thus, by Proposition G.3 the iterates
$\\{X^{s}\\}_{s=0}^{\infty}\subseteq S_{X^{0}}$.

By assumption C.2, without loss of generality we can add the constraint
$\text{deg}(A_{\text{bp}}(X))\geq\eta\mathbf{1}_{R+C}$ to (24). Since
$\eta>0$, the constraint given by (29) implies that at any solution
$||U_{k}||_{2}^{2}\leq\frac{1}{\eta}$ where the inequality is applied element
wise. Thus the iterates $\\{U^{s}\\}_{s=0}^{\infty}\subseteq
S_{\eta}:=\\{U|||U_{k}||_{2}^{2}\leq\frac{1}{\eta},k\in[K]\\}$ which is a
compact set. We conclude $\\{(X^{s},U^{s})\\}_{s=0}^{\infty}\subseteq
S:=S_{X^{0}}\times S_{\eta}$ a compact set and condition (1) holds.

By construction of $\mathcal{A}$ and Proposition G.3, condition (2) holds.

If $(X^{\prime},U^{\prime})=\mathcal{A}(X,U)$ then the constraint sets for the
$U$ and $X$ update problem starting from $(X^{\prime},U^{\prime})$ are non-
empty. Therefore, the $X$ update and $U$ update steps are non-empty that
compose to make $\mathcal{A}$ are non-empty. By the above discussion, the $U$
update always lives in the compact set $S_{\eta}$. Therefore, Lemma G.1 shows
$\mathcal{A}$ is closed and condition (3) holds.

Therefore, by Theorem 2 of Sriperumbudur and Lanckriet, (2009), all the limit
points of $\\{(X^{s},U^{s})\\}_{s=0}^{\infty}$ are the generalized fixed
points of $\mathcal{A}$ and
$\lim_{s\to\infty}\phi(X^{s},U^{s})=\lim_{s\to\infty}\phi(X^{*},U^{*})$ where
$(X^{*},U^{*})$ is some generalized fixed point of $\mathcal{A}$.

The result now follows from Lemma G.2. ∎

###### Proof.

of Proposition C.3 First implication follows from proof by contradiction.
Consider Problem (24) where $U$ is fixed and the constraint set is non-empty.
For all such fixed $U$, the constraints on $X$ are affine thus Linear
Constraint Qualification holds thus the KKT conditions are first order
necessary. ∎

###  G.6 Block diagonal MVMM

###### Proof.

of Proposition E.1 Let $f(x)=-a\log(x+\epsilon)+bx$. Then $f(x)$ is strictly
convex on $(-\epsilon,\infty)$. Setting
$0=f^{\prime}(x)=\frac{-1}{x+\epsilon}+b$ leaves us with
$x=\frac{a}{b}-\epsilon$ for the unique stationary point of $f(x)$. If
$\frac{a}{b}-\epsilon>0$, this must be the minimizer of (46). Otherwise, $x=0$
must be the minimizer. ∎

## References

  * Asteris et al., (2015) Asteris, M., Papailiopoulos, D., Kyrillidis, A., and Dimakis, A. G. (2015).  Sparse pca via bipartite matchings.  In Advances in Neural Information Processing Systems, pages 766–774.
  * Bickel and Scheffer, (2004) Bickel, S. and Scheffer, T. (2004).  Multi-view clustering.  In ICDM, volume 4, pages 19–26.
  * Boyd et al., (2004) Boyd, S., Boyd, S. P., and Vandenberghe, L. (2004).  Convex optimization.  Cambridge university press.
  * Bugdary and Maymon, (2019) Bugdary, S. and Maymon, S. (2019).  Online clustering by penalized weighted gmm.  arXiv preprint arXiv:1902.02544.
  * Carmichael et al., (2019) Carmichael, I., Calhoun, B. C., Hoadley, K. A., Troester, M. A., Geradts, J., Couture, H. D., Olsson, L., Perou, C. M., Niethammer, M., Hannig, J., et al. (2019).  Joint and individual analysis of breast cancer histologic images and genomic covariates.  arXiv preprint arXiv:1912.00434.
  * Chen and Chen, (2008) Chen, J. and Chen, Z. (2008).  Extended bayesian information criteria for model selection with large model spaces.  Biometrika, 95(3):759–771.
  * Chen et al., (2013) Chen, K., Dong, H., and Chan, K.-S. (2013).  Reduced rank regression via adaptive nuclear norm penalization.  Biometrika, 100(4):901–920.
  * Davidson-Pilon et al., (2020) Davidson-Pilon, C. et al. (2020).  lifelines.  https://github.com/CamDavidsonPilon/lifelines.
  * 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.
  * Devijver and Gallopin, (2018) Devijver, E. and Gallopin, M. (2018).  Block-diagonal covariance selection for high-dimensional gaussian graphical models.  Journal of the American Statistical Association, 113(521):306–314.
  * 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, pages 269–274.
  * Diamond and Boyd, (2016) Diamond, S. and Boyd, S. (2016).  Cvxpy: A python-embedded modeling language for convex optimization.  The Journal of Machine Learning Research, 17(1):2909–2913.
  * Domahidi et al., (2013) Domahidi, A., Chu, E., and Boyd, S. (2013).  Ecos: An socp solver for embedded systems.  In 2013 European Control Conference (ECC), pages 3071–3076. IEEE.
  * Fan, (1949) Fan, K. (1949).  On a theorem of weyl concerning eigenvalues of linear transformations i.  Proceedings of the National Academy of Sciences of the United States of America, 35(11):652.
  * Feng et al., (2014) Feng, J., Lin, Z., Xu, H., and Yan, S. (2014).  Robust subspace segmentation with block-diagonal prior.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 3818–3825.
  * Feng et al., (2018) Feng, Q., Jiang, M., Hannig, J., and Marron, J. (2018).  Angle-based joint and individual variation explained.  Journal of multivariate analysis, 166:241–265.
  * Fu and Perry, (2020) Fu, W. and Perry, P. O. (2020).  Estimating the number of clusters using cross-validation.  Journal of Computational and Graphical Statistics, 29(1):162–173.
  * Gabasova et al., (2017) Gabasova, E., Reid, J., and Wernisch, L. (2017).  Clusternomics: Integrative context-dependent clustering for heterogeneous datasets.  PLoS computational biology, 13(10):e1005781.
  * (19) Gao, L., Bien, J., and Witten, D. (2019a).  Are clusterings of multiple data views independent?  Biostatistics (Oxford, England).
  * (20) Gao, L. L., Witten, D., and Bien, J. (2019b).  Testing for association in multi-view network data.  arXiv preprint arXiv:1909.11640.
  * Gavish and Donoho, (2014) Gavish, M. and Donoho, D. L. (2014).  The optimal hard threshold for singular values is 4/sqrt(3).  IEEE Transactions on Information Theory, 60(8):5040–5053.
  * Gaynanova and Li, (2017) Gaynanova, I. and Li, G. (2017).  Structural learning and integrative decomposition of multi-view data.  arXiv preprint arXiv:1707.06573.
  * Gouwens et al., (2020) Gouwens, N. W., Sorensen, S. A., Baftizadeh, F., Budzillo, A., Lee, B. R., Jarsky, T., Alfiler, L., Arkhipov, A., Baker, K., Barkan, E., et al. (2020).  Toward an integrated classification of neuronal cell types: morphoelectric and transcriptomic characterization of individual gabaergic cortical neurons.  BioRxiv.
  * Gouwens et al., (2019) Gouwens, N. W., Sorensen, S. A., Berg, J., Lee, C., Jarsky, T., Ting, J., Sunkin, S. M., Feng, D., Anastassiou, C. A., Barkan, E., et al. (2019).  Classification of electrophysiological and morphological neuron types in the mouse visual cortex.  Nature neuroscience, 22(7):1182–1195.
  * Gu et al., (2014) Gu, S., Zhang, L., Zuo, W., and Feng, X. (2014).  Weighted nuclear norm minimization with application to image denoising.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 2862–2869.
  * Gu and Xu, (2019) Gu, Y. and Xu, G. (2019).  Learning attribute patterns in high-dimensional structured latent attribute models.  Journal of Machine Learning Research, 20(115):1–58.
  * Han et al., (2017) Han, J., Song, K., Nie, F., and Li, X. (2017).  Bilateral k-means algorithm for fast co-clustering.  In Proceedings of the Thirty-First AAAI Conference on Artificial Intelligence, pages 1969–1975.
  * Hellton and Thoresen, (2016) Hellton, K. H. and Thoresen, M. (2016).  Integrative clustering of high-dimensional data with joint and individual clusters.  Biostatistics, 17(3):537–548.
  * Hoadley et al., (2018) Hoadley, K. A., Yau, C., Hinoue, T., Wolf, D. M., Lazar, A. J., Drill, E., Shen, R., Taylor, A. M., Cherniack, A. D., Thorsson, V., et al. (2018).  Cell-of-origin patterns dominate the molecular classification of 10,000 tumors from 33 types of cancer.  Cell, 173(2):291–304.
  * Hoda et al., (2020) Hoda, S. A., Rosen, P. P., Brogi, E., and Koerner, F. C. (2020).  Rosen’s breast pathology.  Lippincott Williams & Wilkins.
  * Hotelling, (1936) Hotelling, H. (1936).  Relations between two sets of variates.  Biometrika, 28(3/4):321–377.
  * Huang et al., (2017) Huang, T., Peng, H., and Zhang, K. (2017).  Model selection for gaussian mixture models.  Statistica Sinica, pages 147–169.
  * Hunter, (2007) Hunter, J. D. (2007).  Matplotlib: A 2d graphics environment.  Computing in science & engineering, 9(3):90–95.
  * Kirk et al., (2012) Kirk, P., Griffin, J. E., Savage, R. S., Ghahramani, Z., and Wild, D. L. (2012).  Bayesian correlated clustering to integrate multiple datasets.  Bioinformatics, 28(24):3290–3297.
  * Klami et al., (2014) Klami, A., Virtanen, S., Leppäaho, E., and Kaski, S. (2014).  Group factor analysis.  IEEE transactions on neural networks and learning systems, 26(9):2136–2147.
  * Kumar et al., (2011) Kumar, A., Rai, P., and Daume, H. (2011).  Co-regularized multi-view spectral clustering.  In Advances in neural information processing systems, pages 1413–1421.
  * Kumar et al., (2019) Kumar, S., Ying, J., Cardoso, J. V. d. M., and Palomar, D. (2019).  A unified framework for structured graph learning via spectral constraints.  arXiv preprint arXiv:1904.09792.
  * Lange et al., (2000) Lange, K., Hunter, D. R., and Yang, I. (2000).  Optimization transfer using surrogate objective functions.  Journal of computational and graphical statistics, 9(1):1–20.
  * Liu et al., (2018) Liu, J., Lichtenberg, T., Hoadley, K. A., Poisson, L. M., Lazar, A. J., Cherniack, A. D., Kovatich, A. J., Benz, C. C., Levine, D. A., Lee, A. V., et al. (2018).  An integrated tcga pan-cancer clinical data resource to drive high-quality survival outcome analytics.  Cell, 173(2):400–416.
  * Lock and Dunson, (2013) Lock, E. F. and Dunson, D. B. (2013).  Bayesian consensus clustering.  Bioinformatics, 29(20):2610–2616.
  * Lock et al., (2013) Lock, E. F., Hoadley, K. A., Marron, J. S., and Nobel, A. B. (2013).  Joint and individual variation explained (jive) for integrated analysis of multiple data types.  The annals of applied statistics, 7(1):523.
  * Lu et al., (2018) Lu, C., Feng, J., Lin, Z., Mei, T., and Yan, S. (2018).  Subspace clustering by block diagonal representation.  IEEE transactions on pattern analysis and machine intelligence, 41(2):487–501.
  * Marlin and Murphy, (2009) Marlin, B. M. and Murphy, K. P. (2009).  Sparse gaussian graphical models with unknown block structure.  In Proceedings of the 26th Annual International Conference on Machine Learning, pages 705–712.
  * Marshall et al., (1979) Marshall, A. W., Olkin, I., and Arnold, B. C. (1979).  Inequalities: theory of majorization and its applications, volume 143.  Springer.
  * McKinney et al., (2010) McKinney, W. et al. (2010).  Data structures for statistical computing in python.  In Proceedings of the 9th Python in Science Conference, volume 445, pages 51–56. Austin, TX.
  * Network et al., (2012) Network, C. G. A. et al. (2012).  Comprehensive molecular portraits of human breast tumours.  Nature, 490(7418):61.
  * Nie et al., (2017) Nie, F., Wang, X., Deng, C., and Huang, H. (2017).  Learning a structured optimal bipartite graph for co-clustering.  In Advances in Neural Information Processing Systems, pages 4129–4138.
  * Nie et al., (2016) Nie, F., Wang, X., Jordan, M. I., and Huang, H. (2016).  The constrained laplacian rank algorithm for graph-based clustering.  In Thirtieth AAAI Conference on Artificial Intelligence.
  * Parker et al., (2009) Parker, J. S., Mullins, M., Cheang, M. C., Leung, S., Voduc, D., Vickery, T., Davies, S., Fauron, C., He, X., Hu, Z., et al. (2009).  Supervised risk predictor of breast cancer based on intrinsic subtypes.  Journal of clinical oncology, 27(8):1160.
  * Pedregosa et al., (2011) Pedregosa, F. et al. (2011).  Scikit-learn: Machine learning in python.  the 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.
  * Razaviyayn et al., (2013) Razaviyayn, M., Hong, M., and Luo, Z.-Q. (2013).  A unified convergence analysis of block successive minimization methods for nonsmooth optimization.  SIAM Journal on Optimization, 23(2):1126–1153.