# Cluster Quilting: Spectral Clustering for Patchwork Learning

Lili Zheng† Corresponding author<EMAIL_ADDRESS>Department of
Statistics, University of Illinois Urbana-Champaign  Andersen Chang†
Department of Neuroscience, Baylor College of Medicine  Genevera I. Allen
Department of Statistics, Columbia University

###### Abstract

Patchwork learning arises as a new and challenging data collection paradigm
where both samples and features are observed in fragmented subsets. Due to
technological limits, measurement expense, or multimodal data integration,
such patchwork data structures are frequently seen in neuroscience,
healthcare, and genomics, among others. Instead of analyzing each data patch
separately, it is highly desirable to extract comprehensive knowledge from the
whole data set. In this work, we focus on the clustering problem in patchwork
learning, aiming at discovering clusters amongst all samples even when some
are never jointly observed for any feature. We propose a novel spectral
clustering method called Cluster Quilting, consisting of (i) patch ordering
that exploits the overlapping structure amongst all patches, (ii) patchwise
SVD, (iii) sequential linear mapping of top singular vectors for patch
overlaps, followed by (iv) k-means on the combined and weighted singular
vectors. Under a sub-Gaussian mixture model, we establish theoretical
guarantees via a non-asymptotic misclustering rate bound that reflects both
properties of the patch-wise observation regime as well as the clustering
signal and noise dependencies. We also validate our Cluster Quilting algorithm
through extensive empirical studies on both simulated and real data sets in
neuroscience and genomics, where it discovers more accurate and scientifically
more plausible clusters than other approaches.

Keywords— Spectral clustering, blockwise missingness, patchwork learning,
mixture models, multi-omics, data integration, matrix completion
††footnotetext: $\dagger$: Equal contribution.

##  1 Introduction

Clustering is a fundamental statistical tool for discovering prominent
patterns from large-scale data sets, such as finding gene pathways (Dalton et
al., , 2009), tumor cell subtypes (Shen et al., , 2009), text topics (Aggarwal
and Zhai, , 2012), communities in social networks (Mishra et al., , 2007), and
many others. Various clustering techniques were proposed in the literature and
widely applied, e.g., k-means (MacQueen, , 1967), hierarchical clustering
(Bridges, , 1966), convex clustering (Chen et al., , 2015), and spectral
clustering (Donath and Hoffman, , 1972). Although clustering for fully
observed data has been extensively studied, in this paper, we consider a novel
setting where one only has access to a collection of data patches, i.e.,
partial observations of both features and samples (see e.g., Fig. 1). This is
motivated by many biomedical applications. For example, in neuroscience,
genomics, and proteomics, simultaneous measurement of all observations of
interest can be extremely expensive or impossible due to technological
constraints; thus, in these situations, data may be collected such that only a
subset of features and observations are recorded within a single experimental
session (MICrONS Consortium, , 2021; Browning, , 2008; Baghfalaki et al., ,
2016). Another example arises in healthcare applications, where data collected
on patients from different sites or hospitals may only cover distinct subsets
of all possible modalities due to clinician discretion or availability of
facilities and thus lead to a lack of measurements for entire blocks of
observations and features (Cai et al., , 2019; Kline et al., , 2022). This
particular data missingness problem was introduced as the new “patchwork
learning paradigm” in Rajendran et al.,  (2023); there, the authors also
outline numerous opportunities and challenges in this area, one of which is
the fragmented observational structure where missingness is non-random. Given
such patchwise observations, naively applying existing clustering techniques
to each patch fails to give a holistic view of the whole population. Instead,
we aim to effectively utilize all data patches and perform a comprehensive
clustering for all samples available. We term this as the “cluster quilting”
problem, following a recent work on “graph quilting” (Vinci et al., , 2019)
which aims to learn graph structure from a similar type of data.

###  1.1 Problem Set-up

To fix ideas, we start by introducing the commonly studied mixture model for
data generation and the patchwork data observation setting. Suppose there are
$K>1$ clusters with $p$-dimensional centroids
$\theta^{*}_{1},\dots,\theta^{*}_{K}$. For each sample $1\leq i\leq n$, let
$z_{i}\in[K]$ be its cluster label, and its associated full feature vector
$X_{i}=\theta_{z_{i}}+\epsilon_{i}\in\mathbb{R}^{p}$ with mean-zero random
noise $\epsilon_{i}$. Let $X=(X_{1},\cdots,X_{n})\in\mathbb{R}^{p\times n}$ be
the data matrix,
$\Theta^{*}=(\theta_{1},\cdots,\theta_{K})\in\mathbb{R}^{p\times K}$ be the
centroid matrix, $E=(\epsilon_{1},\dots,\epsilon_{n})\in\mathbb{R}^{p\times
K}$ be the noise matrix, and $F^{*}\in\mathbb{R}^{n\times K}$ encode the
cluster membership: $F^{*}_{i,k}={{\mathbbm{1}}}\\{z_{i}=k\\}$. The model can
be written in the following matrix format:

$X=X^{*}+E=\Theta^{*}F^{*\top}+E\in\mathbb{R}^{p\times n},$ (1)

where $X^{*}=\mathbb{E}X$ is the expectation of observational data. Our goal
is then to assign the $n$ samples into $K$ clusters, or equivalently, to
estimate the cluster labels $z_{1},\dots,z_{n}$ up to permutation on $[K]$.
However, in the patchwork observation setting, we do not have access to the
full data matrix $X$ but rather only subsets of the rows and columns
structured as patches. More formally, we suppose here that all $p$ features
can be divided into $M$ disjoint subsets $T_{1},\,\dots,\,T_{M}$,
$\cup_{m=1}^{M}T_{m}=[p]$ and that each feature subset $T_{m}$ of size $p_{m}$
is associated with a corresponding observation subset $I_{m}\subset[n]$ of
size $n_{m}$; the observed data thus consists of submatrices of the full data:
$\\{X_{T_{1},I_{1}},\dots,X_{T_{M},I_{M}}\\}$. Following the prior work on
patchwork learning (Rajendran et al., , 2023), we also use a graph to
characterize the observational relationship between different feature blocks
(e.g., modalities). In particular, consider a graph $\mathcal{G}=(V,E)$ where
the node set $V=[M]$ are the $M$ patches, and $(j,k)\in E$ if and only if
$I_{j}\cap I_{k}\neq\emptyset$ (see Fig. 1 for illustration). Throughout the
paper, we assume the graph $\mathcal{G}$ to be connected so that the signals
can be shared across different patches.

(a) Sequential patches motivated by neuronal functional recording via calcium
imaging, where feature subsets are distinct time periods and samples are
neurons.

(b) Mosaic patches motivated by healthcare and genomics, where disjoint
feature subsets cover different modalities.

Figure 1: Two examples of patchwork observational patterns. The graphs above
represent the overlapping relationships between disjoint feature blocks.

###  1.2 Related Literature

Many different approaches have previously been proposed for clustering with
missing data. Classically, two-step approaches have been applied in which
missing cases are either deleted or imputed (Inman et al., , 2015; Aktaş et
al., , 2019); the former can not be applied in the patchwork learning scenario
because there do not exist any fully observed rows or columns, while the
performance of latter does not come with any statistical guarantees about
clustering performance. More recently, several papers have proposed general
likelihood-based optimization approaches for estimating clusters with missing
data. Examples of these include a constrained k-means objective estimated by
majorization-minimization (Chi et al., , 2016), an $\ell_{0}$ fusion-penalty
based optimization method Poddar and Jacob,  (2019), a likelihood-based
approach to estimate mixture models parameters from missing data Hunt and
Jorgensen,  (2003); Lin et al.,  (2006), a modified fuzzy c-means approach
with an adjusted objective to account for data missingness Hathaway and
Bezdek,  (2001), and a Bayesian model-based random labeling point process
approach Boluki et al.,  (2019).

However, one key difference between the setting used in many prior papers and
the one we consider in this paper is in the structure of the missing entries.
Cluster quilting is particularly a non-trivial task because of the imposed
block observation scheme, as the missingness pattern is highly non-random and
very dense. Additionally, many pairs of samples are only observed on disjoint
sets of features, meaning that there is no direct way to compute similarities
between many of the pairs of samples in the data. Although imputation comes up
as a natural idea, how to accurately impute the data or similarity matrix from
non-random observational patterns is still a challenging and active research
area (Chatterjee, , 2020; Yan and Wainwright, , 2024; Choi and Yuan, , 2023),
where most prior works consider missing patterns different from ours; it is
also unclear whether the imputed data suffices for the downstream clustering
task in the patchwork observation setting. Therefore, previous results on
theoretical guarantees and empirical performance of the likelihood-based
approaches outlined above may not be applicable here. In particular, many of
the previous works lack provable theoretical computational convergence, while
others only study the case of random missingness in the data matrix. Some
prior existing methods formulate cluster quilting as various optimization
problems (Wang et al., , 2019; Wen et al., , 2022), but these methods lack
statistical guarantees of convergence, relying on non-convex optimization with
no assurance of global algorithmic convergence.

###  1.3 Contribution and Organization

In this paper, we propose a new spectral-based method for clustering
incomplete data called the Cluster Quilting algorithm; this is introduced in
Section 2. Broadly, our Cluster Quilting procedure consists of the following
steps: (i) order the patches based on their overlapping signals; (ii) learn
the top singular vectors of each data patch; (iii) perform sequential matching
across data patches in the sample space; (iv) obtain an estimate of the full
top singular spaces and apply k-means on the singular value-weighted singular
vectors. While much work has been done regarding the theoretical and empirical
properties of spectral clustering under the stochastic block (Rohe et al., ,
2011; Jin, , 2015; Lei and Rinaldo, , 2015; Zhou and Chen, , 2023) and
Gaussian mixture (Löffler et al., , 2021; Abbe et al., , 2022) data generation
models which we study in this work, we are the first to study spectral
clustering for the general patchwork learning observation setting. Our method
is most similar to those proposed in Bishop and Yu,  (2014), Yu et al.,
(2020), and Chang et al.,  (2022), which also utilize a low-rank assumption in
order to perform matrix completion for downstream analysis and modeling in the
patchwork learning context. However, our work fundamentally differs from these
previous works in terms of the considered data-generative model and the
clustering objective, the general patchwork observations that go beyond
sequential observational pattern (Chang et al., , 2022), as well as our new
method that starts with a theory-motivated, data-driven patch ordering step.
In particular, we provide novel theory for our cluster quilting algorithm in
Section 3, characterizing the conditions and misclustering rate in terms of
separation distance or minimum singular values in the data patches; the result
presents a much tighter bound compared to the previous literature on matrix
completion (Bishop and Yu, , 2014). We also prove that the data-driven
selection of patch ordering guarantees a similar performance to the oracle
ordering, with only a constant factor inflation in the SNR condition and
misclustering rate. Lastly, we validate our method through extensive empirical
studies on simulated and real data sets in Sections 4-5, in which we show that
Cluster Quilting improves upon previous methods proposed for incomplete
spectral clustering in the task of recovering true cluster labels in the
patchwork observation setting.

##  2 Cluster Quilting Algorithm

Inspired by the success of spectral clustering on sub-Gaussian mixture models
(Löffler et al., , 2021), we propose a spectral algorithm called _Cluster
Quilting_ , where the key steps are learning the top right singular vectors of
$X^{*}$. Due to the blockwise missingness, the full right singular vectors of
$X^{*}$ cannot be recovered from any data patch alone; but instead, we need to
combine all pieces together and match them using linear transformations.
Denote by $r>0$ the rank of the ground truth matrix $X^{*}$. We define the SVD
of $X^{*}$ as $X^{*}=U^{*}\Lambda^{*}V^{*\top}$ for orthogonal matrices
$U\in\mathbb{O}^{p\times r}$, $V\in\mathbb{R}^{n\times r}$, and diagonal
matrix $\Lambda\in\mathbb{R}^{r\times r}$. Since we do not assume the cluster
centroid matrix $\Theta^{*}$ to be of full rank, the rank $r\leq K$ does not
necessarily equal $K$. For now, we assume the data patches are appropriately
ordered and each sample subset
$I_{m}\cap(\cup_{m^{\prime}<m}I_{m^{\prime}})\neq\emptyset$; in fact, such an
order always exists under our assumption that the observational graph
$\mathcal{G}$ is connected. We will introduce several theory-inspired
strategies to choose the ordering in Section 2.1.

Given a fixed ordering, our algorithm consists of three main steps. First, we
perform singular value decomposition for each patch $X_{T_{m},I_{m}}$ and take
its top $r$ singular vectors. We then find the best invertible matrix to match
the eigenvectors of the $m$th patch with the previous $m-1$ patches, for
$2\leq m\leq M$. Lastly, we perform k-means clustering on the full estimate of
top singular vectors. The full cluster-quilting algorithm is summarized in
Algorithm 1. We note here that the blockwise SVD and consecutive block
matching ideas also appeared in an earlier work (Bishop and Yu, , 2014) for
PSD matrix completion from block observations. However, we consider a
different model, i.e. a mixture model, with a different objective, i.e.
clustering, and hence our matching step looks for the best linear
transformation between patches through least squares instead of finding the
best rotation matrix as in Bishop and Yu,  (2014).

Input: $\\{X_{T_{m},I_{m}}\in\mathbb{R}^{|T_{m}|\times|I_{m}|},m\in 1,\ldots
M\\}$, $r>0$

Initialize: $\widetilde{H}=\boldsymbol{0}_{p\times r}$,
$\widetilde{V}=\boldsymbol{0}_{n\times r}$.

Find: low-rank solution for first patch:

(1) Let $\widehat{U}_{1}\widehat{\Lambda}_{1}\widehat{V}_{1}^{\top}$ be the
top-$r$ SVD of $X_{T_{1},I_{1}}$

(2) Set $\widetilde{V}_{I_{1},\,:}=\widehat{V}_{1}$,
$\widetilde{H}_{T_{1},\,:}=\widehat{U}_{1}\widehat{\Lambda}_{1}$

for  _$m=2,\ldots M$_ do

(1) Find low-rank solution for the $m$-th patch:

(a) Let $\widehat{U}_{m}\widehat{\Lambda}_{m}\widehat{V}_{m}^{\top}$ be the
top-$r$ SVD of $X_{T_{m},I_{m}}$

(c) Let $\widehat{H}_{m}=\widehat{U}_{m}\widehat{\Lambda}_{m}$.

(2) Merge with previous patches:

(a) Find overlaps $J_{m}^{(1)}=(\bigcup_{j=1}^{m-1}I_{j})\cap
I_{m},\,J_{m}^{(2)}=\\{j:I_{m}[j]\in J^{(1)}\\}$.

(b) Find the best linear transform
$G_{m}=((\widehat{V}_{m})_{J_{m}^{(2)},:}^{\top}(\widehat{V}_{m})_{J_{m}^{(2)},:})^{\dagger}(\widehat{V}_{m})_{J_{m}^{(2)},:}^{\top}\widetilde{V}_{J_{m}^{(1)},:}$;

(c) Set $\widetilde{V}_{I_{m}\backslash
J_{m}^{(1)},\,:}=(\widehat{V}_{m})_{\backslash J_{m}^{(2)},\,:}G_{m}$,
$\widetilde{H}_{T_{m}}=\widehat{H}_{m}(G_{m}^{\top})^{-1}$.

end for

Post-processing of $\widetilde{H}\widetilde{V}^{\top}$: Let the SVD of
$\widetilde{H}\widetilde{V}^{\top}$ be
$\widehat{U}\widehat{\Lambda}\widehat{V}^{\top}$.

Perform $k$-means on the columns of $\widehat{\Lambda}\widehat{V}^{\top}$:
$(\widehat{z},\\{\widehat{c}_{j}\\}_{j=1}^{K})=\arg\min_{z\in[K]^{n},\\{\widehat{c}_{j}\\}_{j=1}^{K}\in\mathbb{R}^{r}}\sum_{i=1}^{n}\|(\widehat{\Lambda}\widehat{V}^{\top})_{:,\,i}-\widehat{c}_{z_{i}}\|_{2}^{2}$
return _Clustering labels $\widehat{z}\in[K]^{n}$, cluster means
$\\{\widehat{U}\widehat{c}_{j}\\}_{j=1}^{K}$_

Algorithm 1 Cluster Quilting

###  2.1 Practical Considerations

Input: $\\{X_{T_{m},I_{m}}\in\mathbb{R}^{|T_{m}|\times|I_{m}|},m\in 1,\ldots
M\\}$, score function $s(\cdot,\cdot):[M]\times 2^{[n]}\rightarrow[0,\infty)$

Initialize: Let $\hat{\pi}(1)\neq\hat{\pi}(2)\in[M]$ be an ordered pair of
patch indices such that
$(\hat{\pi}(1),\,\hat{\pi}(2))=\operatorname*{arg\,max}_{k_{1}\neq
k_{2}\in[M]}s(k_{2},I_{k_{1}})$. for  _$3\leq k\leq M$_ do

  1. 1.

Let $L=\\{\hat{\pi}(1),\,\dots,\,\hat{\pi}(k-1)\\}$.

  2. 2.

Let $\hat{\pi}(k)=\operatorname*{arg\,max}_{j:j\notin L}s(j,\cup_{l\in
L}I_{l})$.

end for

return _Ordering function $\hat{\pi}(\cdot):[M]\rightarrow[M]$._

Algorithm 2 Forward Search for Block Ordering

The application of the Cluster Quilting algorithm involves several practical
choices. Firstly, the ordering of the patches is critical, as for the first
step of Algorithm 1 to be applicable, we need to ensure that
$I_{m}\cap(\cup_{l<m}I_{l})\neq\emptyset$; (ii) for accurate matching in each
step (computation of $G_{m}$), the overlapping signal within
$J^{(1)}_{m}=I_{m}\cap(\cup_{l<m}I_{l})$ also needs to be sufficient.
Motivated by this, we propose to first order the patches appropriately with
the following ordering function $\hat{\pi}:[M]\rightarrow[M]$:

$\hat{\pi}=\operatorname*{arg\,max}_{\pi:[M]\rightarrow[M]}\prod_{m=2}^{M}s(\pi(m),\cup_{l<m}I_{\pi(l)}),$
(2)

where $s(\pi(m),\cup_{l<m}I_{\pi(l)})$ is a generic score function that
measures the overlapping signal between the $\pi(m)$’s patch and prior patches
$\cup_{l<m}I_{\pi(l)}$ in the sample space, which can also be interpreted as
the connection strength in the patch graph $\mathcal{G}$. Our choice to
maximize the multiplication of all overlapping scores is motivated by our
theoretical results in Section 3, and this minimizes the error propagation of
sequential patch matching in Alg. 1. The generic score functions can also take
different forms: one simple idea is to take $s(k,I)=|I_{k}\cap I|$ as the
overlapping size; we will also introduce a theory-motivated score function in
Section 3 and show theoretical guarantees with this data-driven patch
ordering. Given a score function $s(k,I)$, one can solve (2) by searching over
the whole space of feasible orderings such that
$I_{m}\cap(\cup_{l<m}I_{l})\neq\emptyset$. Since the number of patches $M$ is
often not large (under $10$) in our motivating applications, e.g., multiomics
data sets or calcium imaging data, this combinatorial optimization problem is
still computationally feasible. Given $\hat{\pi}(\cdot)$, we can then apply
Alg. 1 on $\\{X_{T_{\hat{\pi}(m)},I_{\hat{\pi}(m)}}\\}_{m=1}^{M}$. For
computational speed up when $M$ is large, we also propose a greedy forward
search algorithm for the patch-ordering, summarized in Alg. 2. In this schema,
Alg. 2 first chooses the $\hat{\pi}(1),\,\hat{\pi}(2)$ such that
$s(\hat{\pi}(2),I_{\hat{\pi}(1)})$ is maximized over all ordered patch pairs;
then it sequentially selects $\hat{\pi}(k)$ that maximizes
$s(\hat{\pi}(k),\cup_{l<k}I_{\hat{\pi}(l)})$.

The Cluster Quilting algorithm also requires the selection of both the number
of clusters as well as the rank of the spectral clustering step. To select
these hyperparameters, we follow the previously proposed prediction validation
approach (Tibshirani and Walther, , 2005) across a grid of potential
hyperparameter values selected a priori. For this approach, the data are split
into a training and a test set and clustered separately using each potential
combination of rank and number of clusters. A random forest model is then
trained on the former to predict the cluster assignments of the observations
in the latter. The final rank and number of clusters are chosen by selecting
the setting which provides the highest agreement between the labels found by
the clustering algorithm and the predictions from the supervised model.

##  3 Theoretical Guarantees

In this section, we provide theoretical characterization for the clustering
performance of our Alg. 1. We first assume that Alg. 1 is applied on the
original data patch ordering $1,\,\dots,\,M$, without making probabilistic
assumptions on the noise distribution; Section 3.1 gives the consequences
under sub-Gaussian noise; Section 3.2 then provides theory for the data-driven
selection of the patch ordering.

Here, we first introduce some key technical quantities we will use throughout
Section 3. Define the misclustering rate as:
$\ell(\widehat{z},z):=\min_{\phi\in\Phi}\frac{1}{n}\sum_{i\in[n]}{{\mathbbm{1}}}\\{\phi(\widehat{z}_{i})\neq
z_{i}\\},$ where $\Phi=\\{\phi:\phi$ is a bijection from $[K]$ to $[K]\\}$. We
will show upper bounds for this misclustering rate $\ell(\widehat{z},z)$ that
hold with high probability. For $1\leq j\leq K$, let
$n_{j}=\sum_{i=1}^{n}{{\mathbbm{1}}}\\{z_{i}=j\\}$ be the sample size in each
cluster $j$. Let $\Delta=\min_{j\neq
j^{\prime}}\|\theta_{j}-\theta_{j^{\prime}}\|_{2}$ be the minimum separation
distance between clusters, and $\Delta_{m}=\min_{j\neq
j^{\prime}}\|(\theta_{j})_{T_{m}}-(\theta_{j^{\prime}})_{T_{m}}\|_{2}$ be the
separation distance within block $m$. Given a patch ordering
$\pi:[M]\rightarrow[M]$, for $m\geq 2$, let
$J_{m}^{(1)}(\pi)=(\cup_{j=1}^{m-1}I_{\pi(j)})\cap I_{\pi(m)}$ be the overlap
of the $m$th block with prior blocks. We use $0<\gamma_{m}(\pi)<1$ to denote
the signal proportion of the overlapping part $J_{m}^{(1)}(\pi)$ in the whole
block $I_{\pi(m)}$:
$\gamma_{m}(\pi)=\frac{\sigma_{r}(X^{*}_{T_{\pi(m)},\,J_{m}^{(1)}(\pi)})}{\|X_{T_{\pi(m)},\,I_{\pi(m)}}^{*}\|}$,
and $\gamma(\pi)=\min_{2\leq m\leq K}\gamma_{m}(\pi)$. When the original patch
ordering is in use, e.g., $\pi$ is the identity mapping, we abbreviate these
notations to $J_{m}^{(1)}$, $\gamma_{m}$, and $\gamma$. Let $\kappa_{m}$ be
the condition number of the $m$th block. Consider the population block SVD
$X_{T_{m},\,I_{m}}^{*}=U_{m}^{*}\Lambda_{m}^{*}V_{m}^{*\top}$, and define the
matching matrix for block $m$ as
$B_{m}^{*}=\Lambda^{*}U_{T_{m},:}^{*}U_{m}^{*}\Lambda_{m}^{*-1}$, which
transforms the whole matrix SVD into the block SVD results:
$V_{m}^{*}=V_{I_{m},\,:}^{*}B_{m}^{*}$. We also define $\beta_{\max}$ and
$\beta_{\min}$ to quantify the heterogeneity across blocks:
$\beta_{\max}=\max_{1\leq m,\,m^{\prime}\leq
M}\|B_{m}^{*-1}B_{m^{\prime}}^{*}\|$, $\beta_{\min}=\min_{1\leq
m,\,m^{\prime}\leq M}\sigma_{r}(B_{m}^{*-1}B_{m^{\prime}}^{*})$. Another
measure of the heterogeneity of signals across blocks is defined as
$\alpha_{m,\,m^{\prime}}=\frac{\|X^{*}_{T_{m},I_{m}}\|}{\sigma_{r}(X_{T_{m^{\prime}},I_{m^{\prime}}}^{*})}$,
and $\alpha_{\max}=\max_{m,\,m^{\prime}}\alpha_{m,\,m^{\prime}}$. The
following regularity and block-wise SNR conditions are needed for our general
theory with the original patch ordering.

######  Assumption 1 (Regularity conditions on centroids).

The clustering centroids $\theta_{1},\dots,\theta_{K}\in\mathbb{R}^{p}$
satisfies $\max_{j}\|\theta_{j}\|_{2}\leq C\Delta=C\min_{j\neq
j^{\prime}}\|\theta_{j}-\theta_{j^{\prime}}\|_{2}$ for some universal constant
$C>0$

Assumption 1 is a regularity condition on clustering centroids, which also
appears in prior works on spectral clustering Abbe et al.,  (2022). It is a
mild condition as long as the data matrix is appropriately centered as a
preprocessing step.

######  Assumption 2 (SNR condition).

For some universal constant $C>0$, the SNR of each block $m$ satisfies

$\frac{\sigma_{r}(X^{*}_{T_{m},I_{m}})}{\|E_{T_{m},I_{m}}\|}\geq\frac{C\alpha_{\max}(\beta_{\max}+1)}{\beta_{\min}^{2}}\prod_{l=2}^{M}\left(\frac{1.1}{\gamma_{l}}+1\right)\sqrt{\frac{rK\max_{j}n_{j}}{\min_{j}n_{j}}}.$

Assumption 2 suggests that the clustering separation signal in each block well
exceeds the noise level. When the number of blocks is larger (i.e. larger
$M$), a stronger SNR condition is also required. This is because our Cluster
Quilting algorithm (Alg. 1) requires a matching step that sequentially
transforms the SVD results of each block to match its previous block. It takes
$M-1$ sequential transformations to obtain the full singular space estimate,
each of which leads to an additional error factor $\frac{1.1}{\gamma_{l}}+1$
depending on the signal ratio between each block and its overlapping part. The
sequential error propagation leads to the factor
$\prod_{l=2}^{M}\left(\frac{1.1}{\gamma_{l}}+1\right)$ in Assumption 2.

######  Theorem 1.

Consider the model set-up described in Section 2 and Alg. 1. Suppose that
Assumptions 1, 2 hold. Then the outputs of Alg. 1 satisfy

$\begin{split}\ell(\widehat{z},z)\leq&B^{2}(X^{*},\mathcal{O})\frac{rK\max_{j}n_{j}}{n}\max_{m}\frac{\|E_{T_{m},I_{m}}\|^{2}}{\sigma_{r}^{2}(X^{*}_{T_{m},I_{m}})},\\\
\min_{\phi\in\Phi}\max_{1\leq j\leq
K}\|\widehat{U}\widehat{c}_{j}-\theta_{\phi(j)}\|_{2}\leq&B(X^{*},\mathcal{O})\sqrt{\frac{rK\max_{j}n_{j}}{\min_{j}n_{j}}}\max_{m}\frac{\|E_{T_{m},I_{m}}\|}{\sigma_{r}(X^{*}_{T_{m},I_{m}})}\max_{j}\|\theta_{j}\|_{2}.\end{split}$
(3)

where
$B(X^{*},\mathcal{O})=\frac{C\alpha_{\max}(\beta_{\max}+1)}{\beta_{\min}^{2}}\prod_{l=2}^{M}\left(\frac{1.1}{\gamma_{l}}+1\right)$
depends on the ground truth matrix $X^{*}$ and the measurement pattern
$\mathcal{O}$.

Theorem 1 suggests that as long as the SNR within each block
$\frac{\sigma_{r}(X^{*}_{T_{m},I_{m}})}{\|E_{T_{m},I_{m}}\|}$ is sufficiently
large, accuracy of both clustering centroid estimation are well controlled by
a function of the SNR, rank $r$, number of clusters $K$, and number of blocks
$M$.

######  Remark 1 (Consistent clustering).

Recall that parameters $\alpha_{\max},\,\beta_{\max}\geq 1$, $\beta_{\min}\leq
1$ are defined in the beginning of Section 3, and they characterize the
heterogeneity of different blocks. When all blocks $X^{*}_{T_{m},I_{m}}$ are
well conditioned and the signal in different blocks are similar, we may safely
assume $\frac{\alpha_{\max}(\beta_{\max}+1)}{\beta_{\min}^{2}}$ to be bounded.
If $K\max_{j}n_{j}\leq Cn$ (different clusters are balanced) also holds, we
know that the misclustering rate $\ell(\widehat{z},z)=o(1)$ as long as
$\sigma_{r}(X^{*}_{T_{m},I_{m}})\gg\sqrt{r}\prod_{l=2}^{M}\left(\frac{1.1}{\gamma_{l}}+1\right)\|E_{T_{m},I_{m}}\|$.

The patchwise observational pattern plays a critical role in the difficulty of
the problem. As illustrated in Assumption 2, Theorem 1, when
$\prod_{l=2}^{M}\left(\frac{1.1}{\gamma_{l}}+1\right)$ is larger, the problem
gets harder as we need a stronger SNR condition and have a larger
misclustering rate. Therefore, we ideally hope for larger $\gamma_{l}$’s,
which are the fraction of the overlapping signal between the $l$th block with
previous blocks, as a proportion of the whole signal in the $l$th block:
$\gamma_{l}=\frac{\sigma_{r}(X^{*}_{T_{l},I_{l}\cap(\cup_{k<l}I_{k})})}{\|X^{*}_{T_{l},I_{l}}\|}$.
Different from prior works for PSD matrix or covariance completion (Bishop and
Yu, , 2014), our theory hinges on a much smaller merging factor
$(1.1\gamma_{l}^{-1}+1)$; we use spectral norm instead of Frobenious norm to
characterize noise level; and we require no condition on the gap between non-
zero singular values. This more refined result is based on a careful induction
proof, and it motivates a data-driven, theory-motivated ordering method for
the patches in Section 3.2.

In Theorem 1, the signal strength is captured by the minimum singular value of
the ground truth matrix within each block. In fact, we can also express the
signal strength in terms of the separation distances between clusters, which
might be more natural in the clustering context. In the following, we present
a corollary that illustrates the effect of seperation distances on the
clustering performance of Alg. 1.

######  Assumption 3 (Separation condition).

Within each block $m$, the minimum separation between clusters
$\Delta_{m}=\min_{j\neq
j^{\prime}}\|\theta_{T_{m},j}-\theta_{T_{m},j^{\prime}}\|_{2}$ satisfies

$\frac{\sqrt{\min_{j}n_{j,m}}\Delta_{m}}{\|E_{T_{m},I_{m}}\|}\geq\frac{C\alpha_{\max}(\beta_{\max}+1)\kappa_{r}(\Theta^{*}_{T_{m},:})}{\beta_{\min}^{2}}\prod_{l=2}^{M}\left(\frac{1.1}{\gamma_{l}}+1\right)\sqrt{\frac{rK\max_{j}n_{j}}{\min_{j}n_{j}}},$

for some universal constant $C>0$, where
$\kappa_{r}(\Theta^{*}_{T_{m},:})=\frac{\sigma_{1}(\Theta^{*}_{T_{m},:})}{\sigma_{r}(\Theta^{*}_{T_{m},:})}$,
the condition number of the $m$th block’s centroid matrix; $n_{j,m}=\sum_{i\in
I_{m}}{{\mathbbm{1}}}\\{z_{i}=j\\}$ is the number of samples within $m$th
block and $j$th cluster.

Here in Assumption 3, the separation SNR
$\frac{\sqrt{\min_{j}n_{j,m}}\Delta_{m}}{\|E_{T_{m},I_{m}}\|}$ can be
approximately viewed as the ratio between the $\ell_{2}$ separation distance
amongst clusters and the noise level for each sample within block $m$. In
particular, in the balanced setting ($\max_{1\leq j\leq K}n_{j,m}\leq
C\min_{1\leq j\leq K}n_{j,m}$),

$\frac{\|E_{T_{m},I_{m}}\|}{\sqrt{\min_{j}n_{j,m}}}\leq
C\frac{\|E_{T_{m},I_{m}}\|}{\sqrt{n_{m}}}\leq C\sqrt{\frac{1}{n_{m}}\sum_{i\in
I_{m}}\|E_{T_{m},i}\|_{2}^{2}},$

which is the average $\ell_{2}$ norm of the noise associated with each sample
within block $m$.

######  Corollary 1 (SNR in the form of separation distance).

Consider the model set-up described in Section 2 and Alg. 1. Suppose that
Assumptions 1 and 3 hold. Then the outputs of Alg. 1 satisfy

$\begin{split}\ell(\widehat{z},z)\leq&2B^{2}(X^{*},\mathcal{O})\frac{rK\max_{j}n_{j}}{n}\max_{m}\frac{\kappa_{r}^{2}(\Theta^{*}_{T_{m},:})\|E_{T_{m},I_{m}}\|^{2}}{\min_{1\leq
j\leq K}n_{j,m}\Delta_{m}^{2}},\end{split}$ (4)

Corollary 1 suggests that in a balanced ($\frac{K\max_{j}n_{j}}{n}\leq C$),
well-conditioned setting with bounded heterogeneity across blocks, as long as
the separation distance within each block satisfies
$\Delta_{m}\gg\sqrt{r}\prod_{l=2}^{M}\left(\frac{1.1}{\gamma_{l}}+1\right)$
$\frac{\|E_{T_{m},I_{m}}\|}{\min_{1\leq j\leq K}\sqrt{n_{j,m}}}$, our cluster
quilting algorithm achieves consistent clustering with misclustering rate
$\ell(\widehat{z},z)=o(1)$. We will provide more comparisons with prior
spectral clustering theory in Section 3.1, under sub-Gaussian noise
assumptions.

###  3.1 Special Case: Sub-Gaussian Noise and Dependent Features

Theorem 1 is a general deterministic result that holds for all noise types. To
make it more concrete and comparable with prior works, here we consider an
example of mean-zero sub-Gaussian noise. Specifically, suppose that there
exist $n^{\prime}\geq n,\,p^{\prime}\geq p$, matrices
$A_{r}\in\mathbb{R}^{n\times n^{\prime}}$ and $A_{c}\in\mathbb{R}^{p\times
p^{\prime}}$, such that

$E=A_{r}ZA_{c}^{\top},\quad Z\in\mathbb{R}^{n^{\prime}\times p^{\prime}},\quad
Z_{j,k}\overset{\text{ind.}}{\sim}\text{sub-Gaussian}(\sigma)\text{ with zero
mean}.$ (5)

A special case of (5) is when $n^{\prime}=n$, $p^{\prime}=p$, and $A_{r}$,
$A_{c}$ are identity matrices; this is equivalent to write
$E_{j,k}\overset{\text{ind.}}{\sim}\text{sub-Gaussian}(\sigma)$. By
considering general matrices $A_{r}$ and $A_{c}$, we also allow linear
dependency across features and samples.

In the following, we present the consequence of Theorem 1 under the noise
model (5); For presentation simplicity, we consider the following regularity
assumption.

######  Assumption 4 (Bounded patch heterogeneity, condition number, balanced
clusters).

$\begin{split}&\alpha_{\max},\beta_{\max},\beta_{\min}^{-1},\frac{\max_{1\leq
j\leq K}n_{j}}{\min_{1\leq j\leq K}n_{j}}\leq C,\\\ &\frac{\max_{1\leq j\leq
K}n_{j,m}}{\min_{1\leq j\leq K}n_{j,m}},\kappa_{r}(\Theta^{*}_{T_{m},:})\leq
C,\quad 1\leq m\leq M.\end{split}$ (6)

Here, $\alpha_{\max},\beta_{\max},\beta_{\min}^{-1}\geq 1$ quantify the
heterogeneity of the signal levels in different blocks;
$\frac{\max_{j}n_{j}}{\min_{j}n_{j}}$,
$\frac{\max_{j}n_{j,m}}{\max_{j}n_{j,m}}$ quantify how unbalanced that
different clusters are; $\kappa_{r}(\Theta^{*}_{T_{m},:})$ gives the condition
number of the centroid matrix within each block.

######  Theorem 2.

Consider the model set-up described in Section 2. Suppose that the noise $E$
satisfies (5) and Assumptions 1, 4 hold. If the minimum separation distance
within each block $1\leq m\leq M$ satisfies

$\Delta_{m}\geq
C\sigma\|A_{r}\|\|A_{c}\|\sqrt{r}K\prod_{l=2}^{M}\left(\frac{1.1}{\gamma_{l}}+1\right)\left(1+\sqrt{\frac{p_{m}}{n_{m}}}\right),$
(7)

then with probability at least $1-C\sum_{m=1}^{M}\exp\\{-c(n_{m}+p_{m})\\}$,
the misclustering rate is bounded as follows: $\ell(\widehat{z},z)\leq
C\sigma^{2}\|A_{r}\|^{2}\|A_{c}\|^{2}rK\prod_{l=2}^{M}\left(\frac{1.1}{\gamma_{l}}+1\right)^{2}\max_{m}\left(1+\frac{p_{m}}{n_{m}}\right)\Delta_{m}^{-2}.$

Theorem 2 suggests that we can obtain consistent clustering as long as

$\Delta_{m}\gg\sigma\|A_{r}\|\|A_{c}\|\sqrt{r}K^{\frac{3}{2}}\prod_{l=2}^{M}\left(\frac{1.1}{\gamma_{l}}+1\right)\left(1+\sqrt{\frac{p_{m}}{n_{m}}}\right),$

where $\sigma\|A_{r}\|\|A_{c}\|$ represents the noise level;
$\prod_{l=2}^{M}\left(\frac{1.1}{\gamma_{l}}+1\right)$ quantifies the effect
of matching $M$ semi-overlapping blocks;
$\left(1+\sqrt{\frac{p_{m}}{n_{m}}}\right)$ reflects the effect of dimension
and sample size within each block. Similar to our discussion for Assumption 2,
the exponential dependence on $M$ is due to the nature of the observational
pattern: the first and last blocks are only connected through the $M-1$
sequential overlaps between consecutive blocks.

######  Remark 2 (Comparison with prior works).

If the noise matrix $E$ has independent sub-Gaussian-$\sigma$ entries
($A_{r}=I_{n\times n},\,A_{c}=I_{p\times p}$) in our setting, then the minimum
separation condition for consistent clustering ($\ell(\widehat{z},z)=o(1)$)
would be
$\Delta_{m}\gg\sigma\sqrt{r}K\prod_{l=2}^{M}\left(\frac{1.1}{\gamma_{l}}+1\right)\left(1+\sqrt{\frac{p_{m}}{n_{m}}}\right)$.
As a comparison, in prior literature (Löffler et al., , 2021) on spectral
clustering for Gaussian mixtures from the full data matrix, they require
$\Delta\gg\sigma\sqrt{r}(1+\sqrt{\frac{p}{n}})$. The additional factor we
require $K\prod_{l=2}^{M}\left(\frac{1.1}{\gamma_{l}}+1\right)$ is due to the
cost of matching $M$ semi-overlapping blocks. One might note that our explicit
bound for the misclustering rate only scales as $\Delta^{-2}$ instead of the
exponential rate $\exp\\{-c\Delta\\}$ as in Löffler et al.,  (2021). This is
due to that our blockwise SVD and sequential matching steps in Alg. 1 creates
complicated dependency structures between our algorithm output and all
samples, and hence the proof techniques for deriving exponential error bounds
are not applicable.

One key quantity in Theorem 2 is $\frac{p_{m}}{n}$, the ratio between number
of features in each block and the sample size. In the high-dimensional setting
where $\frac{p_{m}}{n_{m}}$ tends to infinity, Theorem 2 can still be useful
as long as the coordinate-wise separation distances are lower bounded by a
constant; in this case, the separation distance $\Delta_{m}=\min_{k_{1}\neq
k_{2}}\|(\theta^{*}_{k_{1}})_{T_{m}}-(\theta^{*}_{k_{2}})_{T_{m}}\|_{2}$ will
also grow with the rate $\sqrt{p_{m}}$. On the other hand, when the centroid
difference is sparse, a different clustering method that exploits such
sparsity might be more appropriate.

######  Remark 3 (Effect of dependency).

We note that the effect of dependency is reflected by $\|A_{c}\|$ and
$\|A_{r}\|$ in the separation condition and misclustering rate. Motivated by
the neuronal functional data where firing activities are collected over time,
let’s consider the autoregressive process as an example. Suppose that the
noise $E_{:,i}$ comes from an AR(1) process with coefficient $\rho\in(-1,1)$:
$E_{k,i}=\rho E_{k-1,i}+\varepsilon_{k,i}$, where $\varepsilon_{k,i}$ and
$E_{1,i}$ are independent $\sigma_{\varepsilon}$-sub-Gaussian, then we can let
$A_{c}$ be the identity matrix, and $A_{r}$ be a lower triangular matrix with
$(A_{r})_{i,j}=\begin{cases}1,&i=j;\\\ \rho,&i=j+1\end{cases}.$ In this
setting, $\|A_{c}\|=1$, $\|A_{r}\|\leq 1+|\rho|$; a larger autoregressive
coefficient $|\rho|$ leads to a stronger SNR condition, and a higher
misclustering rate.

###  3.2 Data-driven Selection of Patch Ordering

Recall that in our theoretical results in Section 3, the SNR conditions and
final misclustering error hinge on one key factor:
$\prod_{m=2}^{M}(1.1\gamma_{m}^{-1}(\pi)+1)$, where
$\gamma_{m}(\pi)=\frac{\sigma_{r}(X^{*}_{T_{\pi(m)},J_{m}^{(1)}(\pi)})}{\|X^{*}_{T_{\pi(m)},I_{\pi(m)}}\|}$
is the proportion of the overlapping signals in the $m$th patch. Let
$\pi^{*}=\operatorname*{arg\,min}_{\pi:[M]\rightarrow[M]}\prod_{m=2}^{M}(1.1\gamma_{m}^{-1}(\pi)+1)$
be the minimizer of this factor, i.e., the oracle patch ordering. Since each
$\gamma_{m}(\pi)$ depends on the unknown ground truth matrix $X^{*}$, we
propose to estimate it by its empirical counterpart
$\hat{\gamma}_{m}(\pi)=\frac{\sigma_{r}(X_{T_{\pi(m)},J_{m}^{(1)}(\pi)})}{\|X_{T_{\pi(m)},I_{\pi(m)}}\|}$;
and then find the ordering $\pi$ that minimizes
$\prod_{m=2}^{M}(1.1\hat{\gamma}_{m}^{-1}(\pi)+1)$. This is equivalent to
setting the score function in (2) as follows:

$s(k,I)=\left(\frac{1.1\|X_{T_{k},I_{k}}\|}{\sigma_{r}(X_{T_{k},I_{k}\cap
I})}+1\right)^{-1}.$ (8)

The following theorem gives the clustering guarantee with this data-driven
selection of the patch ordering.

######  Theorem 3.

Suppose that the factor $\prod_{m=2}^{M}(1.1\gamma_{m}^{-1}+1)$ is substituted
by $e^{2}\prod_{m=2}^{M}(1.1\gamma_{m}^{-1}(\pi^{*})+1)$ in Assumptions 2, 3,
Theorems 1, 2, and Corollary 1. Then Theorems 1, 2, and Corollary 1 still hold
for Alg. 1 applied on $\\{X_{T_{\hat{\pi}(m)},I_{\hat{\pi}(m)}}\\}_{m=1}^{M}$,
where the patch ordering $\hat{\pi}$ satisfies (2) with the score function
(8).

Theorem 3 suggests that the data-driven selection of the patch ordering leads
to almost the same performance as the oracle ordering; the slight constant
factor inflation ($e^{2}$) is due to the finite sample approximation of the
score function. The proof of Theorem 3 is included in Section A.4 of the
Supplementary material, and builds on careful analysis of the accumulated
finite sample errors across $M$ patches.

##  4 Simulation Studies

In this section, we study the performance of the proposed Cluster Quilting
method in the mosaic patch data observation contexts using data simulated from
a Gaussian mixture model. We compare our Cluster Quilting approach to several
of the incomplete spectral clustering methods: IMSC-AGL (Wen et al., , 2018),
IMG (Zhao et al., , 2016), DAIMC (Hu and Chen, 2019a, ), and OPIMC (Hu and
Chen, 2019b, ). We also compare our procedure to the application of standard
spectral clustering after imputing the missing portion of the data matrix
using nuclear-norm regularized least squares (Chang et al., , 2022). The
accuracy of each of the methods is quantified using the adjusted Rand index
(Santos and Embrechts, , 2009) with respect to true underlying cluster labels
in the generative models. In all simulation studies, we show results for
oracle hyperparameter tuning. We tune other hyperparameters specific to each
comparison method using a grid search and selection by best adjusted Rand
Index score.

###  4.1 Data Generation Process

For the simulation studies below, we generate the underlying full data from
Gaussian mixture models with low-rank cluster centroid matrices. We first
create the true cluster labels by randomly grouping $n$ samples into $K$
clusters, encoded by matrix $F^{*}\in\mathbb{R}^{n\times K}$ where

$F^{*}_{ik}=\begin{cases}1,&\text{sample }i\text{ belongs to cluster }k;\\\
0,&\text{otherwise}.\end{cases}$

To construct a low-rank cluster centroid matrix, we generate two rank $r$
matrices with $r<K$: a random orthogonal matrix
$\boldsymbol{Z}\in\mathbb{R}^{p\times r}$ and a matrix
$\boldsymbol{W}\in\mathbb{R}^{K\times r}$ whose entries are randomly selected
from $\\{-d,0,d\\}$ for a constant $d$. For the latter matrix, we check that
the rows are not linearly dependent and re-simulate new entries if this
requirement is not met. The cluster centroid matrix
$\Theta^{*}=(\theta_{1},\cdots,\theta_{K})^{\top}\in\mathbb{R}^{K\times p}$ is
formulated as the product of these two rank $r$ matrices, i.e.
$\Theta^{*}=\boldsymbol{W}\boldsymbol{Z}^{\top}$. This generation ensures that
$\Theta^{*}$ has rank $r$, with expected average separation distance
$\mathbb{E}\|(\Theta^{*})_{i,:}-(\Theta^{*})_{j,:}\|_{2}^{2}=\frac{4}{3}rd^{2}$.
Data are then simulated from a multivariate Gaussian distribution with mean
matrix $B=F^{*}\Theta^{*}\in\mathbb{R}^{n\times p}$ and covariance matrix
$I_{p}$.

###  4.2 Simulation Study

We now analyze the effectiveness of the Cluster Quilting procedure in the
mosaic patch observation setting, where multi-view data is simulated. For this
simulation study, instead of creating a single cluster centroid matrix as we
do for the simulations in Section B.2, we generate separate cluster centroid
matrices for each view in the data, applying the procedure outlined in Section
4.1 for each view. To create the missing observation pattern, we first
randomly segment the rows of the full data matrix into $M$ blocks of
observations. Following this, for each observation block, we randomly select
$h$ views that are fully included in the set of observed entries while masking
all other columns for the corresponding rows in the observation block. Cluster
Quilting and the other incomplete spectral clustering methods are then applied
to the masked data and compared based on how well they recover the true
cluster labels. For this simulation study, we use default parameters of 840
observations partitioned into 3 underlying clusters, and 12 views of 50
features each for a total of 600 columns in the full data matrix. The cluster
centroid matrix for each view is generated as a rank 2 matrix with $d=7.5$.
After data generation, we randomly divide the rows into 4 blocks of equal
size. To create the mosaic patch observation structure, for each observation
block, we randomly select 6 views to be part of the set of observed entries.
We ensure in this step that the set of views selected for each observation
block are unique, that each feature is part of the set of observed entries at
least once, and that the graph of connected observation blocks with respect to
overlapping views can be fully traversed. We analyze the performance of
Clustering Quilting and other incomplete spectral clustering methods when
varying the number of unmasked views per observation block, the number of
observation blocks, the number of total features in the data, the number of
underlying clusters, the distance between clusters, and the correlation
between features.

(a) Views observed per block.

(b) Number of blocks.

(c) Number of total features.

(d) Number of true clusters.

(e) Cluster distance.

(f) Feature autocorrelation.

Figure 2: Performance of Cluster Quilting and comparison cluster imputation
methods with oracle tuning on Gaussian mixture model data in the mosaic patch
setting for different simulation parameters, measured by adjusted Rand index.

We show the results of the mosaic patch observation simulation study in Figure
2. The Cluster Quilting method outperforms the comparative incomplete spectral
clustering methods in the vast majority of the simulation settings; we see
that Cluster Quilting performs particularly better when the number of views
seen per block (Figure 2(a)) or the number of features (Figure 2(c)) is
relatively small. Additionally, the performance of the Cluster Quilting is
better when the number of observation blocks is smaller but tends to degrade
more than other methods with an increase in the number of blocks (Figure 2(b)
), which matches with the theoretical results of Section 3. We see that an
increase in the number of true underlying clusters tends to decrease the
performance of all methods consistently except for the IMSC-AGL algorithm
(Figure 2(d)), while the performance of methods is relatively consistent with
changing the distance between clusters (Figure 2(e)). Also, in Figure 2(f), we
find that increasing the correlation between features increases the accuracy
of cluster estimates.

In the Supporting Information, we include several additional simulation
studies. We first show results for a parallel empirical study as above in the
sequential patch observation setting. For both mosaic and sequential patch
observations, show results for data-driven hyperparameter tuning, for Gaussian
mixture model data in the high-dimensional case, as well as when the data are
generated from non-Gaussian copula mixture models. The results of these
additional studies paint a similar picture, i.e. that the Cluster Quilting
method performs similarly well, relative to the comparative incomplete
spectral clustering methods, for recovering cluster labels.

##  5 Multiomics Case Study

Below, we investigate the performance of the Cluster Quilting method in the
mosaic patch setting for identifying cancer types from multi-omics data. In
this case studies, we take complete data sets and create synthetic patch
missingness in the data matrix. As in Section 4, we compare our Cluster
Quilting approach to several other incomplete spectral clustering methods,
namely IMSC-AGL (Wen et al., , 2018), IMG (Zhao et al., , 2016), DAIMC (Hu and
Chen, 2019a, ), and OPIMC (Hu and Chen, 2019b, ). The data analyzed below come
from The Cancer Genome Atlas (TCGA) (Weinstein, J.N. et. al., , 2013), which
contains genetic, transcriptomic, and proteomic data for a wide-ranging sample
of tumors. For this case study, we amalgamate tumor data from 3 different
primary sites, namely kidney, breast, and lung, and we consider features
collected on each tumor from 5 different views: DNA methylation, miRNA
expression, RNAseq gene expression, and reverse phase protein array. In
accordance with the procedures in Lock and Dunson,  (2013) and Wang and Allen,
(2021), we preprocess the data set before analysis by performing a square root
transformation of the DNA methylation features and a log transformation of the
miRNA features, as well as filtering out columns in the gene expression data
with standard deviation less than 1.5. Additionally, we scale the data column-
wise before creating the partially observed data set.

(a) Oracle tuning.

(b) Data driven tuning.

Figure 3: Adjusted Rand index of cluster estimates from each incomplete
spectral clustering methods with respect to true primary sites of tumors in
the TCGA dataset. Results are shown across varying number of observation
blocks and number of views unmasked per block.

From the complete data set, we create synthetic missing data with a mosaic
patch observation pattern by randomly partitioning the rows of the data into
distinct observation blocks, then selecting a subset of views for each block
to include in the set of observed entries. We then apply Cluster Quilting and
the comparison incomplete spectral clustering methods to the synthetic
partially observed data set and measure the performance of each method using
adjusted Rand index with respect to the true primary site labels. The efficacy
of the Cluster Quilting and other incomplete spectral clustering methods are
shown across varying numbers of observation blocks and views seen per block;
for each combination of the two, we perform 50 replications and report the
mean and standard deviation of the adjusted Rand index for each method.
Additionally, we present results for both oracle and data-driven tuning of the
number of clusters and rank of the data matrix; for the latter, we use the
prediction validation method proposed in Tibshirani and Walther,  (2005).
Here, we set the oracle rank to be the rank selected by data-driven tuning on
spectral clustering when applied to the fully observed data.

The results of our study are presented in Figure 3, with oracle hyperparameter
tuning shown in Figure 3(a) and data-driven tuning in Figure 3(b). Across the
different numbers of views observed per patch, we observe that the Cluster
Quilting algorithm does best when the number of observation blocks is
relatively small, while several methods, namely the likelihood imputation and
IMSC-AGL methods, outperform Cluster Quilting, when the number of observation
patches in the data is larger. The degredation in accuracy of the Cluster
Quilting method aligns with the theoretical results presented in Section 3.
However, even when the number of blocks is larger, the Cluster Quilting method
still achieves close to the highest accuracy out of all methods.

Figure 4: Cluster heatmap for each individual view for the (a) true data
matrix and (b) imputed low-rank data matrix estimated from Cluster Quilting.
Cells are colored by scaled value of the matrix. Rows are ordered by true
primary site, with different underlying groupings partitioned; color labels on
the left bar show the estimated cluster labels returned by the Cluster
Quilting procedure. Columns are ordered by similarities from hierarchical
clustering.

We further investigate the cluster label estimates from the Cluster Quilting
procedure in Figure 4. Here, we show cluster heatmaps for each view included
in our TCGA empirical study, with the heatmap of the fully observed data set
in column (a) and the heatmap of the data matrix as imputed via the Cluster
Quilting algorithm ($\widehat{U}\widehat{\Lambda}\widehat{V}^{\top}$ in
Algorithm 1) for the corresponding views in column (b). The rows in the
heatmap are sorted by the estimated cluster labels returned by Cluster
Quilting, while the true primary sites of each tumor are shown in the colored
bar on the left side of the figure. The Cluster Quilting method is able to
correctly cluster together all of the tumors with kidney as the primary site,
but misclassifies several of the tumors from the breast and lung primary sites
as the other. This is reflected in the cluster heatmaps for both the true data
matrix and the imputed data matrix estimated via the Cluster Quilting
algorithm, in which we see that the transcriptomic and genetic information in
the kidney primary site tumors are more easily distinguishable compared those
from the other two primary sites. The misclassification between the lung and
breast cancer tumors can likely be explained by the relative similarity
between the two, which has been shown in previous literature (Gemma et al., ,
2001; Schwartz et al., , 1999); in particular, metastasis between the lung and
breast have been studied in great depth (Montel et al., , 2005; Jin et al., ,
2018; Landemaine et al., , 2008). In the Supplementary Material, we perform a
real-world data study in the sequential patch observation setting, in which we
analyze functional neuronal clustering from calcium imaging data.

##  6 Discussion

In this work, we study the problem of clustering in the patchwork learning
setting. We propose the Cluster Quilting procedure, the first spectral
clustering method on patchwork data with statistical guarantees. We have shown
finite sample misclustering rates based on the separation distance between
clusters and signal-to-noise ratio in the data patches, and we have
demonstrated through empirical studies that the Cluster Quilting method
outperforms other incomplete spectral clustering methods for recovering true
underlying cluster labels in both simulated and real-world data. While other
works in the literature have studied the general problem of cluster estimation
under the presence of missing data, including many techniques proposed for
partial multi-view spectral clustering, we are the first to propose a method
and study its performance from a statistical perspective for the patchwork
learning setting. Additionally, our work significantly differs from while also
complements the matrix completion literature with deterministic missing
patterns: prior works in this area mostly consider a more specific
observational pattern, e.g., causal panel data, sequential observational
blocks, etc., and have no clear implications in the clustering context.

While the Cluster Quilting method has shown promising results for clustering
in the patchwork learning setting, there are several open problems and
potential extensions which require further investigation. Our current
empirical studies have focused on the case where each cluster is well-
represented in each observation patch and where the amount of overlap between
patches is relatively evenly distributed and sufficiently large. Violations of
these assumptions, however, may occur in real-world applications, and it is
unclear how well Cluster Quilting will perform in these situations. We have
also restricted our work to the case where patch merging is done sequentially,
i.e. where we only consider iteratively merging a new patch to ones that have
already been previously consolidated. This does not have to be the case
algorithmically, and in fact it may be beneficial to merge non-sequentially in
order to maximize the total overlapping signal; however, this scheme presents
extra challenges for quantifying theoretical performance. From a broader
perspective, ideas from the Cluster Quilting algorithm could potentially be
extended for clustering patchwork-observed data in the context of multimodal
data integration; however, further work will be required to study whether our
method can be used to meet challenges specific to this setting, such as mixed
data types or for highly erosely measured modalities. Beyond data clustering,
our methods and ideas may also be leveraged for other problems in patchwork
learning, such as PCA, ranking, and community detection in networks.

## Acknowledgements

The authors gratefully acknowledge support from NSF NeuroNex-1707400, NIH
1R01GM140468, and NSF DMS-2210837.

##  Appendix A Proofs

In this section, we present the complete proofs for all our theoretical
results.

###### Proof of Theorem 1.

We first present a lemma that shows that spectral error bound for $X^{*}$
implies clustering consistency:

######  Lemma 1.

Suppose that $\|\widetilde{H}\widetilde{V}^{\top}-X^{*}\|\leq\varepsilon$. Let
$\Phi$ denote the set of all permutations on $[n]$. Then as long as
$\varepsilon<\sqrt{\frac{\min_{1\leq j\leq K}n_{j}}{64r}}\Delta$, we have

$\ell(\widehat{z},z)\leq\frac{32r\varepsilon^{2}}{n\Delta^{2}},\quad\min_{\phi\in\Phi}\max_{1\leq
j\leq
K}\|\widehat{U}\widehat{c}_{j}-\theta^{*}_{\phi(j)}\|_{2}^{2}\leq\frac{16r\varepsilon^{2}}{\min_{1\leq
j\leq K}n_{j}}.$

The proof of Lemma 1 follows similar ideas in Löffler et al.,  (2021). The
main difficulty of the full proof of Theorem 1 lies in showing the spectral
error bound for the output of Algorithm 1 so that we can apply Lemma 1. In the
following, we present a theorem that bounds
$\|\widetilde{H}\widetilde{V}^{\top}-X^{*}\|$ with high probability:

######  Theorem 1.

Consider the model set-up described in Section 2 and Algorithm 1. If
Assumption 2 holds, we have

$\|\widetilde{H}\widetilde{V}-X^{*}\|\leq\frac{C\alpha_{\max}(\beta_{\max}+1)}{\beta_{\min}^{2}}\prod_{l=2}^{M}\left(\frac{1.1}{\gamma_{l}}+1\right)\max_{m}\frac{\|E_{T_{m},I_{m}}\|}{\sigma_{r}(X^{*}_{T_{m},I_{m}})}\|X^{*}\|.$

Combining Theorem 1 and Assumption 2, we have

$\|\widetilde{H}\widetilde{V}-X^{*}\|\leq
C\|X^{*}\|\sqrt{\frac{\min_{j}n_{j}}{rK\max_{j}n_{j}}}.$

By Assumption 1,
$\|X^{*}\|\leq\|X^{*}\|_{F}\leq\max_{j}\|\theta_{j}\|_{2}\sqrt{n}\leq
C\Delta\sqrt{K\max_{j}n_{j}}$, and hence by choosing the constant $C$
appropriately in Assumptions 2, we have
$\|\widetilde{H}\widetilde{V}-X^{*}\|<\sqrt{\frac{\min_{j}n_{j}}{64r}}\Delta$.
Now we can apply Lemma 1, plug in the error bound in Theorem 1 into Lemma 1 to
complete the proof of Theorem 1. ∎

The following lemma for bounding the sum of sequential products will be useful
in many steps of our main proof.

######  Lemma 2.

For any sequence $\\{a_{1},a_{2},a_{n}\\}$, if $a_{l}\geq a_{\min}$ for all
$1\leq l\leq n$ and some $a_{\min}>1$, then we have

$\sum_{i=1}^{n}\prod_{l=1}^{i}a_{l}\leq\frac{\prod_{l=1}^{n}a_{l}a_{\min}-a_{1}}{a_{\min}-1}.$

###### Proof of Lemma 2.

Let $B=\sum_{i=1}^{n}\prod_{l=1}^{i}a_{l}$. We can then write

$\begin{split}Ba_{\min}=&\sum_{i=1}^{n}\prod_{l=1}^{i}a_{l}a_{\min}\\\
\leq&\sum_{i=1}^{n-1}\prod_{l=1}^{i+1}a_{l}+\prod_{l=1}^{n}a_{l}a_{\min}\\\
=&\sum_{i=1}^{n}\prod_{l=1}^{i}a_{l}+\prod_{l=1}^{n}a_{l}a_{\min}-a_{1}\\\
=&B+\prod_{l=1}^{n}a_{l}a_{\min}-a_{1}.\end{split}$

Therefore, $B(a_{\min}-1)\leq\prod_{l=1}^{n}a_{l}a_{\min}-a_{1}$. The proof is
now complete. ∎

###  A.1 Proof of Theorem 1

The proof consists of three main steps: (i) we first upper bound the spectral
errors for the SVD of each block; (ii) then we show that the matching matrices
for consecutive blocks are accurately estimated, given that each block has
small perturbation bound shown in the first step; (iii) combining the first
two results together, we obtain final error bounds for the full matrix factors
$\widetilde{V}\in\mathbb{R}^{n\times r}$ and
$\widetilde{H}\in\mathbb{R}^{p\times r}$.

#### Bounding spectral errors for each block

Consider the block-wise SVD for the population quantity
$X^{*}_{T_{m},I_{m}}=U_{m}^{*}\Lambda_{m}^{*}V_{m}^{*\top}$ with
$U_{m}^{*}\in\mathbb{O}^{p_{m}\times r}$, $\Lambda_{m}\in\mathbb{R}^{r\times
r}$, and $V_{m}\in\mathbb{O}^{n_{m}\times r}$. Then we would like to show
spectral error bounds for $\widehat{H}_{m}$ and $\widehat{V}_{m}$ as estimates
for $U_{m}^{*}\Lambda_{m}^{*}$ and $V_{m}^{*}$ up to rotation. By Weyl’s
inequality, $\sigma_{r+1}(X_{T_{m},I_{m}})\leq\|E_{T_{m},I_{m}}\|$. As long
the universal constant $C>0$ in Assumption 2 is taken to be larger than $2$,
we have
$\sigma_{r+1}(X_{T_{m},I_{m}})\leq\|E_{T_{m},I_{m}}\|\leq\frac{1}{2}\sigma_{r}(X^{*}_{T_{m},I_{m}})$.
Applying the Davis-Kahan’s theorem (see, e.g., Theorem 2.7 in Chen et al., ,
2021), we know that

$\begin{split}\min_{W\in\mathcal{O}^{r\times
r}}\|\widehat{V}_{m}W-V_{m}^{*}\|\leq&\frac{\sqrt{2}\|E_{T_{m},I_{m}}\|}{\sigma_{r}(X^{*}_{T_{m},I_{m}})-\sigma_{r+1}(X_{T_{m},I_{m}})}\\\
\leq&\frac{2\sqrt{2}\|E_{T_{m},I_{m}}\|}{\sigma_{r}(X^{*}_{T_{m},I_{m}})}.\end{split}$
(9)

For the rest of the proof, we use $\widehat{W}_{m}$ to denote the best
rotation matrix that achieves the minimum $\min_{W\in\mathcal{O}^{r\times
r}}\|\widehat{V}_{m}W-V_{m}^{*}\|$. Now we turn to upper bounding
$\widehat{H}_{m}\widehat{W}_{m}-U_{m}^{*}\Lambda_{m}^{*}$. Recognizing the
fact that
$\widehat{H}_{m}\widehat{W}_{m}=X_{T_{m},I_{m}}\widehat{V}_{m}\widehat{W}_{m}$
and $U_{m}^{*}\Lambda_{m}^{*}=X^{*}_{T_{m},I_{m}}V_{m}^{*}$, we can show the
following:

$\begin{split}\|\widehat{H}_{m}\widehat{W}_{m}-U_{m}^{*}\Lambda_{m}^{*}\|=&\|X_{T_{m},I_{m}}\widehat{V}_{m}\widehat{W}_{m}-X^{*}_{T_{m},I_{m}}V_{m}^{*}\|\\\
\leq&\|X^{*}_{T_{m},I_{m}}\|\|\widehat{V}_{m}\widehat{W}_{m}-V_{m}^{*}\|+\|E_{T_{m},I_{m}}\|\\\
\leq&4\kappa_{m}\|E_{T_{m},I_{m}}\|\end{split}$ (10)

where the second line utilizes (9). We have also used
$\kappa_{m}=\frac{\sigma_{1}(X^{*}_{T_{m},I_{m}})}{\sigma_{r}(X^{*}_{T_{m},I_{m}})}$
to denote the condition number of the ground truth matrix in block $m$.

#### Oracle matching matrices for consecutive blocks

Now we would like to establish a correspondence between our block-wise SVD
results and the sub-matrices of the full SVD. By establishing such a
correspondence, we can also find the oracle matching matrices between each
block and its prior blocks. Consider the SVD for the full ground truth matrix
$X^{*}=U^{*}\Lambda^{*}V^{*\top}$. Then for each block $m$, note that
$U_{T_{m},:}^{*}\Lambda^{*}V_{I_{m},:}^{*\top}=U_{m}^{*}\Lambda_{m}^{*}V_{m}^{*\top}$.
Hence we define
$B_{m}^{*}=\Lambda^{*}U_{T_{m},:}^{*\top}U_{m}^{*}\Lambda_{m}^{*-1}\in\mathbb{R}^{r\times
r}$, and one can show that $B_{m}^{*}$ can transform the $m$th block of the
full SVD matrix $V^{*}_{I_{m},:}$ to $V_{m}^{*}$:

$V_{m}^{*}=V^{*}_{I_{m},:}B_{m}^{*}.$ (11)

We can also see that

$U_{m}^{*}\Lambda_{m}^{*}B_{m}^{*\top}=U_{m}^{*}U_{m}^{*\top}U_{T_{m},:}^{*}\Lambda^{*}=U_{T_{m},:}^{*}\Lambda^{*},$
(12)

since $U_{T_{m},:}^{*}$ lies within the space spanned by $U_{m}^{*}$. Let
$\widetilde{V}^{*}=V^{*}B_{1}^{*}\widehat{W}_{1}^{\top}$, and
$\widetilde{H}^{*}=U^{*}\Lambda^{*}(B_{1}^{*-1})^{\top}\widehat{W}_{1}^{\top}$,
then we also have another decomposition for the ground truth matrix
$X^{*}=\widetilde{H}^{*}\widetilde{V}^{*\top}$, and they satisfy:

$(\widetilde{V}^{*})_{I_{1},:}=V_{1}^{*}\widehat{W}_{1}^{\top},\quad(\widetilde{H})^{*}_{T_{1},:}=U_{1}^{*}\Lambda_{1}^{*}\widehat{W}_{1}^{\top}.$

Therefore, our perturbation bounds for blockwise SVD in the first step already
suggest $\widetilde{V}_{I_{1},:}=\widehat{V}_{1}$ and
$\widetilde{H}_{T_{1},:}=\widehat{H}_{1}$ as good estimates for
$(\widetilde{V}^{*})_{I_{1},:}$ and $(\widetilde{H})^{*}_{T_{1},:}$. On the
other hand, for block $m>1$, some transformation is needed to match the first
block so that they can estimate $(\widetilde{V}^{*})_{I_{m},:}$ and
$(\widetilde{H}^{*})_{T_{m},:}$ well. Define

$G_{m}^{*}:=\widehat{W}_{m}B_{m}^{*-1}B_{1}^{*}\widehat{W}_{1}^{\top},$ (13)

then some computations show that

$(\widetilde{V}^{*})_{I_{m},:}=V_{m}^{*}\widehat{W}_{m}^{\top}G_{m}^{*},\quad(\widetilde{H})^{*}_{T_{m},:}=U_{m}^{*}\Lambda_{m}^{*}\widehat{W}_{m}^{\top}(G_{m}^{*})^{\top-1}.$

We will need to show that the oracle $G_{m}^{*}$ can be estimated well by
Algorithm 1 later in the proof, which then suggests all blocks of
$\widetilde{V}$ and $\widetilde{H}$ are accurate estimates of
$\widetilde{V}^{*}$ and $\widetilde{H}^{*}$.

Furthermore, the spectral error for the full low-rank reconstruction satisfies
the following:

$\begin{split}\|\widetilde{H}\widetilde{V}^{\top}-X^{*}\|=&\|\widetilde{H}\widetilde{V}^{\top}-\widetilde{H}^{*}\widetilde{V}^{*\top}\|\\\
\leq&\|\widetilde{H}\|\|\widetilde{V}-\widetilde{V}^{*}\|+\|\widetilde{H}-\widetilde{H}^{*}\|\|\widetilde{V}^{*}\|\\\
\leq&(\|\widetilde{H}^{*}\|+\|\widetilde{H}-\widetilde{H}^{*}\|)\|\widetilde{V}-\widetilde{V}^{*})\|+\|\widetilde{H}-\widetilde{H}^{*}\|\|\widetilde{V}^{*}\|.\end{split}$
(14)

For each block $m$, let $J_{m}^{(1)}=(\bigcup_{k=1}^{m-1}I_{k})\cap
I_{m},\,J_{m}^{(2)}=\\{j:I_{m}[j]\in J_{m}^{(1)}\\}$ be the overlapping set of
block $m$ with prior blocks. When $m=1$, $J_{m}^{(1)}$ and $J_{m}^{(2)}$ are
trivially defined as empty sets. By the definition of the spectral norm, we
know that the spectral error for both $\widetilde{V}$ and $\widetilde{H}$ can
be decomposed as follows:

$\begin{split}\|\widetilde{V}-\widetilde{V}^{*}\|\leq&\sqrt{\sum_{m=1}^{M}\|\widetilde{V}_{I_{m}\backslash
J_{m}^{(1)},:}-\widetilde{V}^{*}_{I_{m}\backslash J_{m}^{(1)},:}\|^{2}},\\\
\|\widetilde{H}-\widetilde{H}^{*}\|\leq&\sqrt{\sum_{m=1}^{M}\|\widetilde{H}_{T_{m},:}-\widetilde{H}^{*}_{T_{m},:}\|^{2}}.\end{split}$
(15)

We will now bound the spectral norm error $\|\widetilde{V}_{I_{m}\backslash
J_{m}^{(1)},:}-\widetilde{V}^{*}_{I_{m}\backslash J_{m}^{(1)},:}\|$ and
$\|\widetilde{H}_{T_{m},:}-\widetilde{H}^{*}_{T_{m},:}\|$ for each block $m$
in the following step.

#### Bounding the matching-induced errors

For the first block, by construction,
$\|\widetilde{V}_{I_{1},:}-\widetilde{V}^{*}_{I_{1},:}\|=\|\widehat{V}_{1}\widehat{W}_{1}-V_{1}^{*}\|$,
$\|\widetilde{H}_{T_{1},:}-\widetilde{H}^{*}_{T_{1},:}\|=\|\widehat{H}_{1}\widehat{W}_{1}-U_{1}^{*}\Lambda_{1}^{*}\|$.
For the $m$th block, let
$\varepsilon_{m}^{(2)}:=\big{\|}\widetilde{V}_{I_{m}\backslash
J_{m}^{(1)},:}-\widetilde{V}^{*}_{I_{m}\backslash J_{m}^{(1)},:}\big{\|}$.
Then one can show that

$\begin{split}\varepsilon_{m}^{(2)}\leq&\|\widehat{V}_{m}G_{m}-\widetilde{V}^{*}_{I_{m},:}\|\\\
=&\|\widehat{V}_{m}G_{m}-V_{m}^{*}B_{m}^{*}B_{1}^{*-1}\widehat{W}_{1}^{\top}\|\\\
=&\|\widehat{V}_{m}G_{m}-V_{m}^{*}\widehat{W}_{m}^{\top}G_{m}^{*}\|\\\
\leq&\|G_{m}-G_{m}^{*}\|+\|\widehat{V}_{m}\widehat{W}_{m}-V_{m}^{*}\|\|G_{m}^{*}\|\\\
=:&\varepsilon_{m}^{(3)}+\varepsilon_{m}^{(1)}\|G_{m}^{*}\|,\end{split}$ (16)

where $G_{m}^{*}$ was defined earlier in (13) as the oracle matching matrix
that aligns block $m$ with the first block, and
$\varepsilon_{m}^{(3)}:=\|G_{m}-G_{m}^{*}\|$,
$\varepsilon_{m}^{(1)}:=\|\widehat{V}_{m}\widehat{W}_{m}-V_{m}^{*}\|$.

While for $\|\widetilde{H}_{T_{m},:}-\widetilde{H}^{*}_{T_{m},:}\|$, we have

$\begin{split}\|\widetilde{H}_{T_{m},:}-\widetilde{H}^{*}_{T_{m},:}\|&=\|\widehat{H}_{m}G_{m}^{\top-1}-U_{m}^{*}\Lambda_{m}^{*}B_{m}^{*\top}(B_{1}^{*-1})^{\top}\widehat{W}_{1}^{\top}\|_{2}\\\
&=\|\widehat{H}_{m}G_{m}^{\top-1}-U_{m}^{*}\Lambda_{m}^{*}\widehat{W}_{m}^{\top}(G_{m}^{*\top})^{-1}\|\\\
&\leq\|X_{T_{m},I_{m}}^{*}\|\|G_{m}^{-1}-G_{m}^{*-1}\|\\\
&\quad+\|\widehat{H}_{m}\widehat{W}_{m}-U_{m}^{*}\Lambda_{m}^{*}\|(\|G_{m}^{*-1}\|+\|G_{m}^{-1}-G_{m}^{*-1}\|).\end{split}$
(17)

This would require an error bound for $\|G_{m}^{-1}-G_{m}^{*-1}\|$.

In the following, we first show bounds for $G_{m}-G_{m}^{*}$, and then use
them to derive error bounds for $\widetilde{V}_{I_{m}\backslash
J_{m}^{(1)},:}-\widetilde{V}^{*}_{I_{m}\backslash J_{m}^{(1)},:}$ and
$\widetilde{H}_{T_{m},:}-\widetilde{H}^{*}_{T_{m},:}$. First we note that
$G_{m}^{*}$ satisfies

$\begin{split}(V_{m}^{*})_{J_{m}^{(2)},:}\widehat{W}_{m}^{\top}G_{m}^{*}&=(V_{m}^{*})_{J_{m}^{(2)},:}B_{m}^{*}B_{1}^{*-1}\widehat{W}_{1}^{\top}\\\
&=V^{*}_{J_{m}^{(1)},:}B_{1}^{*-1}\widehat{W}_{1}^{\top}\\\
&=\widetilde{V}^{*}_{J_{m}^{(1)},:},\end{split}$

while $\widehat{G}_{m}$ is the OLS solution for transforming
$(\widehat{V}_{m})_{J_{m}^{(2)},:}$ to $\widetilde{V}_{J_{m}^{(1)},:}$. The
following lemma shows that the perturbation bound of the OLS solution:

######  Lemma 3.

Suppose $Y^{*}=Z^{*}G^{*}$ for $Z^{*}\in\mathbb{R}^{p\times r}$,
$G^{*}\in\mathbb{R}^{r\times r}$, and $G=(Z^{\top}Z)^{\dagger}Z^{\top}Y$ with
$Z$ and $Y$ of the same dimension as $Z^{*}$ and $Y^{*}$. If
$\|Z-Z^{*}\|\leq\varepsilon_{1}$, $\|Y-Y^{*}\|\leq\varepsilon_{2}$ and
$\sigma_{r}(Z^{*})>\varepsilon_{1}$, then we have

$\|G-G^{*}\|\leq\frac{\|G^{*}\|\varepsilon_{1}+\varepsilon_{2}}{\sigma_{r}(Z^{*})-\varepsilon_{1}}.$

Now we apply Lemma 3 with
$Z^{*}=(V_{m}^{*})_{J_{m}^{(2)},:}\widehat{W}_{m}^{\top}$, $G^{*}=G_{m}^{*}$,
$Y^{*}=\widetilde{V}^{*}_{J_{m}^{(1)},:}$;
$Z=(\widehat{V}_{m})_{J_{m}^{(2),:}}$, $Y=\widetilde{V}_{J_{m}^{(1)},:}$. To
verify the condition in Lemma 3, we note that

$\begin{split}\sigma_{r}((V_{m}^{*})_{J_{m}^{(2)},:}\widehat{W}_{m}^{\top})&=\sigma_{r}((V_{m}^{*})_{J_{m}^{(1)}})\geq\frac{\sigma_{r}(X^{*}_{T_{m},J_{m}^{(1)}})}{\|U_{m}^{*}\Lambda_{m}^{*}\|}\geq\frac{\sigma_{r}(X^{*}_{T_{m},J_{m}^{(1)}})}{\|X^{*}_{T_{m},I_{m}}\|}=\gamma_{m},\end{split}$

where the last equation is due to the definition of $\gamma_{m}$ in Section 3.
Meanwhile, due to our previous error bound (9) for each block, the first error
term in Lemma 3 satisfies

$\begin{split}\varepsilon_{1}=&\|(\widehat{V}_{m})_{J_{m}^{(2),:}}-(V_{m}^{*})_{J_{m}^{(2)},:}\widehat{W}_{m}^{\top}\|\\\
\leq&\|\widehat{V}_{m}\widehat{W}_{m}-V_{m}^{*}\|\\\
\leq&\frac{2\sqrt{2}\|E_{T_{m},I_{m}}\|}{\sigma_{r}(X^{*}_{T_{m},I_{m}})}.\end{split}$

Since the number of blocks $M\geq 2$, as long as the universal constant $C$ in
Assumption 2 is chosen to satisfy $C\geq 20\sqrt{2}$, we have
$\varepsilon_{1}\leq\frac{2\sqrt{2}\|E_{T_{m},I_{m}}\|}{\sigma_{r}(X^{*}_{T_{m},I_{m}})}\leq\frac{\gamma}{11}$,
which implies

$\gamma_{m}\geq\gamma\geq 11\varepsilon_{1}.$

Then, invoking Lemma 3, we obtain the following:

$\begin{split}\varepsilon_{m}^{(3)}&=\|G_{m}-G_{m}^{*}\|\\\
&\leq\frac{1.1}{\gamma_{m}}(\|G_{m}^{*}\|\|\widehat{V}_{m}\widehat{W}_{m}-V_{m}^{*}\|+\|\widetilde{V}_{J_{m}^{(1)},:}-\widetilde{V}^{*}_{J_{m}^{(1)},:}\|)\\\
&\leq\frac{1.1}{\gamma_{m}}(\beta_{\max}\varepsilon_{m}^{(1)}+\sum_{l=1}^{m-1}\varepsilon_{l}^{(2)}),\end{split}$
(18)

where we have applied the fact that
$\|G_{m}^{*}\|=\|B_{m}^{*-1}B_{1}^{*}\|\leq\beta_{\max}$, and
$J_{m}^{(1)}\subset\cup_{l=1}^{m-1}(I_{l}\backslash J_{l}^{(1)})$. We have
also utilized the definitions
$\varepsilon_{m}^{(1)}=\|\widehat{V}_{m}\widehat{W}_{m}-V_{m}^{*}\|$, the
error term for each block SVD; and
$\varepsilon_{l}^{(2)}=\|\widetilde{V}_{I_{l}\backslash
J_{l}^{(1)},:}-\widetilde{V}_{I_{l}\backslash J_{l}^{(1)},:}^{*}\|$, the final
error of the $l$’s block after the matching step. Putting (16) and (18)
together, and using an induction proof, we can show the following claim.

######  Claim 1.

Let $\varepsilon$ denote the maximum block-SVD error $\max_{1\leq m\leq
M}\varepsilon_{m}^{(1)}$, and $\widetilde{\gamma}_{m}=\frac{1.1}{\gamma_{m}}$.
Then $\varepsilon_{1}^{(2)}\leq\varepsilon$, and for $m=2,\dots,M$,

$\begin{split}\varepsilon_{m}^{(2)}\leq&(2\beta_{\max}+1)\widetilde{\gamma}_{m}\prod_{l=2}^{m-1}(\widetilde{\gamma}_{l}+1)\varepsilon,\\\
\sum_{l=1}^{m}\varepsilon_{m}^{(2)}\leq&(2\beta_{\max}+1)\prod_{l=2}^{m}(\widetilde{\gamma}_{l}+1)\varepsilon-\beta_{\max}\varepsilon.\end{split}$
(19)

We will assume this claim to be true for the rest of the arguments while
leaving its proof in the end. Combining Claim 1 and (18), for $m\geq 2$, we
have

$\begin{split}\|G_{m}-G_{m}^{*}\|=\varepsilon_{m}^{(3)}\leq&\widetilde{\gamma}_{m}\beta_{\max}\varepsilon+\tilde{\gamma}_{m}\sum_{l=1}^{m-1}\varepsilon_{l}^{(2)}\\\
\leq&(2\beta_{\max}+1)\widetilde{\gamma}_{m}\prod_{l=2}^{m-1}(\widetilde{\gamma}_{l}+1)\varepsilon.\end{split}$
(20)

Recall (15), the estimation error for the whole $\widetilde{V}$ can be bounded
as follows:

$\begin{split}\|\widetilde{V}-\widetilde{V}^{*}\|&\leq\sum_{m=1}^{M}\varepsilon_{m}^{(2)}\\\
&\leq(2\beta_{\max}+1)\prod_{l=2}^{M}(\widetilde{\gamma}_{l}+1)\varepsilon.\end{split}$
(21)

Now we turn to the left singular space: $\|\widetilde{H}-\widetilde{H}^{*}\|$.
We still consider the contributed error
$\|\widetilde{H}_{T_{m},:}-\widetilde{H}^{*}_{T_{m},:}\|$ of each block first,
which then requires an error bound for $\|G_{m}^{-1}-G_{m}^{*-1}\|$. Recall
our bound for $\|G_{m}-G_{m}^{*}\|$ in (20). We can further show that

$\begin{split}\|G_{m}-G_{m}^{*}\|&\leq
2\sqrt{2}(2\beta_{\max}+1)\prod_{l=2}^{m}\left(\frac{1.1}{\gamma_{l}}+1\right)\max_{m}\frac{\|E_{T_{m},I_{m}}\|}{\sigma_{r}(X^{*}_{T_{m},I_{m}})}\\\
&\leq
4\sqrt{2}(\beta_{\max}+1)\prod_{l=2}^{M}\left(\frac{1.1}{\gamma_{l}}+1\right)\max_{m}\frac{\|E_{T_{m},I_{m}}\|}{\sigma_{r}(X^{*}_{T_{m},I_{m}})}\\\
&\leq\frac{1}{2}\beta_{\min},\end{split}$

where the first inequality is due to (20), the definition
$\varepsilon=\max_{m}\epsilon_{m}^{(1)}$ in Claim 1, (9), and
$\widetilde{\gamma}_{m}=\frac{1.1}{\gamma_{m}}$; the third inequality makes
use of Assumption 2 with a sufficiently large constant $C>0$. Here,
$\beta_{\min}=\min_{1\leq m\leq
M}\sigma_{r}(B_{m}^{*-1}B_{1}^{*})\leq\sigma_{r}(G_{m}^{*})$. Then through
some calculations, one can show that

$\begin{split}\|G_{m}^{-1}-G_{m}^{*-1}\|&=\|G_{m}^{-1}(G_{m}^{*}-G_{m})G_{m}^{*-1}\|\\\
&\leq\frac{\|G_{m}^{*}-G_{m}\|}{\sigma_{r}(G_{m})\sigma_{r}(G_{m}^{*})}\\\
&\leq 2\beta_{\min}^{-2}\|G_{m}-G_{m}^{*}\|\\\
&\leq\beta_{\min}^{-1},\end{split}$

where the third inequality utilizes the fact that
$\sigma_{r}(G_{m})\geq\sigma_{r}(G_{m}^{*})-\|G_{m}-G_{m}^{*}\|\geq\frac{1}{2}\sigma_{r}(G_{m}^{*})$,
and the last inequality is due to
$\|G_{m}-G_{m}^{*}\|\leq\frac{1}{2}\beta_{\min}$. Therefore, we can plug in
the error bound above into (17):

$\begin{split}\|\widetilde{H}_{T_{m},:}-\widetilde{H}_{T_{m},:}^{*}\|&\leq\|U_{m}^{*}\Lambda_{m}^{*}\|\|G_{m}^{-1}-G_{m}^{*-1}\|+2\beta_{\min}^{-1}\|\widehat{H}_{m}\widehat{W}_{m}-U_{m}^{*}\Lambda_{m}^{*}\|\\\
&\leq\frac{2\|\Lambda_{m}^{*}\|}{\beta_{\min}^{2}}\left(\|G_{m}-G_{m}^{*}\|+\frac{4\beta_{\min}\|E_{T_{m},I_{m}}\|}{\sigma_{r}(X^{*}_{T_{m},I_{m}})}\right)\\\
&\leq\frac{C\|\Lambda_{m}^{*}\|}{\beta_{\min}^{2}}(\beta_{\max}+1)\widetilde{\gamma}_{m}\prod_{l=2}^{m-2}(\widetilde{\gamma}_{l}+1)\max_{m}\frac{\|E_{T_{m},I_{m}}\|}{\sigma_{r}(X^{*}_{T_{m},I_{m}})},\end{split}$

where the second inequality is due to (10); the third inequality is due to
(20) and the fact that $\varepsilon\leq
2\sqrt{2}\max_{m}\frac{\|E_{T_{m},I_{m}}\|}{\sigma_{r}(X^{*}_{T_{m},I_{m}})}$.
Combining all the bounds for
$\|\widetilde{H}_{T_{m},:}-\widetilde{H}_{T_{m},:}^{*}\|$ over $m=1,\dots,M$,
we have

$\begin{split}\|\widetilde{H}-\widetilde{H}^{*}\|&\leq\frac{C\max_{m}\|\Lambda_{m}^{*}\|}{\beta_{\min}^{2}}(\beta_{\max}+1)\sqrt{\sum_{m=1}^{M}\prod_{l=2}^{m}(\widetilde{\gamma}_{l}+1)^{2}}\max_{m}\frac{\|E_{T_{m},I_{m}}\|}{\sigma_{r}(X^{*}_{T_{m},I_{m}})}\\\
&\leq
C\beta_{\min}^{-2}(\beta_{\max}+1)\prod_{l=2}^{M}\left(\frac{1.1}{\gamma_{l}}+1\right)\max_{m}\|\Lambda_{m}^{*}\|\max_{m}\frac{\|E_{T_{m},I_{m}}\|}{\sigma_{r}(X^{*}_{T_{m},I_{m}})}.\end{split}$
(22)

In the last line above, we have utilized the fact that

$\begin{split}&\sqrt{\sum_{m=1}^{M}\prod_{l=2}^{m}(\widetilde{\gamma}_{l}+1)^{2}}\\\
\leq&\frac{\prod_{l=2}^{M}(\widetilde{\gamma}_{l}+1)^{2}(\min_{l}\widetilde{\gamma}_{l}+1)^{2}-1}{(\min_{l}\widetilde{\gamma}_{l}+1)^{2}-1}\\\
\leq&C\prod_{l=2}^{M}(\widetilde{\gamma}_{l}+1)^{2},\end{split}$

where the second line is due to Lemma 2, and the last line is due to
$\min_{l}\widetilde{\gamma}_{l}\geq 1.1$.

Now we would like to combine (21) and (22). Note that Assumption 2 suggests
the following:

$\begin{split}\|\widetilde{V}-\widetilde{V}^{*}\|&\leq(2\beta_{\max}+1)\prod_{l=2}^{M}(\widetilde{\gamma}_{l}+1)\varepsilon\\\
&\leq
2\sqrt{2}(2\beta_{\max}+1)\prod_{l=2}^{M}\left(\frac{1.1}{\gamma_{l}}+1\right)\max_{m}\frac{\|E_{T_{m},I_{m}}\|}{\sigma_{r}(X^{*}_{T_{m},I_{m}})}\\\
&\leq C,\end{split}$

where the second inequality is due to the fact that $\varepsilon\leq
2\sqrt{2}\max_{m}\frac{\|E_{T_{m},I_{m}}\|}{\sigma_{r}(X^{*}_{T_{m},I_{m}})}$,
and the last line is due to Assumption 2. Finally, we also note that the
spectral norm of $\widetilde{V}^{*}$ and $\widetilde{H}^{*}$ can be bounded as
follows:

$\begin{split}\|\widetilde{V}^{*}\|&=\|V^{*}B_{1}^{*}\|\leq\|B_{1}^{*}\|=\|\Lambda^{*}U_{T_{1},:}^{*\top}U_{1}^{*}\Lambda_{1}^{*-1}\|\leq\frac{\|X^{*}\|}{\sigma_{r}(\Lambda^{*}_{1})},\\\
\|\widetilde{H}^{*}\|&=\|U^{*}\Lambda^{*}(B_{1}^{*-1})^{\top}\|\leq\|\Lambda^{*}\|\|B_{1}^{*-1}\|=\|\Lambda^{*}\|\|V_{1}^{*\top}V_{I_{1},:}^{*}\|\leq\|X^{*}\|,\end{split}$

where we have applied (11) and (12), which suggest
$B_{1}^{*}=\Lambda^{*}U_{T_{1},:}^{*\top}U_{1}^{*}\Lambda_{1}^{*-1}$,
$B_{1}^{*-1}=V_{1}^{*\top}V_{I_{1},:}^{*}$. Therefore, we can now bound
$\|\widetilde{H}\widetilde{V}-X^{*}\|$ as follows:

$\begin{split}\|\widetilde{H}\widetilde{V}-X^{*}\|&\leq\|\widetilde{H}-\widetilde{H}^{*}\|(\|\widetilde{V}^{*}\|+\|\widetilde{V}-\widetilde{V}^{*}\|)+\|\widetilde{V}-\widetilde{V}^{*}\|\|\widetilde{H}^{*}\|\\\
&\leq
C\beta_{\min}^{-2}(\beta_{\max}+1)\max_{m}\frac{\|\Lambda_{m}^{*}\|}{\sigma_{r}(\Lambda_{1}^{*})}\prod_{l=2}^{M}\left(\frac{1.1}{\gamma_{l}}+1\right)\max_{m}\frac{\|E_{T_{m},I_{m}}\|}{\sigma_{r}(X^{*}_{T_{m},I_{m}})}\|X^{*}\|\\\
&\leq
C\frac{\alpha_{\max}(\beta_{\max}+1)}{\beta_{\min}^{2}}\prod_{l=2}^{M}\left(\frac{1.1}{\gamma_{l}}+1\right)\max_{m}\frac{\|E_{T_{m},I_{m}}\|}{\sigma_{r}(X^{*}_{T_{m},I_{m}})}\|X^{*}\|.\end{split}$

###### Proof of Claim 1.

We first combine (16) and (18) and obtain the following: for $2\leq m\leq M$,

$\begin{split}\varepsilon_{m}^{(2)}\leq\beta_{\max}(\widetilde{\gamma}_{m}+1)\varepsilon_{m}^{(1)}+\widetilde{\gamma}_{m}\sum_{l=1}^{m-1}\varepsilon_{l}^{(2)}.\end{split}$
(23)

Let $\delta_{m}=\sum_{l=1}^{m}\varepsilon_{l}^{(2)}$. Then (23) can be
rewritten as

$\delta_{m}\leq\beta_{\max}(\widetilde{\gamma}_{m}+1)\varepsilon+(\widetilde{\gamma}_{m}+1)\delta_{m-1}.$
(24)

We first note that when $m=1$, $J_{m}^{(1)}=\emptyset$ and hence

$\varepsilon_{1}^{(2)}=\|\widetilde{V}_{I_{1},:}-\widetilde{V}^{*}_{I_{1},:}\|=\|\widehat{V}_{1}\widehat{W}_{1}-V_{1}^{*}\|=\varepsilon^{(1)}_{1}.$

Immediately, we have $\delta_{1}=\varepsilon_{1}^{(2)}\leq\varepsilon$.

In the following, we will prove by induction that for $m\geq 2$.

$\delta_{m}\leq\left[\sum_{k=2}^{m}\prod_{l=k}^{m}(\widetilde{\gamma}_{l}+1)\beta_{\max}+\prod_{l=2}^{m}(\widetilde{\gamma}_{l}+1)\right]\varepsilon.$
(25)

When $m=2$, (24) implies that
$\delta_{2}\leq\beta_{\max}(\widetilde{\gamma}_{2}+1)\varepsilon+(\widetilde{\gamma}_{2}+1)\varepsilon$,
which satisfies (25). Suppose that (25) holds for $2\leq m\leq i$, then when
$m=i+1$, (24) suggests that

$\begin{split}\delta_{i+1}\leq&\beta_{\max}(\widetilde{\gamma}_{i+1}+1)\varepsilon+(\widetilde{\gamma}_{i+1}+1)\delta_{i}\\\
\leq&\beta_{\max}(\widetilde{\gamma}_{i+1}+1)\varepsilon+(\widetilde{\gamma}_{i+1}+1)\left[\sum_{k=2}^{i}\prod_{l=k}^{i}(\widetilde{\gamma}_{l}+1)\beta_{\max}+\prod_{l=2}^{i}(\widetilde{\gamma}_{l}+1)\right]\varepsilon\\\
\leq&\left[\beta_{\max}(\widetilde{\gamma}_{i+1}+1)+\sum_{k=2}^{i}\prod_{l=k}^{i+1}(\widetilde{\gamma}_{l}+1)\beta_{\max}+\prod_{l=2}^{i+1}(\widetilde{\gamma}_{l}+1)\right]\varepsilon\\\
\leq&\left[\sum_{k=2}^{i+1}\prod_{l=k}^{i+1}(\widetilde{\gamma}_{l}+1)\beta_{\max}+\prod_{l=2}^{i+1}(\widetilde{\gamma}_{l}+1)\right]\varepsilon,\end{split}$

and hence (25) holds for all $2\leq m\leq M$.

Now we combine (25) and (23) to show that

$\varepsilon_{m}^{(2)}\leq\left[(\widetilde{\gamma}_{m}+1)\beta_{\max}+\widetilde{\gamma}_{m}\sum_{k=2}^{m-1}\prod_{l=k}^{m-1}(\widetilde{\gamma}_{l}+1)\beta_{\max}+\widetilde{\gamma}_{m}\prod_{l=2}^{m-1}(\widetilde{\gamma}_{l}+1)\right]\varepsilon.$
(26)

We now apply Lemma 2 to show an upper bound for
$\sum_{k=2}^{m}\prod_{l=k}^{m}(\widetilde{\gamma}_{l}+1)$, which can help
simplify both (25) and (26). Let
$\widetilde{\gamma}_{\min}=\min_{l}\widetilde{\gamma}_{l}$. Lemma 2 then
implies

$\sum_{k=2}^{m}\prod_{l=k}^{m}(\widetilde{\gamma}_{l}+1)\leq\widetilde{\gamma}_{\min}^{-1}\left[\prod_{l=2}^{m}(\widetilde{\gamma}_{l}+1)(\widetilde{\gamma}_{\min}+1)-(\widetilde{\gamma}_{m}+1)\right].$

Since $0<\gamma_{l}\leq 1$ for all $2\leq l\leq M$, we have
$\widetilde{\gamma}_{\min}=\min_{l}\frac{1.1}{\gamma_{l}}\geq 1.1$. Therefore,
(25) and (26) can be simplified as follows:

$\begin{split}\delta_{m}\leq&(2\beta_{\max}+1)\prod_{l=2}^{m}(\widetilde{\gamma}_{l}+1)\varepsilon-\beta_{\max}\varepsilon,\\\
\varepsilon_{m}^{(2)}\leq&(2\beta_{\max}+1)\widetilde{\gamma}_{m}\prod_{l=2}^{m-1}(\widetilde{\gamma}_{l}+1)\varepsilon+(\widetilde{\gamma}_{m}+1)\beta_{\max}\varepsilon-\frac{(\widetilde{\gamma}_{m-1}+1)\widetilde{\gamma}_{m}}{\widetilde{\gamma}_{\min}}\beta_{\max}\varepsilon\\\
\leq&(2\beta_{\max}+1)\widetilde{\gamma}_{m}\prod_{l=2}^{m-1}(\widetilde{\gamma}_{l}+1)\varepsilon,\end{split}$

where the last line is due to the fact that
$\frac{(\widetilde{\gamma}_{m-1}+1)\widetilde{\gamma}_{m})}{\widetilde{\gamma}_{\min}}=\frac{\widetilde{\gamma}_{m-1}\widetilde{\gamma}_{m}}{\widetilde{\gamma}_{\min}}+\frac{\widetilde{\gamma}_{m}}{\widetilde{\gamma}_{\min}}\geq\widetilde{\gamma}_{m}+1$.
Therefore, the proof of Claim 1 is complete. ∎

###  A.2 Proofs of Auxiliary Lemmas

###### Proof of Lemma 1.

Our proof follows very similar arguments to Löffler et al.,  (2021), but we
still state most steps for clarity and completeness. Since
$\widetilde{H}\widetilde{V}^{\top}$ and $X^{*}$ are both of rank $r$, the
spectral norm error bound
$\|\widetilde{H}\widetilde{V}^{\top}-X^{*}\|_{2}\leq\varepsilon$ implies the
following Frobenious error bound:

$\|\widetilde{H}\widetilde{V}^{\top}-X^{*}\|_{F}\leq\sqrt{2r}\varepsilon.$

Let $\widehat{X}\in\mathbb{R}^{p\times n}$ be the estimated center matrix with
each column representing the estimated center mean of the corresponding sample
$i$: $\widehat{X}_{:,i}=\widetilde{U}\widehat{c}_{z_{i}}$, where
$\widehat{c}_{j}\in\mathbb{R}^{r}$ and $\widehat{z}_{i}\in[K]$ are the global
minimizer of the k-means step in Algorithm 1. By the definition of the k-means
optimization problem, one has

$\begin{split}\|\widetilde{H}\widetilde{V}^{\top}-\widehat{X}\|_{F}=\|\widetilde{U}^{\top}\widetilde{H}\widetilde{V}^{\top}-\widetilde{U}^{\top}\widehat{X}\|_{F}\leq&\|\widetilde{U}^{\top}\widetilde{H}\widetilde{V}^{\top}-\widetilde{U}^{\top}X^{*}\|_{F}\\\
=&\|\widetilde{H}\widetilde{V}^{\top}-X^{*}\|_{F},\end{split}$

which then implies $\|\widehat{X}-X^{*}\|_{F}\leq 2\sqrt{2r}\varepsilon$.
Define the sample set
$S=\\{i\in[n]:\|\widehat{X}_{:,i}-X^{*}_{:,i}\|_{2}>\frac{\Delta}{2}\\}$. The
Frobenious norm error bound for $\widehat{X}-X^{*}$ then implies an upper
bound for the cardinality of set $S$:

$|S|\leq\frac{\|\widehat{X}-X^{*}\|_{F}^{2}}{\Delta^{2}/4}\leq\frac{32r\varepsilon^{2}}{\Delta^{2}}.$

In addition, since we have assumed $\varepsilon<\sqrt{\frac{\beta
n}{rK}}\frac{\Delta}{8}$, the inequality above also implies
$|S|\leq\frac{\beta n}{2K}$. Following the same arguments as in the proof of
Lemma 4.2 in Löffler et al.,  (2021), we then have
$\ell(\widehat{z},z^{*})\leq\frac{1}{n}|S|\leq\frac{32r\varepsilon^{2}}{n\Delta^{2}}$,
and $\min_{\phi\in\Phi}\max_{1\leq j\leq
K}\|\widehat{U}\widehat{c}_{j}-\theta^{*}_{\phi(j)}\|_{2}^{2}\leq\frac{16Kr\varepsilon^{2}}{\beta
n}$. ∎

###### Proof of Lemma 3.

Since we have assumed $\sigma_{r}(Z^{*})>\varepsilon_{1}$, by Weyl’s
inequality, one has $\sigma_{r}(Z)\geq\sigma_{r}(Z^{*})-\varepsilon_{1}>0$,
which suggests that $Z^{\top}Z\in\mathbb{R}^{r\times r}$ is invertible. Now we
use some rearrangement of the least squares formulation of $G$:

$\begin{split}G&=(Z^{\top}Z)^{-1}Z^{\top}Y\\\
&=(Z^{\top}Z)^{-1}Z^{\top}(Z^{*}G^{*}+Y-Y^{*})\\\
&=G^{*}+(Z^{\top}Z)^{-1}Z^{\top}[Y-Y^{*}-(Z-Z^{*})G^{*}],\end{split}$

Thus we only need to bound
$\|(Z^{\top}Z)^{-1}Z^{\top}[(Z-Z^{*})G^{*}+Y-Y^{*}]\|$. Recognizing the fact
that
$\|(Z^{\top}Z)^{-1}Z^{\top}\|=\sigma_{r}^{-1}(Z)\leq(\sigma_{r}(Z^{*})-\varepsilon_{1})^{-1}$,
one can show that

$\begin{split}&\|(Z^{\top}Z)^{-1}Z^{\top}[(Z-Z^{*})G^{*}+Y-Y^{*}]\|\\\
\leq&\frac{\|(Z-Z^{*})G^{*}+Y-Y^{*}\|}{\sigma_{r}(Z^{*})-\varepsilon_{1}}\\\
\leq&\frac{\|G^{*}\|\varepsilon_{1}+\varepsilon_{2}}{\sigma_{r}(Z^{*})-\varepsilon_{1}}.\end{split}$

The proof is now complete. ∎

###  A.3 Proofs of Corollary 1 and Theorem 2

###### Proof of Corollary 1.

We first establish a connection between the separation distance and the
blockwise SNR in terms of minimum singular value. In particular, for each
block $1\leq m\leq M$,
$X_{T_{m},I_{m}}^{*}=\Theta^{*}_{T_{m},:}F^{*\top}_{I_{m},:}$, where
$\Theta^{*}_{T_{m},:}\in\mathbb{R}^{p_{m}\times K}$,
$F^{*}_{I_{m},:}\in\mathbb{R}^{n_{m}\times K}$. For any $u\in\mathbb{R}^{K}$
with $\|u\|_{2}=1$, we have
$\|F^{*}_{I_{m},:}u\|_{2}=\sqrt{\sum_{j=1}^{K}u_{j}^{2}n_{j,m}}\geq\sqrt{\min_{j}n_{j,m}}$,
and hence $\sigma_{\min}(F^{*}_{I_{m},:})\geq\sqrt{\min_{j}n_{j,m}}$. In
addition, recall that $\Delta_{m}=\min_{j\neq
j^{\prime}}\|\theta_{T_{m},j}-\theta_{T_{m},j^{\prime}}\|_{2}$. Let
$(j_{1},j_{2})=\operatorname*{arg\,min}_{j_{1}\neq
j_{2}}\|\theta_{T_{m},j_{1}}-\theta_{T_{m},j_{2}}\|_{2}$,
$u=\frac{1}{\sqrt{2}}(e_{j_{1}}-e_{j_{2}})$, where $e_{j_{1}},\,e_{j_{2}}$ are
canonical vectors. Then we have $\|u\|_{2}=1$, and
$\|\Theta^{*}_{T_{m},:}u\|_{2}=\frac{1}{\sqrt{2}}\Delta_{m}$, which implies
$\|\Theta^{*}_{T_{m},:}\|\geq\frac{1}{\sqrt{2}}\Delta_{m}$. Therefore, we can
lower bound the $r$th singular value of $X_{T_{m},I_{m}}^{*}$ as follows:

$\begin{split}\sigma_{r}(X^{*}_{T_{m},I_{m}})\geq&\sigma_{\min}(F_{I_{m},:}^{*})\sigma_{r}(\Theta^{*}_{T_{m},:})\\\
\geq&\sqrt{\min_{j}n_{j,m}}\|\Theta^{*}_{T_{m},:}\|\kappa_{r}(\Theta^{*}_{T_{m},:})\\\
\geq&\sqrt{\frac{\min_{j}n_{j,m}}{2}}\frac{\Delta_{m}}{\kappa(\Theta^{*}_{T_{m},:})}.\end{split}$

Thus Assumption 3 implies Assumption 2. Combining the lower bound above with
Theorem 1 leads to Theorem 2. ∎

###### Proof of Theorem 2.

Recall that the noise matrix $E=A_{r}ZA_{c}^{\top}$, where
$A_{r}\in\mathbb{R}^{n\times n^{\prime}}$, $A_{c}\in\mathbb{R}^{p\times
p^{\prime}}$, $Z\in\mathbb{R}^{n^{\prime}\times p^{\prime}}$ has independent
zero-mean sub-Gaussian($\sigma$) entries. We can first bound the spectral norm
of each submatrix $E_{T_{m},I_{m}}=(A_{r})_{T_{m},:}Z(A_{c})_{I_{m},:}^{\top}$
by a function of $\sigma$, spectral norms $\|A_{r}\|$, $\|A_{c}\|$, and its
dimensions. In particular, we apply a technical lemma from Zhou et al.,
(2022) (Lemma 8.2), with $A=\frac{1}{\|(A_{r})_{T_{m},:}\|}(A_{r})_{T_{m},:}$,
$B=\frac{1}{\|(A_{c})_{I_{m},:}\|}(A_{c})_{I_{m},:}^{\top}$,
$t=C(n_{m}+p_{m})$ for a sufficiently large universal constant. Then one can
write

$\begin{split}\mathbb{P}(\|E_{T_{m},I_{m}}\|\geq&2\sigma\|(A_{r})_{T_{m},:}\|\|(A_{c})_{I_{m},:}\|\sqrt{2p_{m}+n_{m}})\\\
\leq&2\cdot
5^{n_{m}}\exp\\{-c_{1}\min\\{\frac{(n_{m}+p_{m})^{2}}{p_{m}},n_{m}+p_{m}\\}\\}\\\
\leq&2\exp\\{-c_{2}(n_{m}+p_{m})\\}.\end{split}$

Therefore, with probability at least $1-2\exp\\{-c(n_{m}+p_{m})\\}$,
$\|E_{T_{m},I_{m}}\|\leq C\sigma\|A_{r}\|\|A_{c}\|(\sqrt{p_{m}}+\sqrt{n_{m}})$
for some universal constants $c,\,C>0$. We then note that (7) together with
Assumption 4 implies Assumption 3; the misclustering rate in (4) and the bound
above for $\|E_{T_{m},I_{m}}\|$ immediately implies Theorem 2. ∎

Figure 5: Illustration of data generation process for Gaussian mixture model
simulations.

###  A.4 Proof of Theorem 3

We first prove that when Algorithm 1 is applied on data patches
$\\{X_{T_{\hat{\pi}_{m}},I_{\hat{\pi}_{m}}}\\}_{m=1}^{M}$ with data-driven
patch ordering $\hat{\pi}$, and when the factor
$\prod_{m=2}^{M}\left(1.1\gamma_{m}^{-1}+1\right)$ in Assumption 2 and Theorem
1 is substituted by $\prod_{m=2}^{M}(1.1\gamma_{m}^{-1}(\pi^{*})+1)$, Theorem
1 still holds. In fact, it suffices to show that the matching error factor for
the data-driven selected ordering $\hat{\pi}$ satisfies

$\prod_{m=2}^{M}\left(1.1\gamma_{m}^{-1}(\hat{\pi})+1\right)\leq
e^{2}\prod_{m=2}^{M}\left(1.1\gamma_{m}^{-1}(\pi^{*})+1\right);$ (27)

then we can simply invoke the original Theorem 1 on data patches
$\\{X_{T_{\hat{\pi}(m)},I_{\hat{\pi}(m)}}\\}_{m=1}^{M}$ to finish the proof.

Let $\hat{S}(\pi)=\prod_{m=2}^{M}s(\pi(m),\cup_{l<m}\pi(l))$ be the finite-
sample score for a given ordering $\pi$; and let
$S^{*}(\pi)=\prod_{m=2}^{M}(1.1\gamma_{m}^{-1}(\pi)+1)^{-1}$ be its population
counterpart. To show (27), it might be tempting to establish an error bound
between the finite sample score function $\hat{S}(\pi)$ and its population
counterpart $S^{*}(\pi)$ for all $\pi$. However, there possibly exist some
patch orderings $\pi$ such that the overlapping set $J_{m}^{(1)}(\pi)$ is
tiny; and for these terrible orderings, $\hat{S}(\pi)$ is not close to
$S^{*}(\pi)$ at all. Therefore, one key step in our proof is to show that the
selected $\hat{\pi}$ cannot be too bad under the given assumptions. For
notational convenience, we let
$\overline{\mathrm{SNR}}(m,\pi)=\frac{\big{\|}X^{*}_{T_{\pi(m)},I_{\pi(m)}}\big{\|}}{\big{\|}E_{T_{\pi(m)},I_{\pi(m)}}\big{\|}}$,
$\underline{\mathrm{SNR}}(m,\pi)=\frac{\sigma_{r}\big{(}X^{*}_{T_{\pi(m)},J_{m}^{(1)}(\pi)}\big{)}}{\big{\|}E_{T_{\pi(m)},J_{m}^{(1)}(\pi)}\big{\|}}$.
We first note that Assumption 2 on the oracle patch ordering implies that

$\frac{\sigma_{r}(X^{*}_{T_{m},I_{m}})}{\|E_{T_{m},I_{m}}\|}\geq\frac{C\alpha_{\max}(\beta_{\max}+1)}{\beta_{\min}^{2}}\prod_{l=2}^{M}\left(1.1\gamma_{l}^{-1}(\pi^{*})+1\right)\sqrt{\frac{rK\max_{j}n_{j}}{\min_{j}n_{j}}},$
(28)

which then implies that for any $1\leq m\leq M$,

$\overline{\mathrm{SNR}}(m,\pi^{*})\geq 2C\times 2.1^{M-1}\sqrt{rK}\geq
4.4CM-2.9C,$ (29)

where we have utilized the fact that $\alpha_{\max},\,\beta_{\max}\geq 1$,
$0<\beta_{\min},\gamma_{m}(\pi^{*})\leq 1$, $K\geq 2$, and the second
inequality holds due to the property of the exponential function of $M\geq 2$.
Similarly, we have

$\underline{\mathrm{SNR}}(m,\pi^{*})\geq\gamma_{m}(\pi^{*})\overline{\mathrm{SNR}}(m,\pi^{*})\geq
2.2C\times 2.1^{M-2}\sqrt{rK}\geq 2.3CM-1.6C.$ (30)

By Weyl’s inequality, we have

$\begin{split}\frac{\|X_{T_{\pi^{*}(m)},I_{\pi^{*}(m)}}\|}{\sigma_{r}(X_{T_{\pi^{*}(m)},J_{m}^{(1)}(\pi^{*})})}\leq&\frac{\|X^{*}_{T_{\pi^{*}(m)},I_{\pi^{*}(m)}}\|+\|E_{T_{\pi^{*}(m)},I_{\pi^{*}(m)}}\|}{\sigma_{r}(X^{*}_{T_{\pi^{*}(m)},J_{m}^{(1)}(\pi^{*})})-\|E_{T_{\pi^{*}(m)},J_{m}^{(1)}(\pi^{*})}\|}\\\
\leq&\frac{\|X^{*}_{T_{\pi^{*}(m)},I_{\pi^{*}(m)}}\|}{\sigma_{r}(X^{*}_{T_{\pi^{*}(m)},J_{m}^{(1)}(\pi^{*})})}\left(1+\frac{1}{\overline{\mathrm{SNR}}(m,\pi^{*})}\right)\left(1-\frac{1}{\underline{\mathrm{SNR}}(m,\pi^{*})}\right)^{-1},\end{split}$
(31)

and the inflated error factor can be further written as

$\begin{split}\left(1+\frac{1}{\overline{\mathrm{SNR}}(m,\pi^{*})}\right)\left(1-\frac{1}{\underline{\mathrm{SNR}}(m,\pi^{*})}\right)^{-1}=1+\left(\frac{1-\overline{\mathrm{SNR}}^{-1}(m,\hat{\pi})}{\overline{\mathrm{SNR}}^{-1}(m,\pi^{*})+\underline{\mathrm{SNR}}^{-1}(m,\pi^{*})}\right)^{-1}.\end{split}$
(32)

By setting the constant $C>0$ sufficiently large in Assumption 2, we can
guarantee that
$\overline{\mathrm{SNR}}(m,\pi^{*}),\,\underline{\mathrm{SNR}}(m,\pi^{*})\geq
10$, and hence (29) and (30) implies

$\begin{split}\left(1+\frac{1}{\overline{\mathrm{SNR}}(m,\pi^{*})}\right)\left(1-\frac{1}{\underline{\mathrm{SNR}}(m,\pi^{*})}\right)^{-1}\leq&1+\frac{1}{M-1}.\end{split}$

Therefore, the finite score of the oracle ordering satisfies

$\begin{split}\hat{S}(\pi^{*})\geq\prod_{m=2}^{M}\left(\frac{1.1\|X^{*}_{T_{\pi^{*}(m)},I_{\pi^{*}(m)}}\|}{\sigma_{r}(X^{*}_{T_{\pi^{*}(m)},J_{m}^{(1)}(\pi^{*})})}+1\right)^{-1}\big{(}1+\frac{1}{M-1}\big{)}^{M-1}\leq
e^{-1}\prod_{m=2}^{M}(1.1\gamma_{m}^{-1}(\pi^{*})+1)^{-1}.\end{split}$

Since the selected patch ordering maximizes the overlapping score function, we
have

$\begin{split}\hat{S}(\hat{\pi})\geq\hat{S}(\pi^{*})\geq
e^{-1}\prod_{m=2}^{M}(1.1\gamma_{m}^{-1}(\pi^{*})+1)^{-1}.\end{split}$ (33)

As mentioned earlier, in order to show (27), we would like to establish a
connection between $\hat{S}(\hat{\pi})$ and $S^{*}(\hat{\pi})$; to achieve
this, we also need to ensure sufficiently large
$\underline{\mathrm{SNR}}(m,\hat{\pi})$. We plug in the definition of
$\hat{S}(\hat{\pi})$ into (33) and apply Assumption 2:

$\begin{split}\hat{S}(\hat{\pi})=\prod_{l=2}^{M}\left(\frac{1.1\|X_{T_{\hat{\pi}(l)},I_{\hat{\pi}(l)}}\|}{\sigma_{r}(X_{T_{\hat{\pi}(l)},J_{l}^{(1)}(\hat{\pi})})}+1\right)^{-1}\geq\frac{2C\sqrt{rK}\|E_{T_{m},I_{m}}\|}{\sigma_{r}(X^{*}_{T_{m},I_{m}})},\end{split}$
(34)

which holds for all $1\leq m\leq M$. Since
$\frac{\|X_{T_{\hat{\pi}(l)},I_{\hat{\pi}(l)}}\|}{\sigma_{r}(X_{T_{\hat{\pi}(l)},J_{l}^{(1)}(\hat{\pi})})}\geq
1$ for all $2\leq l\leq M$, (34) implies that

$\begin{split}&\frac{1.1\|X_{T_{\hat{\pi}(m)},I_{\hat{\pi}(m)}}\|}{\sigma_{r}(X_{T_{\hat{\pi}(m)},J_{m}^{(1)}(\hat{\pi})})}+1\leq\frac{\sigma_{r}(X^{*}_{T_{\hat{\pi}(m)},I_{\hat{\pi}(m)}})}{2C\times
2.1^{M-2}\sqrt{rK}\|E_{T_{\hat{\pi}(m)},I_{\hat{\pi}(m)}}\|},\\\
\Rightarrow&\frac{1.1\left(\|X^{*}_{T_{\hat{\pi}(m)},I_{\hat{\pi}(m)}}\|-\|E_{T_{\hat{\pi}(m)},I_{\hat{\pi}(m)}}\|\right)}{\sigma_{r}(X_{T_{\hat{\pi}(m)},J_{m}^{(1)}(\hat{\pi})})}\leq\frac{\sigma_{r}(X^{*}_{T_{\hat{\pi}(m)},I_{\hat{\pi}(m)}})-\|E_{T_{\hat{\pi}(m)},I_{\hat{\pi}(m)}}\|}{2C\times
2.1^{M-2}\sqrt{rK}\|E_{T_{\hat{\pi}(m)},I_{\hat{\pi}(m)}}\|},\\\
\Rightarrow&\sigma_{r}(X_{T_{\hat{\pi}(m)},J_{m}^{(1)}(\hat{\pi})})\geq
2.2C\times 2.1^{M-2}\sqrt{rK}\|E_{T_{\hat{\pi}(m)},I_{\hat{\pi}(m)}}\|,\\\
\Rightarrow&\sigma_{r}(X^{*}_{T_{\hat{\pi}(m)},J_{m}^{(1)}(\hat{\pi})})\geq\left(2.2C\times
2.1^{M-2}\sqrt{rK}-1\right)\|E_{T_{\hat{\pi}(m)},J_{m}^{(1)}(\hat{\pi})}\|,\end{split}$

where we have utilized $M\geq 2$ in the second line, and the last line is due
to the Weyl’s inequality. Therefore, we have
$\underline{\mathrm{SNR}}(m,\hat{\pi})\geq 2.2\sqrt{2}C\times 2.1^{M-2}-1\geq
2.3CM-1.5C-1$. Since $\overline{\mathrm{SNR}}(m,\hat{\pi})$ also satisfies the
same upper bound as $\overline{\mathrm{SNR}}(m,\pi^{*})$ in (29), we can
follow similar arguments to bounding (32) and obtain:

$\left(1+\frac{1}{\overline{\mathrm{SNR}}(m,\hat{\pi})}\right)\left(1-\frac{1}{\underline{\mathrm{SNR}}(m,\hat{\pi})}\right)^{-1}\leq
1+\frac{1}{M-1}.$ (35)

Therefore, we also have

$\begin{split}S^{*}(\hat{\pi})=&\prod_{m=2}^{M}\left(\frac{1.1\|X^{*}_{T_{\pi^{*}(m)},I_{\pi^{*}(m)}}\|}{\sigma_{r}(X^{*}_{T_{\pi^{*}(m)},J_{m}^{(1)}(\pi^{*})})}+1\right)^{-1}\\\
\geq&\prod_{m=2}^{M}\left(\frac{1.1\|X_{T_{\pi^{*}(m)},I_{\pi^{*}(m)}}\|}{\sigma_{r}(X_{T_{\pi^{*}(m)},J_{m}^{(1)}(\pi^{*})})}+1\right)^{-1}\left(1+\frac{1}{\overline{\mathrm{SNR}}(m,\hat{\pi})}\right)^{-1}\left(1-\frac{1}{\underline{\mathrm{SNR}}(m,\hat{\pi})}\right)\\\
\geq&\hat{S}(\hat{\pi})\left(1+\frac{1}{M-1}\right)^{-(M-1)}\\\
\geq&e^{-1}\hat{S}(\hat{\pi})\\\
\geq&e^{-2}\prod_{m=2}^{M}(1.1\gamma_{m}^{-1}(\pi^{*})+1)^{-1},\end{split}$

where the second line is due to the Weyl’s inequality, the third line is due
to (35), and the last two lines are due to (33). The definition of
$S^{*}(\pi)$ then immediately implies (27), and hence our proof of Theorem 1
is complete.

For the proof of Corollary 1 and Theorem 2, we notice that their proofs do not
depend on the patch ordering, and can directly follow given the proof above of
1.

##  Appendix B Additional Simulation Studies

In this section, we show additional simulation studies from the Gaussian
mixture model in Section 3.1, comparing Cluster Quilting with the same
comparison incomplete spectral clustering methods that we use in our empirical
studies in the main text: MSC-AGL (Wen et al., , 2018), IMG (Zhao et al., ,
2016), DAIMC (Hu and Chen, 2019a, ), and OPIMC (Hu and Chen, 2019b, ).
Performances are compared via the Adjusted Rand index with respect to the true
underlying cluster labels of the underlying mixture model.

###  B.1 Sequential Patch Observed Data, Oracle Tuning

We first investigate the efficacy of Cluster Quilting in the sequential patch
observed data setting. Here, we generate data using the procedure described in
Section 4.1, with one single cluster centroid matrix used for the entire data
set. We then produce the sequential patchwork observation structure by
delineating $M$ sequential blocks in the data matrix each containing $o$ rows,
with overlap between the rows included in consecutive observation blocks, and
$p/M$ non-overlapping columns. These blocks are treated as the observed
portion of the data, while the entries of the data matrix which fall outside
of the designated blocks are masked. In this step, we ensure that all mixture
clusters are represented at least once in each patch. We apply Cluster
Quilting and other incomplete spectral clustering methods to the masked data
set, with the goal of recovering the true underlying cluster labels of the
mixture model. In the simulation study results shown below, we set the
baseline parameters of the simulation to be such that the full data sets
contain 710 observations of 100 features, with each observation belonging to 1
of 3 true underlying clusters. The cluster centroid matrix is rank 2, with
$d=4.5$ used to generate the elements of the $W$ matrix. From this, 4 blocks
are created from the data, each containing 210 observations, which represent
the simulated observed portion of the data. We then alter different elements
of the simulation parameters in order to study how the performance of each
method changes under different conditions; in particular, we present results
for changing the number of observed blocks, the number of observations in each
block, the number of total features in the data, the number of underlying
clusters, the distance between clusters, and the correlation between features.
For each simulation setting, we perform 50 repetitions and show the
aveestimated cluster labels returned by the Cluster Quilting procedure.
Columns are orderedrage and standard deviation of the adjusted Rand index.

(a) Number of blocks.

(b) Size of blocks.

(c) Number of total features.

(d) Number of clusters.

(e) Cluster distance.

(f) Feature autocorrelation.

Figure 6: Performance of Cluster Quilting and comparison cluster imputation
methods with oracle tuning on Gaussian mixture model data in the sequential
patch setting for different simulation parameters, measured by adjusted Rand
index.

The results of applying Cluster Quilting and the comparative cluster
imputation methods are shown in Figure 6. In general, we see that the Cluster
Quilting algorithm outperforms the other competing clustering imputation
methods across all of the different variations of the simulation parameters,
both when using oracle as well as data-driven hyperparameter tuning.
Clustering on the similarity matrix with missing values imputed by the block
singular value decomposition method yields fairly accurate results as well,
close to the performance level of Cluster Quilting, while the other comparison
methods do not recover the underlying clusters well in any scenario. Across
the different sets of simulation parameters, we observe several expected
trends in terms of the impact on accuracy. When the size of each observation
block is decreased, either through directly decreasing the number of rows per
block (Figure 6(a)) or increasing the number of blocks while holding the
number of blocks multiplied by the number of rows per block constant (Figure
6(b)). On the other hand, as the number of features is increased, all methods
perform better (Figure 6(c)); this can likely be attributed to having more
information from having more data to use for estimation, including an increase
in the size of the overlap between each pair of observation blocks. When we
increase the number of clusters (Figure 6(d)) or decrease the distance between
clusters (Figure 6(e)) as measured by the maximum difference between feature
means across different clusters, we see that all methods become less accurate,
which matches what we would expect as this should make the clustering problem
more difficult in general. Also, we find that accuracy of all methods is
relatively consistent across all levels of autocorrelation (Figure 6(f)). In
Section 2 in the Supporting Information, we include additional simulation
results for patchwork observed data from a Gaussian mixture model for the
high-dimensional case, as well as a non-Gaussian simulation study.

###  B.2 Data Driven Tuning

We next present results for the same simulated data as was used for our
empirical studies in Sections 3.2 in the main text for the mosaic patch
observation setting and in Section for the sequential patch observation
setting, but with data-driven tuning of the rank of the data matrix and the
number of clusters. Here, hyperparameter selection is done using the
prediction validation method of Tibshirani and Walther,  (2005), as outlined
in Section 2.1 of the main text. We see that relative performance of all
methods is consistent with the oracle tuning results, i.e. the Cluster
Quilting method performs better than the comparison incomplete spectral
clustering methods in most simulation settings. Additionally, in general,
adjusted Rand Index scores across all methods decreases when compared to
performance, showing that rank and cluster misspecification has a modest
effect on cluster estimation accuracy.

(a) Changing number of blocks.

(b) Changing size of blocks.

(c) Changing number of total features.

(d) Changing number of clusters.

(e) Changing cluster distance.

(f) Autocorrelation between features.

Figure 7: Performance of Cluster Quilting and comparison cluster imputation
methods with data-driven tuning on Gaussian mixture model data in the
sequential patch setting, measured by adjusted Rand index.

(a) Changing views observed per block.

(b) Changing number of blocks.

(c) Changing number of total features.

(d) Changing number of true clusters.

(e) Changing cluster distance.

(f) Autocorrelation between features.

Figure 8: Performance of Cluster Quilting and comparison cluster imputation
methods with data-driven tuning on Gaussian mixture model data in the mosaic
patch setting, measured by adjusted Rand index.

###  B.3 High Dimensional Studies

Here, we provide additional simulation studies in a high-dimensional setting.
For the sequential patch observation setting, the simulated data contain 710
observations spread across three true underlying clusters; 4 synthetic
observation blocks of 210 observations are then created from the data. In
Figures 9(a) and 9(b) we present results when changing the number of features
in the data between 1000 and 2000, as opposed to the 100 features in the
simulations in the main text. As we would expect from the pattern from the
main text simulations, we find that the increasing the number of features is
correlated with increased accuracy of cluster estimates, although there seems
to be an upper limit of the performance overall at an adjusted Rand Index of
around 0.7.

For the mosaic patch observation setting in the high-dimensional setting, we
generate a full data matrix of 840 observations split into 3 underlying
cluster, and 12 views of 50 features each. The cluster centroid matrix used to
generate the data for each view is a rank 2 matrix with $d=7.5$. We then
create 4 observation blocks, for which 6 views per observation block are
selected to be part of the set of observed entries. Figures 9(c) and 9(d) show
the results for varying the number of features per view between 100 and 500.
Overall, we see that the Cluster Quilting method performs better than any of
the other incomplete spectral clustering methods in the simulation study
across the varying simulation parameters.

(a) Sequential patch data,
oracle tuning.

(b) Sequential patch data,
data-driven tuning.

(c) General patchwork data,
oracle tuning.

(d) General patchwork data,
data-driven tuning.

Figure 9: Performance of Cluster Quilting and comparison cluster imputation
methods on Gaussian mixture model data with different missingness patterns in
a high-dimensional setting, measured by Adjusted Rand index with respect to
the underlying cluster membership.

###  B.4 t-Distribution Mixture Model, Sequential Patch Observation

We now test the Cluster Quilting and comparison methods in the sequential
patch observation setting on a mixture model with a heavy tailed non-sub-
Gaussian distribution. To do this, we generate data from the Gaussian mixture
model as described in Section 5a of the main text, then apply a copula
transform to each univariate marginal distribution to a t-distribution with 5
degrees of freedom. We then create the sequential patch observation pattern by
partitioning the data matrix into $M$ blocks with $o$ rows, with overlap
between the rows included in consecutive observation blocks, and $p/M$ non-
overlapping columns. We then estimate clusters from partially observed data
matrix and measure performance against the objective of recovering the cluster
labels used to generate the original Gaussian mixture model. As baseline
simulation parameters, we set mixture model to have 710 observations of 100
features in 3 clusters. 4 synthetic observation blocks are then created from
the data, each containing 210 observations. As in the main text, we study the
effects of changing the number of observed blocks, the number of observations
in each block, the number of features, the number of underlying clusters, and
the distance between clusters.

The simulation results are shown in Figure 10. Similar to the Gaussian mixture
model simulations, the Cluster Quilting method demonstrates superior
performance across almost all of the different simulation parameter settings,
while most of the other methods tend to be substantially less accurate. In
terms of the effects of varying simulation parameters on performance, we also
observe similar patterns as in the Gaussian mixture model results, in which
increasing the number of features, size of blocks, and distance between
clusters improves the accuracy of all methods, while increasing the number of
blocks between which the observations are divided negatively impacts recovery
of cluster labels. One interesting deviation from previous results is that,
for the t-distribution mixture model, increasing the number of clusters while
holding the number of total rows in the data constant often increases the
accuracy of all methods. This may be due to the heavy-tailed shape of the
t-distribution, which increases the probability of an observations in one
cluster to be incorrectly classified empirically due to being closer to points
in a different cluster which does not match the true label.

(a) Changing number of blocks;
oracle tuning.

(b) Changing number of blocks;
data-driven tuning.

(c) Changing size of observed blocks;
oracle tuning.

(d) Changing size of observed blocks;
data-driven tuning.

(e) Changing number of total features;
oracle tuning.

(f) Changing number of total features;
data-driven tuning.

(g) Changing number of underlying clusters;
oracle tuning.

(h) Changing number of underlying clusters;
data-driven tuning.

(i) Changing cluster distance;
oracle tuning.

(j) Changing cluster distance;
data-driven tuning.

Figure 10: Performance of Cluster Quilting and comparison cluster imputation
methods on sequential patch observe data from a multivariate t-distribution
with 5 degrees of freedom, measured by Adjusted Rand index with respect to the
underlying cluster membership.

###  B.5 Copula Mixture Model, Mosaic Patch Observation

We also evaluate Cluster Quilting and other incomplete spectral clustering
methods in the mosaic patch observation setting for data generated from a
copula mixture model in which we set each view to follow a different
multivariate distribution. To do this, we generate data from Gaussian mixture
models from each view as described in Section 5b of the main text, then apply
a copula transform to each column. We then create the mosaic patch observation
pattern by randomly segmenting the rows of the full data matrix into $M$
blocks of observations and then randomly selecting $h$ views for each
observation block to leave in the set of observed entries. We use default
parameters of 840 observations, 3 underlying clusters, 12 views, and 50
features per view. The cluster centroid matrix is a rank 2 matrix with
$d=7.5$. After data generation, we randomly divide the rows into 4 blocks of
equal size and select 6 views to be part of the set of observed entries. For
this simulation, each view is randomly selected to follow either a skew
normal, Gaussian, generalized normal, or beta distribution, with randomized
shape parameters; all columns within a view follow the same distribution.

We show the results of the simulation study in Figure 11. In general, the
Cluster Quilting method is able to more accurately recover the true underlying
cluster labels used to generate the copula mixture model compared to the other
incomplete spectral clustering methods for both oracle and data-driven
hyperparameter tuning. The notable exceptions to this are when the number of
observation blocks is relatively large, which follows what is presented in the
theoretical results of the main text, and when the proportion of the data
matrix that is observed is relatively small, when other methods tend to
achieve better cluster recovery. When varying the simulation parameters, we
find that all methods tend to perform better when the cluster centroids are
further apart and when more data is observed; this matches the results from
the Gaussian mixture model simulations in Section 5.

(a) Changing views observed per block;
oracle tuning.

(b) Changing views observed per block;
data-driven tuning.

(c) Changing number of blocks;
oracle tuning.

(d) Changing number of blocks;
data-driven tuning.

(e) Changing number of total features;
oracle tuning.

(f) Changing number of total features;
data-driven tuning.

(g) Changing number of true clusters;
oracle tuning.

(h) Changing number of true clusters;
data-driven tuning.

(i) Changing cluster distance;
oracle tuning.

(j) Changing cluster distance;
data-driven tuning.

Figure 11: Performance of Cluster Quilting and comparison cluster imputation
methods on copula mixture model data in the mosaic patch observation setting,
measured by adjusted Rand index.

##  Appendix C Calcium Imaging Case Study

Here, we apply Cluster Quilting to the problem of estimating functional
neuronal clusters, i.e. how individual neurons can be organized into broad
groups based on correlations in firing activity, using data collected via two-
photon calcium imaging. Due to current technological limitations, only a small
fraction of a brain region can be recorded with enough granularity to
differentiate individual neurons. Neural activity data from calcium imaging
for entire experiment is therefore generally collected via a series of scans
at sequentially increasing depths, with some overlap between the neurons
recorded between consecutive scans. This leads to the sequential patch
observation scheme as described in Section 1, and thus requires the
implementation of an incomplete data clustering technique in order to estimate
function clustering of neurons for a full brain region of interest. In this
section, we study the applicability of the Cluster Quilting algorithm to
calcium imaging data. Below, we analyze a data set from the Allen Institute
containing fluorescence traces for 227 neurons measured across approximately
8000 observations during spontaneous activity. For this case study, we first
apply spectral clustering to the complete, unmasked data set, with data-driven
tuning of the rank and number of clusters using the prediction validation
method described in Section 2; we treat this as the ground truth clustering to
which we compare the results of Cluster Quilting and competing methods. We
then simulate the missingness pattern found in calcium imaging data collected
in multiple patches by designating synthetic sequential patch observation
blocks with overlap in the rows but not the columns as part of the empirical
observed set of entries in the data matrix, and we use this as the input to
the Cluster Quilting procedure and the comparison incomplete spectral
clustering methods. We compare the procedures based on how well they recover
the cluster assignments estimated by spectral clustering on the fully observed
data for varying block sizes and for different numbers of sequential blocks
while holding the total size of the set of observed entries in the data matrix
constant. For each combination of the number of blocks and block size, we
perform 50 replications and report the mean and standard deviation of the
adjusted Rand index with respect to the cluster labels found using the fully
observed data set for each method. We consider both oracle and data-driven
tuning of the rank and the number of clusters; for this particular case study,
we treat the rank and number of clusters found via data driven tuning with
spectral clustering applied to the fully observed data set as the oracle rank
and number of clusters to use as input to the incomplete spectral clustering
methods.

(a) Changing number of blocks;
oracle tuning.

(b) Changing number of blocks;
data-driven tuning.

(c) Changing size of observed blocks;
oracle tuning.

(d) Changing size of observed blocks;
data-driven tuning.

Figure 12: Performance of Cluster Quilting and comparison cluster imputation
methods on Allen Institute data with synthetic simulated block-missingness,
measured by adjusted Rand index with respect to the cluster labels found with
the fully observed data. In general, Cluster Quilting recovers cluster labels
that most closely match those from the full data compared to other methods.

The performance of Cluster Quilting and the comparison incomplete spectral
clustering methods in recovering the cluster labels found with the fully
observed data in Figure 12, with oracle hyperparameter tuning results in the
left column and data-driven hyperparameter tuning in the right column. Across
the varying block sizes and number of blocks, we find that the Cluster
Quilting is the best performing method in terms of recovering the same
clustering assignments that are assigned by applying spectral clustering to
the fully observed data. Cluster labels from the imputed similarity matrix
from the block singular value decomposition also matches relatively well on
average with those estimated from the fully observed data. On the other hand,
cluster assignments from applying spectral clustering after zero imputation of
the incomplete similarity matrix is by far the worst performing method. These
results show that Cluster Quilting may be the most appropriate method in order
to estimate functional clusters in neurons recorded across multiple calcium
imaging scans. As in the simulation studies in the main text, we see that
increasing both the proportion of the data matrix that is observed and the
size of the intersections between each pair of consecutive blocks, achieved by
increasing the number of neurons in each observation block (Figures 12(a) and
12(b)) or decreasing the total number of observation blocks amongst which the
neurons are divided (Figures 12(c) and 12(d)), increases the concordance
between the clusters predicted by the imputation methods and the cluster
labels found from applying spectral clustering on the full data.

(a) Full data.

(b) DAIMC.

(c) Likelihood data imputation.

(d) Cluster Quilting.

Figure 13: Spatial distribution of functional clusters estimated from neuronal
activity data collected via two-photon calcium imaging. We find that the
spatial patterns of the functional clusters found via Cluster Quilting (d) are
most closely aligned with those derived from the full data (a), as compared to
the functional clusters estimated with other incomplete spectral clustering
techniques (b) and nuclear-norm-regularized-likelihood imputation (c) on the
missing data. Calcium imaging provided by the Allen Institute.

We also analyze the spatial distribution of the functional neuronal clusters
found with full data to those estimated via Cluster Quilting, alongside the
estimates from the data with missing values imputed using nuclear-norm
regularization in Figure 13. In order to visualize the spatial patterns of the
different clustering estimates, we show the cluster labels with respect to the
location of the neurons in the brain slice. Figure 13(a) shows the functional
clusters estimated on the full data set. Here, we see that the estimated
functional clusters appear to have a relationship with the spatial location of
the neurons, as the neurons in the same function clusters tend to be grouped
spatially into several distinct spatial areas; in particular, the clusters
seem to be divided spatially along the horizontal axis of the brain slice.
This finding follows what has been found previously in the study of functional
neuronal clustering from calcium imaging data (Poliakov and Schieber, , 1999;
Dombeck et al., , 2009; Wilson et al., , 2016). In Figures 13(b) and 13(c), we
show an example of the functional clusters estimated from imputation of the
masked data via the DAIMC and nuclear-norm-regularized likelihood methods,
respectively. We see in these cases that the spatial patterns that underlie
the functional clusters estimated from the full data are much less apparent in
the functional clustering assignments from the data imputed using the low-rank
factorization method and zero imputation, instead having much less structured
spatial distributions. On the other hand, in Figure 13(d), which shows the
functional clusters identified by applying Cluster Quilting to the same masked
data, we see more separable spatial patterns as well as cluster assignments
that more closely match the cluster labels found from fully observed data.
From this case study, we find that Cluster Quilting provides imputed
functional neuronal clustering estimates from calcium imaging data that better
match what would be found if the full data were to be observed and which
correspond more closely with what we would expect in the specific neuroscience
context. specifically with respect to the relationship to the spatial
distributions of the functional clusters.

## References

  * Abbe et al., (2022) Abbe, E., Fan, J., and Wang, K. (2022).  An $\ell_{p}$ theory of pca and spectral clustering.  The Annals of Statistics, 50(4):2359–2385.
  * Aggarwal and Zhai, (2012) Aggarwal, C. and Zhai, C. (2012).  A survey of text clustering algorithms.  Springer.
  * Aktaş et al., (2019) Aktaş, M., Kaplan, S., Abacı, H., Kalipsiz, O., Ketenci, U., and Turgut, U. (2019).  Data imputation methods for missing values in the context of clustering.  In Big Data and Knowledge Sharing in Virtual Organizations, pages 240–274. IGI Global.
  * Baghfalaki et al., (2016) Baghfalaki, T., Ganjali, M., and Berridge, D. (2016).  Missing value imputation for rna-sequencing data using statistical models: a comparative study.  Journal of Statistical Theory and Applications, 15(3):221–236.
  * Bishop and Yu, (2014) Bishop, W. E. and Yu, B. M. (2014).  Deterministic symmetric positive semidefinite matrix completion.  Advances in Neural Information Processing Systems, 27.
  * Boluki et al., (2019) Boluki, S., Zamani Dadaneh, S., Qian, X., and Dougherty, E. R. (2019).  Optimal clustering with missing values.  BMC bioinformatics, 20:1–10.
  * Bridges, (1966) Bridges, C. (1966).  Hierarchical cluster analysis.  Psychological Reports, 18(3):851–854.
  * Browning, (2008) Browning, S. (2008).  Missing data imputation and haplotype phase inference for genome-wide association studies.  Human Genetics, 124:439–450.
  * Cai et al., (2019) Cai, Q., Wang, H., Li, Z., and Liu, X. (2019).  A survey on multimodal data-driven smart healthcare systems: approaches and applications.  IEEE Access, 7:133583–133599.
  * Chang et al., (2022) Chang, A., Zheng, L., and Allen, G. (2022).  Low-rank covariance completion for graph quilting with applications to functional connectivity.  arXiv preprint arXiv:2209.08273.
  * Chatterjee, (2020) Chatterjee, S. (2020).  A deterministic theory of low rank matrix completion.  IEEE Transactions on Information Theory, 66(12):8046–8055.
  * Chen et al., (2015) Chen, G., Chi, E., Ranola, J., and Lange, K. (2015).  Convex clustering: An attractive alternative to hierarchical clustering.  PLoS Computational Biology, 11(5):e1004228.
  * Chen et al., (2021) Chen, Y., Chi, Y., Fan, J., Ma, C., et al. (2021).  Spectral methods for data science: A statistical perspective.  Foundations and Trends® in Machine Learning, 14(5):566–806.
  * Chi et al., (2016) Chi, J. T., Chi, E. C., and Baraniuk, R. G. (2016).  k-pod: A method for k-means clustering of missing data.  The American Statistician, 70(1):91–99.
  * Choi and Yuan, (2023) Choi, J. and Yuan, M. (2023).  Matrix completion when missing is not at random and its applications in causal panel data models.  arXiv preprint arXiv:2308.02364.
  * Dalton et al., (2009) Dalton, L., Ballarin, V., and Brun, M. (2009).  Clustering algorithms: on learning, validation, performance, and applications to genomics.  Current Genomics, 10(6):430–445.
  * Dombeck et al., (2009) Dombeck, D. A., Graziano, M. S., and Tank, D. W. (2009).  Functional clustering of neurons in motor cortex determined by cellular resolution imaging in awake behaving mice.  Journal of Neuroscience, 29(44):13751–13760.
  * Donath and Hoffman, (1972) Donath, W. and Hoffman, A. (1972).  Algorithms for partitioning of graphs and computer logic based on eigenvectors of connection matrices.  IBM Technical Disclosure Bulletin, 15(3):938–944.
  * Gemma et al., (2001) Gemma, A., Hosoya, Y., Seike, M., Uematsu, K., Kurimoto, F., Hibino, S., Yoshimura, A., Shibuya, M., Kudoh, S., and Emi, M. (2001).  Genomic structure of the human mad2 gene and mutation analysis in human lung and breast cancers.  Lung Cancer, 32(3):289–295.
  * Hathaway and Bezdek, (2001) Hathaway, R. J. and Bezdek, J. C. (2001).  Fuzzy c-means clustering of incomplete data.  IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics), 31(5):735–744.
  * (21) Hu, M. and Chen, S. (2019a).  Doubly aligned incomplete multi-view clustering.  arXiv preprint arXiv:1903.02785.
  * (22) Hu, M. and Chen, S. (2019b).  One-pass incomplete multi-view clustering.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 33, pages 3838–3845.
  * Hunt and Jorgensen, (2003) Hunt, L. and Jorgensen, M. (2003).  Mixture model clustering for mixed data with missing information.  Computational statistics & data analysis, 41(3-4):429–440.
  * Inman et al., (2015) Inman, D., Elmore, R., and Bush, B. (2015).  A case study to examine the imputation of missing data to improve clustering analysis of building electrical demand.  Building Services Engineering Research and Technology, 36(5):628–637.
  * Jin, (2015) Jin, J. (2015).  Fast community detection by score.  THE ANNALS, 43(1):57–89.
  * Jin et al., (2018) Jin, L., Han, B., Siegel, E., Cui, Y., Giuliano, A., and Cui, X. (2018).  Breast cancer lung metastasis: Molecular biology and therapeutic implications.  Cancer biology & therapy, 19(10):858–868.
  * Kline et al., (2022) Kline, A., Wang, H., Li, Y., Dennis, S., Hutch, M., Xu, Z., Wang, F., Cheng, F., and Luo, Y. (2022).  Multimodal machine learning in precision health: A scoping review.  npj Digital Medicine, 5(1):171.
  * Landemaine et al., (2008) Landemaine, T., Jackson, A., Bellahcene, A., Rucci, N., Sin, S., Abad, B. M., Sierra, A., Boudinet, A., Guinebretiere, J.-M., Ricevuto, E., and Nogues, C. (2008).  A six-gene signature predicting breast cancer lung metastasis.  Cancer research, 68(15):6092–6099.
  * Lei and Rinaldo, (2015) Lei, J. and Rinaldo, A. (2015).  Consistency of spectral clustering in stochastic block models.  The Annals of Statistics, 43(1):215–237.
  * Lin et al., (2006) Lin, T. I., Lee, J. C., and Ho, H. J. (2006).  On fast supervised learning for normal mixture models with missing information.  Pattern Recognition, 39(6):1177–1187.
  * Lock and Dunson, (2013) Lock, E. F. and Dunson, D. B. (2013).  Bayesian consensus clustering.  Bioinformatics, 29(20):2610–2616.
  * Löffler et al., (2021) Löffler, M., Zhang, A. Y., and Zhou, H. H. (2021).  Optimality of spectral clustering in the gaussian mixture model.  The Annals of Statistics, 49(5):2506–2530.
  * MacQueen, (1967) MacQueen, J. (1967).  Some methods for classification and analysis of multivariate observations.  In Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability, Volume 1, pages 281–297. University of California Press.
  * MICrONS Consortium, (2021) MICrONS Consortium (2021).  Functional connectomics spanning multiple areas of mouse visual cortex.  BioRxiv, pages 2021–07.
  * Mishra et al., (2007) Mishra, N., Schreiber, R., Stanton, I., and Tarjan, R. (2007).  Clustering social networks.  In International Workshop on Algorithms and Models for the Web-Graph, pages 56–67, Berlin, Heidelberg. Springer Berlin Heidelberg.
  * Montel et al., (2005) Montel, V., Huang, T.-Y., Mose, E., Pestonjamasp, K., and Tarin, D. (2005).  Expression profiling of primary tumors and matched lymphatic and lung metastases in a xenogeneic breast cancer model.  The American journal of pathology, 166(5):1565–1579.
  * Poddar and Jacob, (2019) Poddar, S. and Jacob, M. (2019).  Clustering of data with missing entries using non-convex fusion penalties.  IEEE transactions on signal processing, 67(22):5865–5880.
  * Poliakov and Schieber, (1999) Poliakov, A. V. and Schieber, M. H. (1999).  Limited functional grouping of neurons in the motor cortex hand area during individuated finger movements: a cluster analysis.  Journal of Neurophysiology, 82(6):3488–3505.
  * Rajendran et al., (2023) Rajendran, S., Pan, W., Sabuncu, M. R., Zhou, J., and Wang, F. (2023).  Patchwork learning: A paradigm towards integrative analysis across diverse biomedical data sources.  arXiv preprint arXiv:2305.06217.
  * Rohe et al., (2011) Rohe, K., Chatterjee, S., and Yu, B. (2011).  Spectral clustering and the high-dimensional stochastic blockmodel.  Annals of statistics, 39(4):1878–1915.
  * Santos and Embrechts, (2009) Santos, J. and Embrechts, M. (2009).  On the use of the adjusted rand index as a metric for evaluating supervised classification.  In International Conference on Artificial Neural Networks, pages 175–184, Berlin, Heidelberg. Springer Berlin Heidelberg.
  * Schwartz et al., (1999) Schwartz, A. G., Siegfried, J. M., and Weiss, L. (1999).  Familial aggregation of breast cancer with early onset lung cancer.  Genetic Epidemiology: The Official Publication of the International Genetic Epidemiology Society, 17(4):274–284.
  * Shen et al., (2009) Shen, R., Olshen, A., and Ladanyi, M. (2009).  Integrative clustering of multiple genomic data types using a joint latent variable model with application to breast and lung cancer subtype analysis.  Bioinformatics, 25(22):2906–2912.
  * Tibshirani and Walther, (2005) Tibshirani, R. and Walther, G. (2005).  Cluster validation by prediction strength.  Journal of Computational and Graphical Statistics, 14(3):511–528.
  * Vinci et al., (2019) Vinci, G., Dasarathy, G., and Allen, G. I. (2019).  Graph quilting: graphical model selection from partially observed covariances.  arXiv preprint arXiv:1912.05573.
  * Wang and Allen, (2021) Wang, M. and Allen, G. I. (2021).  Integrative generalized convex clustering optimization and feature selection for mixed multi-view data.  Journal of Machine Learning Research, 22(55):1–73.
  * Wang et al., (2019) Wang, S., Li, M., Hu, N., Zhu, E., Hu, J., Liu, X., and Yin, J. (2019).  K-means clustering with incomplete data.  IEEE Access, 7:69162–69171.
  * Weinstein, J.N. et. al., (2013) Weinstein, J.N. et. al. (2013).  The cancer genome atlas pan-cancer analysis project.  Nature Genetics, 45(10):1113–1120.
  * Wen et al., (2018) Wen, J., Xu, Y., and Liu, H. (2018).  Incomplete multiview spectral clustering with adaptive graph learning.  IEEE Transactions on Cybernetics, 50(4):1418–1429.
  * Wen et al., (2022) Wen, J., Zhang, Z., Fei, L., Zhang, B., Xu, Y., Zhang, Z., and Li, J. (2022).  A survey on incomplete multiview clustering.  IEEE Transactions on Systems, Man, and Cybernetics: Systems, 53(2):1136–1149.
  * Wilson et al., (2016) Wilson, D. E., Whitney, D. E., Scholl, B., and Fitzpatrick, D. (2016).  Orientation selectivity and the functional clustering of synaptic inputs in primary visual cortex.  Nature neuroscience, 19(8):1003–1009.
  * Yan and Wainwright, (2024) Yan, Y. and Wainwright, M. J. (2024).  Entrywise inference for causal panel data: A simple and instance-optimal approach.  arXiv preprint arXiv:2401.13665.
  * Yu et al., (2020) Yu, G., Li, Q., Shen, D., and Liu, Y. (2020).  Optimal sparse linear prediction for block-missing multi-modality data without imputation.  Journal of the American Statistical Association, 115(531):1406–1419.
  * Zhao et al., (2016) Zhao, H., Liu, H., and Fu, Y. (2016).  Incomplete multi-modal visual data grouping.  In IJCAI, pages 2392–2398.
  * Zhou and Chen, (2023) Zhou, Y. and Chen, Y. (2023).  Heteroskedastic tensor clustering.  arXiv preprint arXiv:2311.02306, 3.
  * Zhou et al., (2022) Zhou, Y., Zhang, A. R., Zheng, L., and Wang, Y. (2022).  Optimal high-order tensor svd via tensor-train orthogonal iteration.  IEEE transactions on information theory, 68(6):3991–4019.