# Statistically Optimal $K$-means Clustering via Nonnegative Low-rank
Semidefinite Programming

Yubo Zhuang*, Xiaohui Chen*†, Yun Yang*, Richard Y. Zhang*
University of Illinois at Urbana-Champaign*
University of Southern California†
{yubo2, xhchen, yy84<EMAIL_ADDRESS>

###### Abstract

$K$-means clustering is a widely used machine learning method for identifying
patterns in large datasets. Semidefinite programming (SDP) relaxations have
recently been proposed for solving the $K$-means optimization problem that
enjoy strong statistical optimality guarantees, but the prohibitive cost of
implementing an SDP solver renders these guarantees inaccessible to practical
datasets. By contrast, nonnegative matrix factorization (NMF) is a simple
clustering algorithm that is widely used by machine learning practitioners,
but without a solid statistical underpinning nor rigorous guarantees. In this
paper, we describe an NMF-like algorithm that works by solving a _nonnegative_
low-rank restriction of the SDP relaxed $K$-means formulation using a
nonconvex Burer–Monteiro factorization approach. The resulting algorithm is
just as simple and scalable as state-of-the-art NMF algorithms, while also
enjoying the same strong statistical optimality guarantees as the SDP. In our
experiments, we observe that our algorithm achieves substantially smaller mis-
clustering errors compared to the existing state-of-the-art.

##  1 Introduction

Clustering remains a common unsupervised learning technique, where its basic
objective is to assign similar data points to the same group. For data in the
Euclidean space, a widely used clustering method is the $K$-means clustering
that computes similarity in terms of the distances from data to the learned
clustering centers (i.e., centroids) (MacQueen, 1967). Given data points
$X_{1},\dots,X_{n}\in\mathbb{R}^{p}$ to be divided into $K$ groups, $K$-means
clustering aims to minimize the following cost function:

$\min_{\beta_{1},\dots,\beta_{K}\in\mathbb{R}^{p}}\sum_{i=1}^{n}\min_{k\in[K]}\|X_{i}-\beta_{k}\|_{2}^{2},$
(1)

where $\beta_{k}$ is the centroid of the $k$-th cluster for
$k\in[K]:=\\{1,\dots,K\\}$. Solving the $K$-means problem in (1) is equivalent
to solving mixed-integer linear program (MILP) that is known to be worst-case
NP hard (Dasgupta, 2007; Aloise et al., 2009). Various computationally
tractable approximation algorithms and relaxed formulations have been
extensively studied in literature. Notable examples include Lloyd’s algorithm
(Lloyd, 1982), spectral clustering (von Luxburg, 2007; Ng et al., 2001),
nonnegative matrix factorization (NMF) (He et al., 2011; Kuang et al., 2015;
Wang and Zhang, 2012), and semidefinite programming (SDP) (Peng and Wei, 2007;
Mixon et al., 2017; Royer, 2017; Fei and Chen, 2018; Giraud and Verzelen,
2018).

Among those popular relaxations, the SDP approach enjoys the strongest
statistical guarantees under the standard Gaussian mixture model in that it
achieves an information-theoretic sharp threshold for the purpose for exact
recovery of the true cluster partition (Chen and Yang, 2021). Unfortunately,
the SDP and its strong statistical guarantees remain completely inaccessible
to real-world datasets, owing to the prohibitively high costs of solving the
resulting SDP relaxation. Given $n$ data points, the SDP is a matrix
optimization problem, over a dense $n\times n$ membership matrix, that is
constrained to be both positive semidefinite as well as elementwise
nonnegative. Even ignoring the constraints, a basic but fundamental difficulty
is the need to store and optimize over the $n^{2}$ individual elements of the
matrix. Even a small dataset like the banknote authentication dataset from the
UCI collection (Dua and Graff, 2017) with $n\approx 1000$ translates into an
SDP with $n\approx 10^{6}$ optimization variables, which is right at the very
limit of state-of-the-art SDP solvers like SDPNAL+ (Yang et al., 2015).

On the other hand, NMF remains one of the simplest and practically useful
approaches to clustering due to its scalability (He et al., 2011; Kuang et
al., 2015). When the problem at hand exhibits latent low-dimensional
structures (e.g., clustering problem is naturally of low-rank), the
computational savings are significant because it is able to impose elementwise
nonnegativity over the $n\times n$ membership matrix using just $O(n)$
explicit constraints imposed on the factor matrix. While highly scalable,
there remains unfortunately very little statistical underpinning behind NMF-
based algorithms.

Figure 1: Log-scale trade-off plot of CPU time versus mis-clustering error,
under an increasing number of data points $n$. Here, BM corresponds to our
proposed method, KM corresponds to $K$-means++ (Arthur and Vassilvitskii,
2007), SC corresponds to spectral clustering (Ng et al., 2001), SDP (Peng and
Wei, 2007) uses the SDPNAL+ solver (Yang et al., 2015), and NMF corresponds to
non-negative factorization (Ding et al., 2005). We follow the same setting as
the first experiment in Section 5, where the theoretically optimal mis-
clustering error decays to zero as $n$ increases to infinity.

Our contributions. In this paper, we propose an efficient, large-scale, NMF-
like algorithm for the $K$-means clustering problem, that _meanwhile_ enjoys
the same sharp exact recovery guarantees provided by SDP relaxations. Our
motivating insight is that SDP and NMF (as well as spectral clustering) solve
the same $K$-means objective function subject to different constraint sets.
This gives us the hope to break the computational and statistical bottleneck
via simultaneously leveraging the implicit positive-semidefinite (psd)
structure of the membership matrix and nonnegativity constraint.

At its core, our proposed algorithm is a primal-dual gradient descent-ascent
algorithm to optimize over the nonnegative factor matrix with respect to an
augmented Lagrangian of the SDP. The resulting iterations closely resemble the
projected gradient descent algorithms widely used for NMF in the existing
literature; in fact, we show that the latter can be recovered from our
algorithm by relaxing two constraints. We prove that the new algorithm enjoys
primal-dual local linear convergence within a neighborhood of the SDP
solution, which is a unique rank-$K$ block diagonal matrix whenever the
centroids satisfy a well-separation condition (Chen and Yang, 2021). In
practice, we observe that the algorithm converges globally at a linear rate.
As shown in Figure 1, our algorithm achieves substantially smaller mis-
clustering errors compared to the existing state-of-the-art.

Organization. The rest of the paper is organized as follows. In Section 2, we
present some background on several equivalent and relaxed formulations of the
$K$-means clustering problem. In Section 3, we introduce a nonnegative low-
rank SDP for solving the Burer–Monteiro formulation of $K$-means problem. In
Section 4, we establish the linear convergence guarantee for our proposed
primal-dual gradient descent-ascent algorithm in the exact recovery regime.
Numerical experiments are reported in Section 5. Proof details are deferred to
the Supplementary Material.

##  2 Background on $K$-means and Related Computational Methods

Due to the parallelogram law in Euclidean space, the centroid-based
formulation of the $K$-means clustering problem in (1) has an equivalent
_partition-based_ formulation (cf. Zhuang et al. (2022)) as

$\max_{G_{1},\dots,G_{K}}\Big{\\{}\sum_{k=1}^{K}{1\over|G_{k}|}\sum_{i,j\in
G_{k}}\langle X_{i},X_{j}\rangle:\bigsqcup_{k=1}^{K}G_{k}=[n]\Big{\\}},$ (2)

where the Euclidean inner product $\langle X_{i},X_{j}\rangle=X_{i}^{T}X_{j}$
represents the similarity between the vectors $X_{i}$ and $X_{j}$, the
clusters $(G_{k})_{k=1}^{K}$ form a partition of the data index $[n]$,
$|G_{k}|$ denotes the cardinality of $G_{k}$ and $\sqcup$ denotes disjoint
union. The objective function in (2) is the log-likelihood function of the
cluster labels (modulo constant) by profiling out the centroids in (1) as
nuisance parameters under the standard Gaussian mixture model with a common
isotropic covariance matrix across all the $K$ components (Zhuang et al.,
2023). Using the one-shot encoding of the partition $G_{1},\dots,G_{K}$, we
can uniquely associate it (up to label permutation) with a binary _assignment
matrix_ $H=(h_{ik})\in\\{0,1\\}^{n\times K}$ such that $h_{ik}=1$ if $i\in
G_{k}$ and $h_{ik}=0$ otherwise. Then the $K$-means clustering problem in (2)
can be written as an MILP as follows:

$\min_{H}\Big{\\{}\langle A,HBH^{T}\rangle:H\in\\{0,1\\}^{n\times
K},\;H\mathbf{1}_{k}=\mathbf{1}_{n}\Big{\\}},$ (3)

where $A=-X^{T}X$ is the $n\times n$ negative Gram matrix of the data
$X_{p\times n}=(X_{1},\dots,X_{n})$,
$B=\text{diag}(|G_{1}|^{-1},\dots,|G_{K}|^{-1})$ is a normalization matrix and
the constraint $H\mathbf{1}_{k}=\mathbf{1}_{n}$ with $\mathbf{1}_{n}$ being
the $n$-dimensional vector of all ones reflects the _row sum constraint_ of
assignment, i.e., each row of $H$ contains exactly nonzero entry with value
one.

It is known that the 0-1 integer constraint makes the MILP problem in (3) NP-
hard (Dasgupta, 2007; Aloise et al., 2009), and various relaxations for the
$K$-means constraint are formulated in literature to turn problem (3)
tractable with the same objective function. Popular methods include: (i)
spectral clustering (von Luxburg, 2007; Ng et al., 2001) which only maintains
the weakened orthonormal constraint $\tilde{H}^{T}\tilde{H}=I_{K}$ for
$\tilde{H}:=HB^{1/2}$; (ii) nonnegative matrix factorization (NMF) (He et al.,
2011; Kuang et al., 2015) which only enforces the elementwise nonnegativity
constraint on $\tilde{H}\geqslant 0$; (iii) semidefinite programming (SDP)
(Peng and Wei, 2007; Royer, 2017; Giraud and Verzelen, 2018) which
reparameterizes the assignment matrix $H$ as the psd membership matrix
$Z:=HBH^{T}\succeq 0$ and additionally preserves the constraints
$\operatorname{tr}(Z)=K$, $Z\mathbf{1}_{n}=\mathbf{1}_{n}$ and $Z\geqslant 0$,
namely we solve

$\min_{Z\in\mathbb{S}^{n}_{+}}\Big{\\{}\langle A,Z\rangle:\mbox{\rm
tr}(Z)=K,\,Z\mathbf{1}_{n}=\mathbf{1}_{n},\,Z\geqslant 0\Big{\\}},$ (4)

where $\mathbb{S}^{n}_{+}$ stands for the convex cone of the $n\times n$ real
symmetric psd matrices.

Among those popular relaxations, the SDP approach enjoys the strongest
statistical guarantees under the Gaussian mixture model (16). It is shown by
Chen and Yang (2021) that the above SDP achieves an information-theoretic
limit for the purpose for exact recovery of the true cluster partition.
Precisely, the sharp threshold on the centroid-separation for phase transition
is given by

$\overline{\Theta}^{2}=4\sigma^{2}\left(1+\sqrt{1+{Kp\over
n\log{n}}}\right)\log{n}$ (5)

in the sense that for any $\alpha>0$ and $K=O(\log(n)/\log\log(n))$, if the
minimal centroid-separation $\Theta_{\min}:=\min_{1\leqslant j\neq k\leqslant
K}\|\mu_{j}-\mu_{k}\|_{2}\geqslant(1+\alpha)\overline{\Theta}$, then with high
probability solution of (4) is unique and it perfectly recovers the cluster
partition structure; while if $\Theta\leqslant(1-\alpha)\overline{\Theta}$,
then the maximum likelihood estimator (and thus any other estimator) fails to
exactly recover the cluster partition. In simpler words, the SDP relaxed
$K$-means optimally solves clustering problem with zero mis-clustering error
as soon as it is possible to do so, i.e., there is no relaxation gap.

We shall note that $\Theta_{\min}\to\infty$ is necessarily needed to achieve
exact recovery with high probability. When the centroid-separation diverges,
spectral clustering is shown to achieve the minimax mis-clustering error rate
$\exp(-(1+o(1))\Theta_{\min}^{2}/8)$ under the condition $p=o(n\Theta_{\min})$
(or $p=o(n\Theta_{\min}^{2})$ with extra assumptions on singular values of the
population matrix $\operatorname{\mathds{E}}[X]$) (Löffler et al., 2021). This
result implies that, in the _low-dimensional_ setting, spectral clustering
asymptotically attains the sharp threshold (5) as
$\liminf_{n\to\infty}{\Theta_{\min}^{2}\over 8\log{n}}>1$. However, it remains
unclear how spectral cluster performs in theory under the high-dimensional
regime when $n\Theta_{\min}^{2}=O(p)$ because running $K$-means on the eigen-
embedded data would also depend on the structure of the population singular
values (Han et al., 2023). Thus, in practice it is not recommended to use
spectral clustering if the data is high dimensional and one often considers
more robust alternatives such as the SDP relaxation in (4).

Despite the $K$-means clustering problem solved via SDP enjoys statistical
optimality, it optimizes over an $n\times n$ dense psd matrix to make it
completely inaccessible to practical datasets with even moderate size of a few
thousands data points. By contrast, the NMF approach using a typical
workstation hardware can handle much larger scale datasets as in the documents
and images clustering applications (Cichocki and Phan, 2009; Kim et al., 2014;
Xu et al., 2003). Nonetheless, little is known in literature about the
provable statistical guarantee for a general NMF approach. For the clustering
problem, our goal is to break the computational and statistical bottleneck by
simultaneously leverage the implicit psd structure of the membership matrix
and nonnegativity constraint. In Section 3, we propose a novel NMF-like
algorithm that achieves the statistical optimality by solving a nonnegative
low-rank restricted SDP formulation.

##  3 Our Proposal: $K$-means via Nonnegative Low-rank SDP

The Burer–Monteiro (BM) is a nonconvex approach for solving a SDP problem when
its solution is expected to be low rank (Burer and Monteiro, 2003). Standard
SDP problem with only equality constraints on the psd matrix variable
$Z\in\mathbb{S}_{+}^{n}$ can be factorized as $Z=UU^{T}$ where $U$ is an
$n\times r$ rectangular matrix with $r\ll n$. In our SDP relaxed $K$-means
formulation (4), the true membership matrix $Z^{*}$ is a block diagonal matrix
containing $K$ blocks, where each block has size $n_{k}\times n_{k}$ with all
entries equal to $n_{k}^{-1}$ and $n_{k}=|G_{k}^{*}|$ is the size of the true
cluster $G_{k}^{*}$. Thus we can exploit this low-rank structure and recast
the SDP in (4) as a nonconvex optimization problem over
$U\in\mathbb{R}^{n\times r}$:

$\min_{U\in\mathbb{R}^{n\times r}}\Big{\\{}\langle
A,UU^{T}\rangle:\|U\|_{F}^{2}=K,\,UU^{T}\mathbf{1}_{n}=\mathbf{1}_{n},\,U\geqslant
0\Big{\\}},$ (6)

where we replaced the constraint $UU^{T}\geqslant 0$ with the stronger
constraint $U\geqslant 0$ that is easier to enforce as in the NMF setting
(Kulis et al., 2007). Note that the BM reformulation (6) can be viewed as a
_restriction_ of the rank-constrained SDP formulation since $U\geqslant 0$
implies $Z=UU^{T}\geqslant 0$. Even though the BM method turns a convex SDP to
a nonconvex problem, it leverages the low-rank structure to achieve
substantial computational savings. In addition, the solution $\hat{U}$ has a
natural statistical interpretation: rows of $\hat{U}$ can be interpreted as
latent positions of $n$ data points when projected into $\mathbb{R}^{r}$, and
thus can be used to conduct other downstream data analyses in reduced
dimensions. Under this latent embedding perspective, formulation (6) can be
viewed as a _constrained PCA_ respecting the clustering structure through
constraint $\|U\|_{F}^{2}=K$.

The BM formulation (6) also has a close connection to a class of clustering
algorithms based on NMF, which instead solve the following optimization
problem

$\min_{U\in\mathbb{R}^{n\times r}}\Big{\\{}\|A+UU^{T}\|_{F}^{2}:U\geqslant
0\Big{\\}}.$ (7)

The NMF formulation (7) is advocated in a number of papers (Ding et al., 2005;
Kuang et al., 2015; Zhu et al., 2018) for data clustering in the NMF
literature due to its remarkable compuational scalability, and is equivalent
to a simpler version of (6) by dropping the two equality constraints. Our
empirical results in Section 5 and Figure 1 show that keeping these two
equality constraints are beneficial: the resulting method significantly
improves the classification accuracy while maintaining comparable
computational scalability. More importantly, we have the theoretical
certification that the resulting method from (6) achieves the information-
theoretic limit for exact recovery of the true cluster labels (see Section 4).

Algorithm 1 Primal-dual algorithm for solving BM formulation (6) of $K$-means
clustering

Input. Dissimilarity matrix $A=-XX^{T}$, clustering parameter $K\geqslant 1$,
rank parameter $r\geqslant 1$, augmentation parameter $\beta>0$, step size
$\alpha>0$.

Output. A second-order locally optimal point $U$.

Algorithm. Initialize $y=0$. Do the following:

  1. 1.

(Primal descent; Projected GD steps) Recursively run the following until
convergence: $U^{0}=U$; for $t\geqslant 0$,

$U^{t+1}=\Pi_{\Omega}\big{(}U^{t}-\alpha\,\nabla_{U^{t}}\mathcal{L}_{\beta}(U^{t},y)\big{)}\quad\text{with
}\ \Pi_{\Omega}(V)=\sqrt{K}\cdot(V)_{+}/\|(V)_{+}\|_{F},$

where
$\nabla_{U^{t}}\mathcal{L}_{\beta}(U^{t},y)=(2A+2L\cdot\text{Id}+\mathbf{1}_{n}\overline{y}^{T}+\overline{y}\mathbf{1}_{n}^{T})U^{t}$
and $\overline{y}=y+\beta(U^{t}{U^{t}}^{T}\mathbf{1}_{n}-\mathbf{1}_{n})$.
Upon convergence at iteration count $t_{0}$, set $U_{\mathrm{new}}=U^{t_{0}}$.

  2. 2.

(Dual ascent; Augmented Lagrangian step) Update dual variable via

$y_{\mathrm{new}}=y+\beta(U_{\mathrm{new}}U_{\mathrm{new}}^{T}\mathbf{1}_{n}-\mathbf{1}_{n}).$

  3. 3.

(Stopping criterion) If
$\max\\{\|U_{\mathrm{new}}-U\|_{F},\,\|U_{\mathrm{new}}U_{\mathrm{new}}^{T}\mathbf{1}_{n}-\mathbf{1}_{n}\|_{F}\\}$
falls below some tolerance threshold, then return $U_{\mathrm{new}}$.
Otherwise, set $U\leftarrow U_{\mathrm{new}}$ and $y\leftarrow
y_{\mathrm{new}}$, and repeat Steps 1 and 2.

We will now derive a simple primal-dual gradient descent-ascent algorithm to
solve (6). To begin, note that we can first turn the nonsmooth inequality
constraint $U\geqslant 0$ together with the trace constraint to a manifold-
like subset

$\Omega:=\\{U\in\mathbb{R}^{n\times r}:\|U\|_{F}^{2}=K,\;U\geqslant 0\\}.$ (8)

The projection operator to $\Omega$ can be easily computed as

$\Pi_{\Omega}(V):=\arg\min_{U\in\Omega}\|U-V\|_{F}={\sqrt{K}\cdot(V)_{+}\over\|(V)_{+}\|_{F}},$
(9)

where $(V)_{+}=\max\\{V,0\\}$ denotes elementwise positive part of $V$. Then
the BM can be converted to the an equality-constrained problem over $\Omega$:

$\min_{U\in\Omega}\Big{\\{}\langle
A,UU^{T}\rangle:UU^{T}\mathbf{1}_{n}=\mathbf{1}_{n}\Big{\\}}.$ (10)

Using the standard augmented Lagrangian method (Nocedal and Wright, 2006,
Chapter 17), we see that the above (6) is also equivalent to

$\min_{U\in\Omega}\Big{\\{}\langle
L\cdot\text{Id}_{n}+A,UU^{T}\rangle+{\beta\over
2}\|UU^{T}\mathbf{1}_{n}-\mathbf{1}_{n}\|_{2}^{2}:UU^{T}\mathbf{1}_{n}=\mathbf{1}_{n}\Big{\\}},$
(11)

where $\beta>0$ is a penalty parameter and $L\in(0,\lambda)$ for some proper
choice of $\lambda$ motivated from the optimal dual variable for the trace
constraint $\|U\|_{F}^{2}=\operatorname{tr}(Z)=K$. The augmented Lagrangian
for problem (11) is

$\mathcal{L}_{\beta}(U,y):=\langle L\cdot\text{Id}_{n}+A,UU^{T}\rangle+\langle
y,UU^{T}\mathbf{1}_{n}-\mathbf{1}_{n}\rangle+{\beta\over
2}\|UU^{T}\mathbf{1}_{n}-\mathbf{1}_{n}\|_{2}^{2}.$ (12)

Therefore, we consider the augmented Lagrangian iterations

$U_{\mathrm{new}}=\arg\min_{U\in\Omega}\mathcal{L}_{\beta}(U,y),\qquad
y_{\mathrm{new}}=y+\beta(U_{\mathrm{new}}U_{\mathrm{new}}^{T}\mathbf{1}_{n}-\mathbf{1}_{n}).$
(13)

Now, we consider minimizing $\mathcal{L}_{\beta}(U,y)$ over $U$, with a fixed
$\beta$ and $y$. Here, we observe that $\mathcal{L}_{\beta}(U,y)$ is simply a
quadratic polynomial over $U$, and therefore its gradient can be easily
derived as
$\nabla_{U}\mathcal{L}_{\beta}(U,y)=(2A+2L\cdot\text{Id}+\mathbf{1}_{n}\overline{y}^{T}+\overline{y}\mathbf{1}_{n}^{T})U$
where $\overline{y}=y+\beta(UU^{T}\mathbf{1}_{n}-\mathbf{1}_{n})$. Combining
this with the insight that it is easy to project onto $\Omega$, we can perform
the following projected gradient descent iterations

$\displaystyle U_{\mathrm{new}}$
$\displaystyle=\Pi_{\Omega}(U-\alpha\nabla_{U}\mathcal{L}_{\beta}(U,y))$

until $\mathcal{L}_{\beta}(U,y)$ is minimized. The complete algorithm is
summarized in Algorithm 1, whose per iteration time and space complexity are
both $O(nr)$. Since the NMF formulation (7) is a simpler version of our BM
formulation (6), a projected gradient descent algorithm for the solving the
former corresponds to the primal decent step of Algorithm 1 with $y=0$,
$\beta=0$ and a simpler projection operator $\Pi_{\mathbb{R}^{n\times
r}_{+}}(V)=(V)_{+}$.

##  4 Theoretical Analysis

In this section, we establish the local linear convergence rate of the BM
algorithm (Algorithm 1) for solving the $K$-means clustering. We shall work
with the standard Gaussian mixture model (GMM) that assumes the data
$X_{1},\dots,X_{n}$ are generated from the following mechanism: if $i\in
G_{k}^{*}$, then

$X_{i}=\mu_{k}+\varepsilon_{i},$ (14)

where $G_{1}^{*},\dots,G_{K}^{*}$ is a true (unknown) partition of $[n]$ we
wish to recover, $\mu_{1},\dots,\mu_{K}\in\mathbb{R}^{p}$ are the cluster
centers and $\varepsilon_{i}\sim N(0,I_{p})$ are i.i.d. standard Gaussian
noises. Let $n_{k}=|G_{k}^{\ast}|$ denote the size of $G_{k}^{\ast}$ and by
convention $n_{0}=0$. Since in the exactly recovery regime (Assumption A
below), there is no relaxation gap — any global optimum $U^{\ast}$ of the BM
problem (4) at $r=K$ corresponds to the unique global optimum $Z^{\ast}$ of
the SDP problem (6) through the relation $Z^{\ast}=U^{\ast}{U^{\ast}}^{T}$, it
is sufficient to focus our analysis on the $r=K$ case. We note that the $r=K$
case corresponds to an exact parameterization of the matrix rank, which is
standard in the analysis of algorithms based on the BM factorization (Ge et
al., 2017; Chi et al., 2019).

Algorithm 1 is formulated as an exact augmented Lagrangian method, in which
the primal subproblem $\min_{U\in\Omega}\mathcal{L}_{\beta}(U,y)$ is solved to
sufficiently high accuracy as to be viewed as an exact solution. To be
precise, Step 1 in the algorithm is simply projected gradient descent
iterations performed on the primal problem, and the dual step in Step 2 is
only taken once projected gradient descent has achieved sufficiently high
accuracy. We will soon prove that projected gradient descent converges at a
rapid linear rate, and therefore it is very reasonable to run Step 1 until it
reaches the numerical floor.

Under exact minimization, it is well-known that the dual multipliers converge
at a linear rate under second-order sufficiency conditions, and that this in
turn implies linear convergence of the primal minimizers towards the true
solution. The typical proof establishes that, once the penalty $\beta$ is set
sufficiently large, the primal subproblem has a unique global minimizer
$U(y,\beta)$, whose quality is comparable to the quality of the multipliers
$y$ used to compute it. In turn, the dual multipliers $y$ converge linearly,
as the dual updates are simply gradient ascent steps on an underlying dual
functional. The following results are very classical; their wording is taken
from Proposition 1 and 2 of Bertsekas (1976), though a more modern proof
without constraint qualification assumptions is given in Fernández and Solodov
(2012). (See also Bertsekas (2014, Proposition 2.7))

######  Proposition 1 (Existence and quality of primal minimizer).

Let $U^{*}$ denote a local minimizing point that satisfies the second-order
sufficient conditions for an isolated local minimum with respect to
multipliers $y^{*}$. Then, there exists a scalar $\beta^{*}\geqslant 0$ such
that, for every $\beta>\beta^{*}$ the augmented Lagrangian
$\mathcal{L}_{\beta}(U,y)$ has a unique minimizing point $U(y,\beta)$ within
an open ball centered at $U^{*}$. Furthermore, there exists some scalar $M>0$
such that $\|U(y,\beta)-U^{*}\|_{F}\leqslant(M/\beta)\|y-y^{*}\|.$

Proposition 1 is directly applicable to problem (6) when $r=K$, because the
global minimum $U^{*}$ is made isolated by the nonnegativity constraint
$U\geqslant 0$; for further details about solution uniqueness, refer to the
discussion surrounding (15) below. In fact, it is possible to generalize
Proposition 1 to all $r\geqslant K$, even when the global minimum $U^{*}$ is
no longer isolated, although this would require further specializing
Proposition 1 to our specific problem. The key insight is that, in the $r>K$
case, the _closest_ global minimum $U^{*}$ to the current iterate $U$
continues to satisfy all the properties satisfied by the isolated global
minimum when $r=K$. This is essentially the same idea as (Ge et al., 2017,
Definition 6) and (Chi et al., 2019, Lemma 4). Nevertheless, we leave the
details of this extension as future work and focus on proving linear
convergence in the $r=K$ case in this paper.

######  Proposition 2 (Linear convergence of dual multipliers).

Under the same assumptions as Proposition 1, define the sequence

$y_{k+1}=y_{k}+\beta(U_{k}U_{k}^{T}\mathbf{1}_{n}-\mathbf{1}_{n})\quad\text{
where }U_{k}=U(y_{k},\beta)\equiv\arg\min_{U}\mathcal{L}_{\beta}(U,y_{k}).$

Then, there exists a radius $R>0$ such that, if $\|y_{0}-y^{*}\|\leqslant R$,
then $y_{k}\to y^{*}$ converges linearly

$\limsup_{k\to\infty}\frac{\|y_{k+1}-y^{*}\|}{\|y_{k}-y^{*}\|}\leqslant\frac{M}{\beta}.$

Therefore, almost all of the difficulty in our analysis amounts to showing
that projected gradient descent is able to solve the primal subproblem
$\min_{U\in\Omega}\mathcal{L}_{\beta}(U,y)$ efficiently. Below is our main
theorem for establishing the local linear rate of convergence of the proposed
BM solution when $r\geqslant K$ and the dual $y=y^{\ast}$ for the $K$-means
clustering under the GMM. Recall that in the exact recovery regime, the BM
formulation (6) of the $K$-means problem admits a global minimum $U^{\ast}$,
where $U^{\ast}_{is}=1/\sqrt{n_{s}}$ if $i\in G_{s}^{\ast}$, $s\in[K]$ and
$U^{\ast}_{is}=0$ otherwise, so that $Z^{\ast}=U^{\ast}{U^{\ast}}^{T}$.

When $r>K$, optimal solution $U$ to the BM formulation (6) is not unique since
it is possible to construct $U=U^{\ast}Q$ for some orthogonal transformation
$Q\in\mathbb{R}^{r\times r}$ corresponds to the same cluster membership matrix
$Z=UU^{T}=U^{\ast}{U^{\ast}}^{T}$ while maintaining the nonnegative constraint
$U\geqslant 0$. In fact, any feasible solution $U\in\mathbb{R}^{n\times r}$
such that $UU^{T}=Z^{\ast}$ must take the following block diagonal
form111Without loss of generality, we assume that data indices are reordered
so that $G_{k}^{\ast}=\\{n_{k-1}+1,\ldots,n_{k}\\}$. (up to column
permutations, see Theorem 4 in Appendix A), due to the nonnegative constraint
$U\geqslant 0$:

$\displaystyle
U=\begin{bmatrix}a_{1,1}\mathbf{1}_{n_{1}},\dots,a_{1,m_{1}}\mathbf{1}_{n_{1}}&0&\cdots&0\\\
0&a_{2,1}\mathbf{1}_{n_{2}},\dots,a_{2,m_{2}}\mathbf{1}_{n_{2}}&\cdots&0\\\
\vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&a_{K,1}\mathbf{1}_{n_{K}},\dots,a_{K,m_{K}}\mathbf{1}_{n_{K}}\end{bmatrix},$
(15)

for some block partition sizes $\mathbf{m}=(m_{1},m_{2},\ldots,m_{K})$
satisfying $m_{k}\geqslant 1$, $\sum_{k}m_{k}=r$ and within block weights
$\mathbf{a}_{k}=(a_{k,1},a_{k,2},\ldots,a_{k,m_{k}})$ satisfying
$a_{k,\ell}\geqslant 0$, $\sum_{\ell}a_{k,\ell}^{2}=\frac{1}{n_{k}}$ for all
$k\in[K]$. We will call any matrix $U$ whose sparsity/nonzero pattern
coincides with (15) (up to column permutations) as an $\mathbf{m}$-block
diagonal matrix222More precisely, an $\mathbf{m}$-block diagonal matrix is any
$n\times r$ matrix whose support (locations of non-zero entries) is a subset
of that of (15).. Let $\mathcal{G}_{\mathbf{m}}$ denote the set of all
$\mathbf{m}$-block diagonal matrices with nonnegative entries. As a direct
consequence of this characterization of optimal solutions, the non-uniqueness
issue disappears when $r=K$ (unique solution with $m_{1}=\cdots=m_{K}=1$);
then the primal subproblem has an isolated local minimum (which is also
global) with respect to $y^{\ast}$ and Propositions 1, 2 can be applied.

Before formally introducing the theoretical results, let us introduce some
more notations. We denote the optimal solution (15) as $U^{\mathbf{a},\ast}$
to emphasize its dependence on $\mathbf{a}$ (which also determines
$\mathbf{m}$). For any matrix $A\in\mathbb{R}^{n\times r}$, we use
$A_{S(\mathbf{a})}$ to denote the matrix of the same size as the restriction
of $A$ onto the support $S(\mathbf{a})$ of $U^{\mathbf{a},\ast}$, that is,
$[A_{S(\mathbf{a})}]_{ij}=A_{ij}$ if $[U^{\mathbf{a},\ast}]_{ij}\neq 0$ and
$[A_{S(\mathbf{a})}]_{ij}=0$ otherwise; and let
$A_{S(\mathbf{a})^{c}}=A-A_{S(\mathbf{a})}$. For any $\mathbf{m}$, we denote
the set of all orthogonal transformations preserving the $\mathbf{m}$-block
diagonal structure as $\mathcal{F}_{\mathbf{m}}$. We denote the initialization
of the algorithm as $U^{0}$, and define $U^{t}$ recursively by projected
gradient descent updating formula

$\displaystyle U^{t}$
$\displaystyle=\Pi_{\Omega}\big{(}U^{t-1}-\alpha\nabla_{U}\mathcal{L}_{\beta}(U^{t-1},y^{\ast})\big{)},\quad
t\in\mathbb{Z}_{+}.$

Assumption A (Exact recovery regime): For the equal cluster size case, the
minimal separation $\Theta_{\min}:=\min_{1\leqslant j\neq k\leqslant
K}\|\mu_{j}-\mu_{k}\|_{2}$ satisfies
$\Theta_{\min}\geqslant(1+\tilde{\alpha})\overline{\Theta}$ for some
$\tilde{\alpha}>0$, where $\overline{\Theta}$ is the information-theoretic
optimal threshold given in (5). For unbalanced cluster case, we assume the
minimal separation satisfies the conditions in Proposition 5 in Appendix B.2.

######  Theorem 3 (Convergence of projected gradient descent).

Suppose Assumption A holds. If there exists some block partition sizes
$\mathbf{m}$ and an associated optimal solution $U^{\mathbf{a},\ast}$ with
some within block weights
$\mathbf{a}_{k}=(a_{k,1},a_{k,2},\ldots,a_{k,m_{k}})$ satisfying
$\bar{a}\leqslant c\underline{a}$333Here,
$\underline{a}=\min_{k\in[K]}\min_{\ell\in[m_{k}]}\\{a_{k,\ell}\\}$ and
$\bar{a}=\max_{k\in[K]}\max_{\ell\in[m_{k}]}\\{a_{k,\ell}\\}$., $c>0$, such
that the initialization discrepancy $\Delta^{0}=U^{0}-U^{\mathbf{a},\ast}$
satisfies

$\displaystyle\|\Delta^{0}_{S(\mathbf{a})^{c}}\|_{\infty}\leqslant\frac{C_{1}\underline{a}\underline{n}}{p/\Theta^{2}_{\min}+m}\quad\mbox{and}\quad\|\Delta^{0}\|_{F}\leqslant\frac{C_{3}}{[1+C_{2}\alpha(L+n\beta+\Theta^{2}_{\max})]^{I}}$
$\displaystyle\cdot\,\min\left\\{\frac{\underline{n}}{n},\frac{{\underline{a}\underline{n}\Theta_{\min}^{2}}}{rn^{2}\bar{a}^{3}(L+n\beta+n\Theta^{2}_{\max})},\frac{(1-\gamma)\underline{a}\underline{n}\Theta_{\min}^{2}/K}{\Theta_{\max}\sqrt{K\log
n}+\sqrt{n}(p+K\log n)},(1-\gamma)\underline{a}\sqrt{\underline{n}}\right\\},$

where $m=\min_{k\neq l}\frac{2n_{k}n_{l}}{n_{k}+n_{l}}$,
$\underline{n}=\min\\{n_{k}\\}$, $\Theta_{\max}=\max_{1\leqslant j\neq
k\leqslant K}\|\mu_{j}-\mu_{k}\|_{2}$, and
$I=(4\alpha)^{-1}(p+{C_{1}}m\Theta^{2}_{\min})^{-1}$. If we take
$\beta=\frac{C_{1}\underline{a}\underline{n}\Theta_{\min}^{2}}{16rn^{2}\bar{a}^{3}},\;L=p+\frac{C_{1}\Theta^{2}_{\min}}{4}(m-\frac{\underline{a}\underline{n}^{2}}{8rn^{2}\bar{a}^{3}})$
in (12) and choose step size
$\alpha=\frac{C_{4}\underline{a}\underline{n}^{2}\Theta_{\min}^{2}}{rn^{2}\bar{a}^{3}(L+n\beta)^{2}}$,
then it holds with probability at least $1-n^{-C_{5}}$ that for any
$t\geqslant I$,

$\displaystyle
U^{t}\in\mathcal{G}_{\mathbf{m}}\quad\mbox{and}\quad\inf_{Q\in\mathcal{F}_{\mathbf{m}}}\|U^{t+1}-U^{\ast}Q\|_{F}\leqslant\gamma\,\inf_{Q\in\mathcal{F}_{\mathbf{m}}}\|U^{t}-U^{\ast}Q\|_{F},$

where
$\gamma^{2}=1-\frac{C_{6}(\underline{a}\underline{n}^{2}\Theta_{\min}^{2})^{2}}{(rn^{2}\bar{a}^{3})^{2}(L+n\beta)^{2}}$
and $C_{1},\dots,C_{6}$ are universal constants.

Proof sketches. Our proof strategy is to divide the convergence of the
projected gradient descent for solving the primal subproblem into two phases.
In phase one, we show that as long as the initialization $U^{0}$ falls into
the basin of attraction of some optimal point $U^{\mathbf{a},\ast}$ taking the
form of (15) with block diagonal sizes $\mathbf{m}$, then after at most $I$
iterations, the iterate will become block diagonal, that is, falls into set
$\mathcal{G}_{\mathbf{m}}$. This is because for all $U$ near the global
minimum $U^{\mathbf{a},\ast}$, the entries of the gradient
$\nabla_{U}\mathcal{L}(U,y^{\ast})$ located outside the support
$S(\mathbf{a})$ of $U^{\mathbf{a},\ast}$ are strictly lower bounded away from
zero. As a consequence, the projected gradient descent algorithm in phase one
keeps making progress in identifying the correct sparsity pattern from
$U^{\mathbf{a},\ast}$ — it will quickly drag the positive entries of $U$
located outside $S(\mathbf{a})$ towards zero (see Lemma 12 in Appendix B.4).
We then show that once the iterate becomes $\mathbf{m}$-block diagonal, the
algorithm enters phase two. In phase two, the iterate remains
$\mathbf{m}$-block diagonal. Moreover, the projected gradient descent attains
a linear convergence rate, since the objective function
$\mathcal{L}(\,\cdot\,,\,y^{\ast})$ is restricted strongly convex at
$U^{\ast}$ within $\mathcal{G}_{\mathbf{m}}$. More precisely, there exists
some $\tilde{\beta}>0$, such that for any $U\in\mathcal{G}_{\mathbf{m}}$, we
have (see Lemma 17 in Appendix B.6)

$\displaystyle\langle\nabla^{2}_{U}\mathcal{L}(U^{\ast},y^{\ast})[\Delta],\Delta\rangle\geqslant\tilde{\beta}\,\|\Delta\|_{F}^{2},\
\mbox{for}\ \Delta=UQ_{U}^{T}-U^{\ast},$

where
$Q_{U}=\operatorname{argmin}_{Q\in\mathcal{F}_{\mathbf{m}}}\|U-U^{\ast}Q\|_{F}$.
A complete proof of this theorem can be found in Appendix B.

Implications of the theorem: According to the theorem, under the equal cluster
size setting with moderate $p=O(n\log n)$ and $\Theta_{\min}^{2}=\Omega(\log
n)$ (corresponding to the exact recovery regime), if we optimally choose
$L=O(n\Theta_{\min}^{2}/K)$, $\beta=O(\Theta_{\min}^{2}/K^{3})$ and
$\alpha^{-1}=O(K^{2}n\Theta_{\min}^{2})$, then phase one will last at most
$I=O(K^{3})$ iterations; and after entering phase two, the contraction rate of
the projected gradient descent is $\gamma=1-O(K^{-6})$. These choices make the
iteration complexity of the projected GD for solving the primal subproblem of
order $O(K^{6})$ (modulo logarithmic factors), and the overall time complexity
of the primal-dual algorithm becomes of order $O(K^{6}nr)$. In addition, the
conditions on the initialization $U^{0}$ simplify to
$\|\Delta^{0}_{S(\mathbf{a})^{c}}\|_{\infty}=O(K/\sqrt{nr})$ and
$\|\Delta^{0}\|_{F}=O\big{(}r^{-0.5}K^{-5.5}\min\\{1,\,K^{-2.5}\Theta_{\min}^{2}/\log
n)\\}\big{)}$. Given that the optimal solution $U^{\ast}$ satisfies
$\|U^{\ast}\|_{\infty}=\sqrt{K/n}$ and $\|U^{\ast}\|_{F}=\sqrt{K}$, the
initialization condition only demands a constant relative error with respect
to $n$, making it a reasonable requirement. For example, under mild
conditions, rapid $K$-means algorithms like Lloyd’s algorithm (Lloyd, 1982; Lu
and Zhou, 2016) can be employed to construct an initialization that meets this
condition with high probability.

##  5 Numerical Experiments

Figure 2: Log-scale plots with errorbars of mis-clustering error (left most)
and time cost (in the middle) as sample size $n$ increases and the convergence
of BM over iterations (right most). The first two plots are partial due to
huge space (time) complexity for SC (SDP) when sample size is large.

Figure 3: Boxplots of mis-clustering error (with means) for CyTOF dataset (on
the left) and CIFAR-10 (on the right) among five different methods. The plots
are partial due to huge space (time) complexity for SC (SDP) when sample size
is large.

We present numerical results to assess the effectiveness of the proposed BM
method. We first conduct two simulation experiments using GMM to evaluate the
convergence and compare the performance of BM with other methods. Then, we
perform a comparison using two real datasets. One of the competing methods in
our comparison is clustering based on solving the NMF formulation (7).
Specifically, we employ the projected gradient descent algorithm, which is a
simplified version of Algorithm 1 discussed in Section 3, to implement this
method. We adopt random initialization and set the same $r$ for both NMF and
BM to ensure a fair comparison.

Performance for BM for GMM. Our goal is to compare the statistical performance
and time complexity of BM, NMF, the SDPNAL+ solver (Yang et al., 2015) for
solving SDP, spectral clustering (SC), and $K$-means++ (KM) under Gaussian
mixture models (GMM). In our setting, we choose cluster numbers $K=4$ and
place the centers of Gaussian distributions at the vertices of a simplex such
that $\Theta_{\min}^{2}=\gamma\,\overline{\Theta}^{2}$ with $\gamma=0.64$,
where $\overline{\Theta}$ is the sharp threshold defined in Equation (5). The
sample size ranges from $n=400$ to $n=57,600$, the dimension is $p=20$, and we
set the rank of the BM relaxation to be $r=2K$. The results are summarized in
the first two plots of Figure 2, with each case repeated $50$ times. From the
first plot, we observe that the mis-clustering error of SDP and BM coincides;
the error remains stable and decreases as the sample size $n$ increases.
However, NMF, KM, and SC exhibit large variance and are far from achieving
exact recovery. The second plot indicates that SDP and SC have super-linear
time complexity, while the log-scale curves of our BM approach, KM, and NMF
are nearly parallel, indicating that they all achieve linear time complexity.
In particular, SDP becomes time-consuming when $n$ is larger than $2000$; when
$n=40,000$, the space complexity issue makes SC infeasible.

Linear convergence of BM. To analyze the convergence of Algorithm 1, we
maintain the same setting as described earlier. However, we now fix $n=1000$
and consider minimum distance between centers
$\Theta_{\min}^{2}=\gamma\,\overline{\Theta}^{2}$ with $\gamma=1.44$. We
explore three cases: $r=K$, $r=2K$, and $r=20K$, and set the number of
iterations for Step 1 (primal update) to be $100$. The convergence of our BM
approach is shown in the third plot of Figure 2, with each result based on
$30$ replicates. We measure the relative difference of $U$ compared to optimal
solution $U^{\ast}$ using the Frobenius norm
$\|UU^{T}-U^{\ast}(U^{\ast})^{T}\|_{F}/\|U^{\ast}(U^{\ast})^{T}\|_{F}$. From
the third plot, we observe that our BM approach achieves linear convergence in
this setting. The curves for $r=K$, $r=2K$, and $r=20K$ overlap closely,
indicating that the choice of $r$ does not significantly affect the
convergence rate.

CyTOF dataset. This mass cytometry (CyTOF) dataset consists of protein
expression levels for $N=265,627$ cells with $p=32$ protein markers. The
dataset contains $14$ clusters or gated cell populations. In our experiment,
we select the labeled data for individual $H_{1}$. From these labeled data, we
uniformly sample $n$ data points from $K=4$ unbalanced clusters, specifically
clusters 2, 7, 8, and 9, which together contain a total of $46,258$ samples.
The results are presented in Figure 3, where we conduct all methods for $100$
replicates for each value of $n$. We consider three cases: case 1 corresponds
to $n=1800$, case 2 corresponds to $n=20,000$, and case 3 corresponds to
$n=46,258$. From the plot, we observe that BM and SDP remain stable, while KM,
SC exhibit significant variance, and NMF produces many outliers. In case 1,
the mis-clustering error of BM is comparable to that of SDP, indicating that
the performance of BM can be as good as SDP. However, when $n$ is on the order
of $1000$, the time complexity of SDP becomes dramatically high, which is at
least $O(n^{3.5})$. Additionally, when $n=46,258$, the space complexity of at
least $O(n^{2})$ becomes an issue for SC.

CIFAR-10 dataset. We conduct a comparison of all methods on the CIFAR-10
dataset, which comprises 60,000 colored images of size $32\times 32\times 3$,
divided into 10 clusters. Firstly, we apply the Inception v3 model with
default settings and perform PCA to reduce the dimensionality to $p=50$. The
test set consists of 10,000 samples with 10 clusters. In our experiment, we
focus on four clusters ("bird", "deer", "ship", and "truck") with $K=4$. We
consider two cases: Case 1 involves random sub-sampling with $n=1800$ to
compare BM with SDP, and Case 2 uses $n=4000$. The results are displayed in
the second plot of Figure 3, based on 100 replicates. Similar to the previous
experiment, we observe that BM and SDP exhibit similar behavior and achieve
superior and more consistent performance compared to KM, SC, and NMF. SC and
NMF display significant variance, while KM produces many outliers.

## References

  * Aloise et al. [2009] Daniel Aloise, Amit Deshpande, Pierre Hansen, and Preyas Popat.  Np-hardness of euclidean sum-of-squares clustering.  _Machine learning_ , 75(2):245–248, 2009.
  * Arthur and Vassilvitskii [2007] David Arthur and Sergei Vassilvitskii.  K-means++ the advantages of careful seeding.  In _Proceedings of the eighteenth annual ACM-SIAM symposium on Discrete algorithms_ , pages 1027–1035, 2007.
  * Bertsekas [1976] Dimitri P Bertsekas.  Multiplier methods: A survey.  _Automatica_ , 12(2):133–145, 1976.
  * Bertsekas [2014] Dimitri P Bertsekas.  _Constrained optimization and Lagrange multiplier methods_.  Academic press, 2014.
  * Burer and Monteiro [2003] Samuel Burer and Renato D. C. Monteiro.  A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization.  _Mathematical Programming_ , 95(2):329–357, 2003\.  doi: 10.1007/s10107-002-0352-8.  URL https://doi.org/10.1007/s10107-002-0352-8.
  * Chen and Yang [2021] Xiaohui Chen and Yun Yang.  Cutoff for exact recovery of gaussian mixture models.  _IEEE Transactions on Information Theory_ , 67(6):4223–4238, 2021.
  * Chi et al. [2019] Yuejie Chi, Yue M Lu, and Yuxin Chen.  Nonconvex optimization meets low-rank matrix factorization: An overview.  _IEEE Transactions on Signal Processing_ , 67(20):5239–5269, 2019.
  * Cichocki and Phan [2009] Andrzej Cichocki and Anh-Huy Phan.  Fast Local Algorithms for Large Scale Nonnegative Matrix and Tensor Factorizations.  _IEICE Transactions on Fundamentals of Electronics Communications and Computer Sciences_ , 92(3):708–721, January 2009.  doi: 10.1587/transfun.E92.A.708.
  * Dasgupta [2007] Sanjoy Dasgupta.  The hardness of $k$-means clustering.  _Technical Report CS2007-0890, University of California, San Diego_ , 2007.
  * Ding et al. [2005] Chris Ding, Xiaofeng He, and Horst D Simon.  On the equivalence of nonnegative matrix factorization and spectral clustering.  In _Proceedings of the 2005 SIAM international conference on data mining_ , pages 606–610. SIAM, 2005.
  * Dua and Graff [2017] Dheeru Dua and Casey Graff.  UCI machine learning repository, 2017.  URL http://archive.ics.uci.edu/ml.
  * Fei and Chen [2018] Yingjie Fei and Yudong Chen.  Hidden integrality of sdp relaxations for sub-gaussian mixture models.  In Sébastien Bubeck, Vianney Perchet, and Philippe Rigollet, editors, _Proceedings of the 31st Conference On Learning Theory_ , volume 75 of _Proceedings of Machine Learning Research_ , pages 1931–1965. PMLR, 06–09 Jul 2018.  URL https://proceedings.mlr.press/v75/fei18a.html.
  * Fernández and Solodov [2012] Damián Fernández and Mikhail V Solodov.  Local convergence of exact and inexact augmented lagrangian methods under the second-order sufficient optimality condition.  _SIAM Journal on Optimization_ , 22(2):384–407, 2012.
  * Ge et al. [2017] Rong Ge, Chi Jin, and Yi Zheng.  No spurious local minima in nonconvex low rank problems: A unified geometric analysis.  In _International Conference on Machine Learning_ , pages 1233–1242. PMLR, 2017.
  * Giraud and Verzelen [2018] Christophe Giraud and Nicolas Verzelen.  Partial recovery bounds for clustering with the relaxed $k$means.  _Math. Stat. Learn._ , (3/4):317–374, 2018.
  * Han et al. [2023] Xiao Han, Xin Tong, and Yingying Fan.  Eigen selection in spectral clustering: A theory-guided practice.  _Journal of the American Statistical Association_ , 118(541):109–121, 2023.  doi: 10.1080/01621459.2021.1917418.  URL https://doi.org/10.1080/01621459.2021.1917418.
  * He et al. [2011] Zhaoshui He, Shengli Xie, Rafal Zdunek, Guoxu Zhou, and Andrzej Cichocki.  Symmetric nonnegative matrix factorization: Algorithms and applications to probabilistic clustering.  _IEEE Transactions on Neural Networks_ , 22(12):2117–2131, 2011.  doi: 10.1109/TNN.2011.2172457.
  * Kim et al. [2014] Jingu Kim, Yunlong He, and Haesun Park.  Algorithms for nonnegative matrix and tensor factorizations: a unified view based on block coordinate descent framework.  _Journal of Global Optimization_ , 58(2):285–319, 2014.  doi: 10.1007/s10898-013-0035-4.  URL https://doi.org/10.1007/s10898-013-0035-4.
  * Kuang et al. [2015] Da Kuang, Sangwoon Yun, and Haesun Park.  Symnmf: nonnegative low-rank approximation of a similarity matrix for graph clustering.  _Journal of Global Optimization_ , 62:545–574, 2015.
  * Kulis et al. [2007] Brian Kulis, Arun C. Surendran, and John C. Platt.  Fast low-rank semidefinite programming for embedding and clustering.  In Marina Meila and Xiaotong Shen, editors, _Proceedings of the Eleventh International Conference on Artificial Intelligence and Statistics_ , volume 2 of _Proceedings of Machine Learning Research_ , pages 235–242, San Juan, Puerto Rico, 21–24 Mar 2007. PMLR.  URL https://proceedings.mlr.press/v2/kulis07a.html.
  * Lloyd [1982] Stuart Lloyd.  Least squares quantization in pcm.  _IEEE Transactions on Information Theory_ , 28:129–137, 1982.
  * Löffler et al. [2021] Matthias Löffler, Anderson Y. Zhang, and Harrison H. Zhou.  Optimality of spectral clustering in the Gaussian mixture model.  _The Annals of Statistics_ , 49(5):2506 – 2530, 2021.  doi: 10.1214/20-AOS2044.  URL https://doi.org/10.1214/20-AOS2044.
  * Lu and Zhou [2016] Yu Lu and Harrison H Zhou.  Statistical and computational guarantees of lloyd’s algorithm and its variants.  _arXiv preprint arXiv:1612.02099_ , 2016.
  * MacQueen [1967] J.B. MacQueen.  Some methods for classification and analysis of multivariate observations.  _Proc. Fifth Berkeley Sympos. Math. Statist. and Probability_ , pages 281–297, 1967.
  * Mixon et al. [2017] Dustin G Mixon, Soledad Villar, and Rachel Ward.  Clustering subgaussian mixtures by semidefinite programming.  _Information and Inference: A Journal of the IMA_ , 6(4):389–415, 03 2017.  ISSN 2049-8764.  doi: 10.1093/imaiai/iax001.  URL https://doi.org/10.1093/imaiai/iax001.
  * Ng et al. [2001] Andrew Y. Ng, Michael I. Jordan, and Yair Weiss.  On spectral clustering: Analysis and an algorithm.  In _Advances in Neural Information Processing Systems_ , pages 849–856. MIT Press, 2001.
  * Nocedal and Wright [2006] Jorge Nocedal and Stephen J. Wright.  _Numerical optimization_.  Springer series in operations research and financial engineering. Springer, New York, NY, 2. ed. edition, 2006.  ISBN 978-0-387-30303-1.  URL http://gso.gbv.de/DB=2.1/CMD?ACT=SRCHA&SRT=YOP&IKT=1016&TRM=ppn+502988711&sourceid=fbw_bibsonomy.
  * Peng and Wei [2007] Jiming Peng and Yu Wei.  Approximating $k$-means-type clustering via semidefinite programming.  _SIAM J. OPTIM_ , 18(1):186–205, 2007.
  * Royer [2017] Martin Royer.  Adaptive clustering through semidefinite programming.  In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, _Advances in Neural Information Processing Systems 30_ , pages 1795–1803. Curran Associates, Inc., 2017.
  * von Luxburg [2007] Ulrike von Luxburg.  A tutorial on spectral clustering.  _Statistics and Computing_ , 17(4):395–416, 2007\.
  * Wang and Zhang [2012] Yu-Xiong Wang and Yu-Jin Zhang.  Nonnegative matrix factorization: A comprehensive review.  _IEEE Transactions on knowledge and data engineering_ , 25(6):1336–1353, 2012.
  * Xu et al. [2003] Wei Xu, Xin Liu, and Yihong Gong.  Document clustering based on non-negative matrix factorization.  In _Proceedings of the 26th Annual International ACM SIGIR Conference on Research and Development in Informaion Retrieval_ , SIGIR ’03, page 267–273, New York, NY, USA, 2003. Association for Computing Machinery.  ISBN 1581136463.  doi: 10.1145/860435.860485.  URL https://doi.org/10.1145/860435.860485.
  * Yang et al. [2015] Liuqin Yang, Defeng Sun, and Kim-Chuan Toh.  Sdpnal+: a majorized semismooth newton-cg augmented lagrangian method for semidefinite programming with nonnegative constraints.  _Mathematical Programming Computation_ , 7(3):331–366, 2015.  doi: 10.1007/s12532-015-0082-6.  URL https://doi.org/10.1007/s12532-015-0082-6.
  * Zhu et al. [2018] Zhihui Zhu, Xiao Li, Kai Liu, and Qiuwei Li.  Dropping symmetry for fast symmetric nonnegative matrix factorization.  _Advances in Neural Information Processing Systems_ , 31, 2018.
  * Zhuang et al. [2022] Yubo Zhuang, Xiaohui Chen, and Yun Yang.  Wasserstein $k$-means for clustering probability distributions.  In _Advances in Neural Information Processing Systems_ , 2022.
  * Zhuang et al. [2023] Yubo Zhuang, Xiaohui Chen, and Yun Yang.  Likelihood adjusted semidefinite programs for clustering heterogeneous data.  In _International Conference on Machine Learning (ICML)_ , 2023.

## Supplementary Material

This supplementary material contains proofs of all theoretical results in the
paper. We shall work with the standard Gaussian mixture model (GMM) that
assumes the data $X_{1},\dots,X_{n}$ are generated from the following
mechanism: if $i\in G_{k}$, then

$X_{i}=\mu_{k}+\varepsilon_{i},$ (16)

where $G_{1},\dots,G_{K}$ is a true (unknown) partition of $[n]$ we wish to
recover, $\mu_{1},\dots,\mu_{K}\in\mathbb{R}^{p}$ are the cluster centers and
$\varepsilon_{i}\sim N(0,I_{p})$ are i.i.d. standard Gaussian noises. Let
$n_{k}=|G_{k}|$ denote the size of $G_{k}$. Define
$Z^{\ast}:=U^{\ast}(U^{\ast})^{T},$ where

$\displaystyle
U^{\ast}=\begin{bmatrix}\frac{1}{\sqrt{n_{1}}}\mathbf{1}_{n_{1}}&0&\cdots&0&0&\cdots&0\\\
0&\frac{1}{\sqrt{n_{2}}}\mathbf{1}_{n_{2}}&\cdots&0&0&\cdots&0\\\
\vdots&\vdots&\ddots&\vdots&0&\vdots&\vdots\\\
0&0&\cdots&\frac{1}{\sqrt{n_{K}}}\mathbf{1}_{n_{K}}&0&\cdots&0\end{bmatrix}\in\mathbb{R}^{n\times
r}.$ (17)

##  Appendix A Characterization of Global Optima of BM Problem

######  Theorem 4 (feasible solutions).

For any $U\in\mathbb{R}^{n\times r}$ such that $UU^{T}=Z^{\ast}$ and
$U\geqslant 0,$ $U$ must take the following block form (up to column
permutations), due to the nonnegative constraint $U\geqslant 0$:

$\displaystyle
U=\begin{bmatrix}a_{1,1}\mathbf{1}_{n_{1}},\dots,a_{1,m_{1}}\mathbf{1}_{n_{1}}&0&\cdots&0\\\
0&a_{2,1}\mathbf{1}_{n_{2}},\dots,a_{2,m_{2}}\mathbf{1}_{n_{2}}&\cdots&0\\\
\vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&a_{K,1}\mathbf{1}_{n_{K}},\dots,a_{K,m_{K}}\mathbf{1}_{n_{K}}\end{bmatrix},$
(18)

for some block partition sizes $\mathbf{m}=(m_{1},m_{2},\ldots,m_{K})$
satisfying $m_{k}\geqslant 1$, $\sum_{k}m_{k}=r$ and within block weights
$\mathbf{a}_{k}=(a_{k,1},a_{k,2},\ldots,a_{k,m_{k}})$ satisfying
$a_{k,\ell}\geqslant 0$, $\sum_{\ell}a_{k,\ell}^{2}=\frac{1}{n_{k}}$ for all
$k\in[K]$.

_Proof._ Recall

$\displaystyle
U^{\ast}=\begin{bmatrix}\frac{1}{\sqrt{n_{1}}}\mathbf{1}_{n_{1}}&0&\cdots&0&0&\cdots&0\\\
0&\frac{1}{\sqrt{n_{2}}}\mathbf{1}_{n_{2}}&\cdots&0&0&\cdots&0\\\
\vdots&\vdots&\ddots&\vdots&0&\vdots&\vdots\\\
0&0&\cdots&\frac{1}{\sqrt{n_{K}}}\mathbf{1}_{n_{K}}&0&\cdots&0\end{bmatrix}.$

Then for any $Q=[q_{1},\dots,q_{r}]\in\mathcal{O}^{r\times r},$ we have

$U:=U^{\ast}Q^{T}=\left[\frac{1}{\sqrt{n_{1}}}q_{1}\mathbf{1}_{n_{1}}^{T},\dots,\frac{1}{\sqrt{n_{K}}}q_{K}\mathbf{1}_{n_{K}}^{T}\right]^{T}.$

$U$ is feasible hence $U^{\ast}Q^{T}\geqslant 0\implies q_{i}\geqslant
0,\forall i.$ Note that $\langle q_{i},q_{i}\rangle=1,\;\langle
q_{i},q_{j}\rangle=0,\forall i\neq j,$ then $U$ must have the form (18) (up to
column permutations). $\blacksquare$

##  Appendix B Proof of Theorem 3

First we will briefly mention about the BM formulation. Then we provide some
properties of dual construction for SDP [Chen and Yang, 2021] and introduce
some notations. Finally we provide the main theorems for phase 1 $\&$ 2 and
the necessary lemmas.

###  B.1 BM method

Assume $y$ is fixed. Denote

$f(U):=\mathcal{L}(U,y)=\langle L\cdot\text{Id}_{n}+A,UU^{T}\rangle+\langle
y,UU^{T}\mathbf{1}_{n}-\mathbf{1}_{n}\rangle+\frac{\beta}{2}\|UU^{T}\mathbf{1}_{n}-\mathbf{1}_{n}\|_{F}^{2},$
(19)

where $L\in(0,\lambda),\beta>0,\;\Omega:=\\{U\in\mathbb{R}^{n\times
r}:\|U\|_{F}^{2}=K,U\geqslant 0\\}.$

Then we can get the gradient

$\nabla
f(U)=(2A+2L\cdot\text{Id}+y\mathbf{1}_{n}^{T}+\mathbf{1}_{n}y^{T})U+\beta[\mathbf{1}_{n}(UU^{T}\mathbf{1}_{n}-\mathbf{1}_{n})^{T}+(UU^{T}\mathbf{1}_{n}-\mathbf{1}_{n})\mathbf{1}_{n}^{T}]U.$

###  B.2  Properties for SDP

First we list our main assumptions below.

Assumption 1. If there exist constants $C_{1}\in(0,1),$ $C_{3}>0$ such that

$\log
n\geqslant\frac{(1-C_{1})^{2}}{C_{1}^{2}}\frac{D_{1}n}{m},\;C_{3}\leqslant\frac{C_{1}^{2}}{(1-C_{1})^{2}}\frac{D_{2}}{K},\;m\geqslant\frac{4(1+C_{3})^{2}}{C_{3}^{2}},$
$\Theta_{\min}^{2}\geqslant\overline{\Theta}_{a}^{2}:=\frac{4\sigma^{2}(1+2C_{3})}{(1-C_{1})^{2}}\left(1+\sqrt{1+\frac{(1-C_{1})^{2}}{1+C_{3}}}\frac{p}{m\log
n}+D_{3}R_{n}\right)\log n,$

with

$R_{n}=\frac{(1-C_{1})^{2}}{(1+C_{3})\log n}\left(\frac{\sqrt{p\log
n}}{\underline{n}}+\frac{\log n}{\underline{n}}\right),$

where $C_{1}\in(0,1),C_{3}>0,$,$\underline{n}:=\min_{k\in[K]}\\{n_{k}\\}$ and
$m=\min_{k\neq l}\frac{2n_{k}n_{l}}{n_{k}+n_{l}}$. $\Theta^{2}_{\min}$ is the
minimum squared separation of order $O(\log n).$ And we denote
$\bar{n}:=\max_{k\in[K]}\\{n_{k}\\}$,$\Theta^{2}_{\max}$ to be the maximum
squared separation. Here $D_{1},D_{2},D_{3}$ are universal constants.

Assumption 2. Choose $(\lambda,y,B)$ to be the dual of SDP [Chen and Yang,
2021]. In particular, $\lambda=p+\frac{C_{1}}{4}m\Theta^{2}_{\min}.$

######  Proposition 5.

If Assumption 1 holds, then SDP achieves exact recovery with high probability
at least $1-D_{4}K^{2}n^{-C_{3}},$ for some constant $D_{4}.$

######  Proposition 6.

If Assumption 1 and 2 hold, then from Lemma VI.1 [Chen and Yang, 2021] we have
$\forall j\in G_{l},\;k\neq l,(k,l)\in[K]^{2},$

$\displaystyle
D_{k,l}(j):=[(2A+1y^{T}+y1^{T})_{G_{l},G_{k}}\mathbf{1}_{n_{k}}]_{j}\in[C_{1}\underline{n}\Theta^{2}_{\min},\bar{n}\Theta^{2}_{\max}],$

with probability at least $1-12n^{-C_{3}}$.

######  Proposition 7.

For any $v\in\mathbb{R}^{n},$ define $S(v):=v^{T}Av$,
$\Gamma_{K}:=\text{span}\\{\mathbf{1}_{G_{k}}:k\in[K]\\}^{\perp},$ which is
the orthogonal complement of the linear subspace of $\mathbb{R}^{n}$ spanned
by the vectors $\mathbf{1}_{G_{1}},\dots,\mathbf{1}_{G_{K}}$. Then Lemma VI.1
[Chen and Yang, 2021] implies

$\mathbb{P}\left(\max_{v\in\Gamma_{K},\;\|v\|=1}|S(v)|\geqslant(\sqrt{n}+\sqrt{p}+\sqrt{2t})^{2}\right)\leqslant
e^{-t},\;\forall t>0.$

######  Proposition 8.

Under Assumption 2, define

$W_{n}:=\lambda\text{Id}+\frac{1}{2}(y\mathbf{1}_{n}^{T}+\mathbf{1}_{n}y^{T})-B+A$

for some $B\in\\{\mathbb{R}^{n\times n}:\;B\geqslant
0,\;B_{G_{k},G_{k}}=0,\;\forall k\in[K]\\}.$ Then from Lemma IV.2 [Chen and
Yang, 2021] we have

$\|W_{n}\|_{op}\leqslant\lambda+C_{2}(n+\sqrt{p\log n}),\;W_{n}\succ 0,$

for some constant $C_{2}>0.$

###  B.3  Notations

For some block partition sizes $\mathbf{m}=(m_{1},m_{2},\ldots,m_{K})$
satisfying $m_{k}\geqslant 1$, $\sum_{k}m_{k}=r$, we define the form the set
$\mathcal{F}_{\mathbf{m}}$ that consists of the matrices $Q$ such that
$U=U^{\ast}Q$ has the form

$\displaystyle\begin{bmatrix}a_{1,1}\mathbf{1}_{n_{1}},\dots,a_{1,m_{1}}\mathbf{1}_{n_{1}}&0&\cdots&0\\\
0&a_{2,1}\mathbf{1}_{n_{2}},\dots,a_{2,m_{2}}\mathbf{1}_{n_{2}}&\cdots&0\\\
\vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&a_{K,1}\mathbf{1}_{n_{K}},\dots,a_{K,m_{K}}\mathbf{1}_{n_{K}}\end{bmatrix},$
(20)

where $\sum_{l}a_{k,l}^{2}=\frac{1}{n_{k}},\forall k\in[K],\;\sum_{k}m_{k}=r.$
And define the set $\mathcal{G}_{\mathbf{m}}$ that consists of the matrices
$U$ of the above form. Note that $Q$ exists for any $U$ of the form since
$UU^{T}=(U^{\ast})(U^{\ast})^{T}$.

Now we denote $\mathbf{a}_{k}=(a_{k,1},a_{k,2},\ldots,a_{k,m_{k}})$ and
$U^{\mathbf{a},\ast}$ as (20). Define
$\underline{a}:=\min_{k,l}\\{a_{k,l}\\}$, $\bar{a}:=\max_{k,l}\\{a_{k,l}\\}.$
For any matrix $A\in\mathbb{R}^{n\times r}$, we use $A_{S(\mathbf{a})}$ to
denote the matrix of the same size as the restriction of $A$ onto the support
$S(\mathbf{a})$ of $U^{\mathbf{a},\ast}$, that is,
$[A_{S(\mathbf{a})}]_{ij}=A_{ij}$ if $[U^{\mathbf{a},\ast}]_{ij}\neq 0$ and
$[A_{S(\mathbf{a})}]_{ij}=0$ otherwise; and let
$A_{S(\mathbf{a})^{c}}=A-A_{S(\mathbf{a})}$.

Let $\Pi_{\Omega}$ be the projection map on to the set $\Omega.$ i.e.,

$v=\Pi_{\Omega}(u)\Longleftrightarrow
v\in\Omega,\|v-u\|\leqslant\|\tilde{u}-u\|,\forall\tilde{u}\in\Omega.$

Recall that $\Omega:=\\{U\in\mathbb{R}^{n\times r}:\|U\|_{F}^{2}=K,U\geqslant
0\\}.$ We define

$V=\Pi_{\Omega}(W),\;W=U-\alpha\nabla f(U),\;\alpha>0.$

And denote $Q,Q^{t}\in\mathcal{F}_{\mathbf{m}}$ as

$Q=\text{argmin}_{\tilde{Q}\in\mathcal{F}_{\mathbf{m}}}\|U-U^{\ast}\tilde{Q}\|_{F},\;Q^{t}=\text{argmin}_{\tilde{Q}\in\mathcal{F}_{\mathbf{m}}}\|U^{t}-U^{\ast}\tilde{Q}\|_{F},\;\forall
t\geqslant 0.$

###  B.4  Main theorems for linear convergence near optimal point

Here we will state the main theorems formally with concrete forms for phases 1
and 2.

######  Theorem 9 (Convergence of Phase 1).

Define $U^{0}$ to be the initialization of BM method and define
$\Delta^{t}:=U^{t}-U^{\mathbf{a},\ast}$. Suppose $\bar{a}\leqslant
c\underline{a},$ for some constant $c.$ Then there exist
$(I,\eta)\in\mathbb{Z}_{+}\times\mathbb{R}_{+}$ s.t.
$U^{I}\in\mathcal{G}_{\mathbf{m}},\;\|U^{t}-U^{\mathbf{a},\ast}\|_{F}\leqslant\epsilon_{c},\forall
t\leqslant I,$ provided that

$\|\Delta^{0}_{S(\mathbf{a})^{c}}\|_{\infty}\leq\frac{C_{1}\underline{a}\underline{n}\Theta_{\min}^{2}}{32\lambda},\;\;\|\Delta^{0}\|_{F}\leqslant\frac{1}{\eta^{I}}\min\left\\{\epsilon_{c},\epsilon_{1},\frac{3}{4}\underline{a}\sqrt{\underline{n}}\right\\}.$

And for the choice of $(L,\beta)\in~{}(\ref{equ:BM}):$

$\frac{2(\lambda-L)}{\underline{n}}<\beta\leqslant\frac{C_{1}\underline{a}\underline{n}\Theta_{\min}^{2}}{12rn^{2}\bar{a}^{3}}.$

Here $\epsilon_{1}:=$

$\min\left\\{\frac{C_{1}\underline{a}\underline{n}\Theta_{\min}^{2}}{64K(\Theta_{\max}\sqrt{2K\log
n}+\sqrt{n}(p+2\sqrt{pK\log n}+2K\log
n))},\sqrt{\frac{C_{1}\underline{a}\underline{n}\Theta_{\min}^{2}}{96\beta
r\bar{a}\bar{n}}},\frac{C_{1}\underline{a}\underline{n}\Theta_{\min}^{2}}{96\beta
r\bar{a}^{2}\bar{n}\sqrt{n}}\right\\},$

$\epsilon_{c}:=\frac{(\lambda-L)\sqrt{K}}{2[(1+\sqrt{K})R_{1}+2(\lambda-L)\sqrt{K}+\sqrt{nr}\cdot\bar{a}\bar{n}\Theta^{2}_{\max}]},$

with probability at least $1-C_{4}n^{-C_{3}},$ where
$R_{1}=2C_{2}(n+\sqrt{p\log n})+2L+2/C_{1}\Theta^{2}_{\max}+12n\beta$. Recall
that $\lambda=p+\frac{C_{1}}{4}m\Theta^{2}_{\min},\;m=\min_{k\neq
l}\frac{2n_{k}n_{l}}{n_{k}+n_{l}}$ and $C_{1},C_{2},C_{3},C_{4}$ are some
constants.

_Proof._ This can be implied from Corollary 13 and Lemma 14. $\blacksquare$

######  Theorem 10 (Linear convergence of Phase 2).

Suppose $U^{0}\in\mathcal{G}_{\mathbf{m}},$ then with probability at least
$1-C_{4}n^{-C_{3}},$

$\|U^{t+1}-U^{\ast}Q^{t+1}\|_{F}\leqslant\gamma\|U^{t}-U^{\ast}Q^{t}\|_{F},\
\;\forall k,$

given
$\|U^{0}-U^{0}Q^{0}\|_{F}\leqslant\epsilon,\|U^{0}-U^{\mathbf{a},\ast}\|_{F}\leqslant\epsilon_{0},$
for some $\epsilon>0,\;\epsilon_{0}>0,\;\gamma\in(0,1)$; and $\alpha,L$ are
appropriately chosen.

_Proof._ This can be implied from Lemma 15 and Lemma 16. $\blacksquare$

######  Remark 11.

The exact values for $\gamma,\;\epsilon_{0},\;\epsilon$ are defined as
following: $\gamma^{2}=(1-\alpha\tilde{\beta}/2),$

$\epsilon=\min\left\\{\frac{\tilde{\beta}}{36\beta
n},\frac{(\lambda-L)\sqrt{K}}{2[\sqrt{K}R_{1}+2(\lambda-L)\sqrt{K}+n\Theta^{2}_{\max}]}\right\\},\;\epsilon_{0}:=\frac{1-\gamma}{2}\min\left\\{\epsilon_{1},\frac{3}{4}\underline{a}\sqrt{\underline{n}}\right\\},$

where step size should be upper bounded:
$\alpha<\frac{\tilde{\beta}}{2(2C_{2}(n+\sqrt{p\log
n})+2L+12n\beta)^{2}},\;\tilde{\beta}=\min\\{2[L-(\sqrt{n}+\sqrt{p}+\sqrt{2\log
n})^{2}],-2(\lambda-L)+\beta\underline{n}\\}.$ $\blacksquare$

If we combine phase 1 and 2 together, we will get the following theorem, which
is stronger version than Theorem 3 in the main paper.

Theorem 3 (Local linear rate of convergence) Recall
$\Delta^{t}:=U^{t}-U^{\mathbf{a},\ast}$. Suppose $\bar{a}\leqslant
c\underline{a},$ for some constant $c.$ Then there exists
$(I,\eta)\in\mathbb{Z}_{+}\times\mathbb{R}_{+}$ s.t. with probability at least
$1-C_{4}n^{-C_{3}},$

$\forall k\geqslant
I,\;U^{t}\in\mathcal{G}_{\mathbf{m}},\;\;\|U^{t+1}-U^{\ast}Q^{t+1}\|_{F}\leqslant\gamma\|U^{t}-U^{\ast}Q^{t}\|_{F},$

provided that

$\|\Delta^{0}_{S(\mathbf{a})^{c}}\|_{\infty}\leq\frac{C_{1}\underline{a}\underline{n}\Theta_{\min}^{2}}{32\lambda},\;\;\|\Delta^{0}\|_{F}\leqslant\frac{1}{\eta^{I}}\min\left\\{\epsilon,\frac{1-\gamma}{2}\epsilon_{1},\frac{1-\gamma}{2}\frac{3}{4}\underline{a}\sqrt{\underline{n}}\right\\}.$

And for $L,t$:

$\frac{2(\lambda-L)}{\underline{n}}<\beta\leqslant\frac{C_{1}\underline{a}\underline{n}\Theta_{\min}^{2}}{12rn^{2}\bar{a}^{3}},\;\;\;\alpha<\frac{\tilde{\beta}}{2(2C_{2}(n+\sqrt{p\log
n})+2L+12n\beta)^{2}}.$

Here $I=\frac{1}{4\alpha\lambda},\;\eta=1+\alpha R_{1}$.

###  B.5  Linear convergence near optimal point (Phase 1)

First we will provide the lemmas needed for Theorem 9, which will imply that
after at most $I$ iterations, the iterate will become block diagonal and fall
into set $\mathcal{G}_{\mathbf{m}}.$ Without loss of generality and for
notation simplicity, in the rest of the proofs we assume
$m_{1}=\cdots=m_{K}=s=r/K,$ $a_{k,1}=\cdots=a_{k,m_{k}}=a_{k},\forall k$, and
we denote $\mathcal{G}_{\mathbf{m}}$ as $\mathcal{G}$,
$\mathcal{F}_{\mathbf{m}}$ as $\mathcal{F}$.

######  Lemma 12.

Under certain conditions on initialization and $L,$ we have

$[U]_{i,\tau}-\alpha[\nabla
f(U)]_{i,\tau}\leqslant-\alpha\frac{C_{1}\underline{a}\underline{n}\Theta_{\min}^{2}}{8},\forall
i\in G_{l},s(k-1)<\tau\leqslant sk,\;l\neq k.$

_Proof._ The proof can then be direct to check for general unequal cases. By
definition we have

$\nabla
f(U)=(2A+2L\cdot\text{Id}+y\mathbf{1}_{n}^{T}+\mathbf{1}_{n}y^{T})U+\beta[\mathbf{1}_{n}(UU^{T}\mathbf{1}_{n}-\mathbf{1}_{n})^{T}+(UU^{T}\mathbf{1}_{n}-\mathbf{1}_{n})\mathbf{1}_{n}^{T}]U.$

WLOG, we suppose $\tau=1.$ Now suppose $\Delta=U-U^{\mathbf{a},\ast},$
$\Delta=[{d}_{1},\dots d_{r}].$ Then

$\displaystyle[\nabla f(U)]_{i,1}$
$\displaystyle=2L\cdot({d}_{1})_{i}+e_{i}^{T}(2A+y\mathbf{1}_{n}^{T}+\mathbf{1}_{n}y^{T})a_{1,1}\mathbf{1}_{G_{1}}+e_{i}^{T}(2A+y\mathbf{1}_{n}^{T}+\mathbf{1}_{n}y^{T}){d}_{1}$
$\displaystyle+\beta
e_{i}^{T}[\mathbf{1}_{n}(UU^{T}\mathbf{1}_{n}-\mathbf{1}_{n})^{T}+(UU^{T}\mathbf{1}_{n}-\mathbf{1}_{n})\mathbf{1}_{n}^{T}]a_{1,1}\mathbf{1}_{G_{1}}+r_{1}$
$\displaystyle=2L\cdot(d_{1})_{i}+a_{1,1}D_{1,l}(i)+e_{i}^{T}(2A+y\mathbf{1}_{n}^{T}+\mathbf{1}_{n}y^{T}){d}_{1}$
$\displaystyle+\beta
e_{i}^{T}[\mathbf{1}_{n}(UU^{T}\mathbf{1}_{n}-\mathbf{1}_{n})^{T}+(UU^{T}\mathbf{1}_{n}-\mathbf{1}_{n})\mathbf{1}_{n}^{T}]a_{1,1}\mathbf{1}_{G_{1}}+r_{1},$

where

$r_{1}=\beta
e_{i}^{T}[\mathbf{1}_{n}(UU^{T}\mathbf{1}_{n}-\mathbf{1}_{n})^{T}+(UU^{T}\mathbf{1}_{n}-\mathbf{1}_{n})\mathbf{1}_{n}^{T}]{d}_{1}.$

Further write ${d}_{\tau}=\sum_{k}x_{k,\tau}\mathbf{1}_{G_{k}}+w_{\tau},$ for
$x_{k,\tau}\in\mathbb{R},\;\tau=1,\dots,r,\;w_{\tau}\in\Gamma_{K}^{\perp}.$
Then we have

$\displaystyle e_{i}^{T}(2A+y\mathbf{1}_{n}^{T}+\mathbf{1}_{n}y^{T}){d}_{1}$
$\displaystyle=\sum_{k}x_{k,1}D_{k,l}(i)+e_{i}^{T}(2A+y\mathbf{1}_{n}^{T}+\mathbf{1}_{n}y^{T})w_{1}$
$\displaystyle\geqslant x_{1,1}D_{1,l}(i)-2\lambda\cdot
x_{l,1}+e_{i}^{T}(2A+y\mathbf{1}_{n}^{T}+\mathbf{1}_{n}y^{T})w_{1},$

since $x_{k,1}\geqslant 0,\;\forall k\neq 1,\;D_{l,l}=-2\lambda n_{l}.$ If we
let $X=[X_{1},\dots,X_{n}],\;X_{i}=\varepsilon_{i}+\mu_{l},\;\forall i\in
G_{l}.$ Then $\forall i,$

$\displaystyle|e_{i}^{T}(2A+y\mathbf{1}_{n}^{T}+\mathbf{1}_{n}y^{T})w_{1}|$
$\displaystyle=|2[Aw_{1}]_{i}-\sum_{k}\frac{2}{n_{k}}\mathbf{1}_{n_{k}}^{T}A_{G_{k},G_{k}}w_{G_{k}}|$
$\displaystyle=\left|2\sum_{k}\left[(\mu_{l}-\mu_{k})+(\varepsilon_{i}-\bar{\varepsilon}_{k})\right]^{T}\left[\sum_{j\in
G_{k}}\varepsilon_{j}(w_{1})_{j}\right]\right|$ $\displaystyle\leqslant
2K(\Theta_{\max}\sqrt{2K\log n}+\sqrt{n}(p+2\sqrt{pK\log n}+2K\log
n))\|w_{1}\|,$

with probability $\geqslant 1-C_{4}/n,$ for some constant $C_{4}>0.$ On the
other hand,

$\displaystyle\;\;\;\;\;\beta
e_{i}^{T}[\mathbf{1}_{n}(UU^{T}\mathbf{1}_{n}-\mathbf{1}_{n})^{T}+(UU^{T}\mathbf{1}_{n}-\mathbf{1}_{n})\mathbf{1}_{n}^{T}]a_{1}\mathbf{1}_{G_{1}}$
$\displaystyle=\beta
a_{1}e_{i}^{T}[\mathbf{1}_{n}\mathbf{1}_{n}^{T}((U^{\mathbf{a},\ast})\Delta^{T}+\Delta(U^{\mathbf{a},\ast})^{T})+((U^{\mathbf{a},\ast})\Delta^{T}+\Delta(U^{\mathbf{a},\ast})^{T})\mathbf{1}_{n}\mathbf{1}_{n}^{T}]\mathbf{1}_{G_{1}}$
$\displaystyle=\beta
a_{1}[\mathbf{1}_{n}^{T}(U^{\mathbf{a},\ast}\Delta^{T}+\Delta(U^{\mathbf{a},\ast})^{T})\mathbf{1}_{G_{1}}+n_{1}e_{i}^{T}(U^{\mathbf{a},\ast}\Delta^{T}+\Delta(U^{\mathbf{a},\ast})^{T})\mathbf{1}_{n}]$
$\displaystyle=\beta
a_{1}\sum_{\tau=1}^{r}[\mathbf{1}_{n}^{T}(u_{\tau}^{0}{d}_{\tau}^{T}+{d}_{\tau}(u_{\tau}^{0})^{T})\mathbf{1}_{G_{1}}+n_{1}e_{i}^{T}(u_{\tau}^{0}{d}_{\tau}^{T}+{d}_{\tau}(u_{\tau}^{0})^{T})\mathbf{1}_{n}]$
$\displaystyle=\beta
a_{1}[(\sum_{g=1}^{K}n_{g}\sum_{\tau=s(g-1)+1}^{sg}a_{g}{d}_{\tau}^{T}\mathbf{1}_{G_{1}})+a_{1}n_{1}\sum_{\tau=1}^{s}\mathbf{1}_{n}^{T}{d}_{\tau}$
$\displaystyle\;\;\;\;\;+n_{1}\sum_{\tau=s(l-1)+1}^{sl}a_{l}\mathbf{1}_{n}^{T}{d}_{\tau}+n_{1}(\sum_{g=1}^{K}n_{g}\sum_{\tau=s(g-1)+1}^{sg}a_{g}({d}_{\tau})_{i})]$
$\displaystyle\geqslant\beta
a_{1}[2ra_{1}n_{1}\sqrt{n}(-\|\Delta\|_{F})+rn_{1}\bar{a}\sqrt{n_{l}}(-\|\Delta\|_{F})$
$\displaystyle\;\;\;\;\;+n_{1}(a_{1}n_{1}(d_{1})_{i}+\bar{a}n_{l}\sum_{\tau=s(l-1)+1}^{sl}(\mathbf{d}_{\tau})_{i}))]$
$\displaystyle\geqslant\beta
a_{1,1}[3r\bar{a}n_{1}\sqrt{n}(-\|\Delta\|_{F})+a_{1}n_{1}^{2}(d_{1})_{i}-\bar{a}n_{1}n_{l}r\bar{a}].$

Similarly we can show that $|r_{1}|\leqslant\beta
r\bar{a}(3n^{2}\bar{a}^{2}+\bar{n}\|\Delta\|_{F}^{2}).$ Hence

$\displaystyle[\nabla f(U)]_{i,1}$ $\displaystyle\geqslant
2L\cdot(d_{1})_{i}+(a_{1,1}+x_{1,1})D_{1,l}(i)-2\lambda\|\Delta_{S^{c}}\|_{\infty}$
$\displaystyle-2K(\Theta_{\max}\sqrt{2K\log n}+\sqrt{n}(p+2\sqrt{pK\log
n}+2K\log n))\|w_{1}\|$ $\displaystyle+\beta
a_{1,1}[3r\bar{a}n_{1}\sqrt{n}(-\|\Delta\|_{F})+a_{1,1}n_{1}^{2}(d_{1})_{i}-\bar{a}n_{1}n_{l}r\bar{a}]+r_{1}$
$\displaystyle\geqslant(2L+\beta
a_{1,1}^{2}n_{1}^{2})(d_{1})_{i}+\frac{C_{1}\underline{a}\underline{n}\Theta_{\min}^{2}}{8},$

provided that $\Delta:$

$x_{1,1}\geqslant-3/4a_{1},$
$\|\Delta_{S(\mathbf{a})^{c}}\|_{\infty}\leqslant\frac{C_{1}\underline{a}\underline{n}\Theta_{\min}^{2}}{32\lambda}:=\epsilon_{f_{1}},\;\|\Delta\|_{F}\leqslant\epsilon_{1}.$

Here $\epsilon_{1}:=$

$\min\left\\{\frac{C_{1}\underline{a}\underline{n}\Theta_{\min}^{2}}{64K(\Theta_{\max}\sqrt{2K\log
n}+\sqrt{n}(p+2\sqrt{pK\log n}+2K\log
n))},\sqrt{\frac{C_{1}\underline{a}\underline{n}\Theta_{\min}^{2}}{96\beta
r\bar{a}\bar{n}}},\frac{C_{1}\underline{a}\underline{n}\Theta_{\min}^{2}}{96\beta
r\bar{a}^{2}\bar{n}\sqrt{n}}\right\\}.$

It is also sufficient to assume

$\|\Delta_{S(\mathbf{a})^{c}}\|_{\infty}\leq\epsilon_{f_{1}},$
$\|\Delta\|_{F}\leq\min\left\\{\epsilon_{1},\frac{3}{4}\underline{a}\sqrt{\underline{n}}\right\\}.$

And for $(L,\beta):$

$\frac{2(\lambda-L)}{\underline{n}}\leqslant\beta\leqslant\frac{C_{1}\underline{a}\underline{n}\Theta_{\min}^{2}}{12rn^{2}\bar{a}^{3}}.$

Thus,

$[U]_{i,1}-\alpha[\nabla
f(U)]_{i,1}\leqslant-\alpha\frac{C_{1}\underline{a}\underline{n}\Theta_{\min}^{2}}{8},\forall
i\in G_{l},l\neq 1.$

$\blacksquare$

######  Corollary 13.

Following the same conditions, suppose $\|\Delta\|_{F}\leqslant\epsilon_{c},$
where $\epsilon_{c}$ is defined in Lemma 19, then we have

$[V]_{i,\tau}\leqslant\max\\{[U]_{i,\tau}-\alpha\frac{C_{1}\underline{a}\underline{n}\Theta_{\min}^{2}}{8},\;0\\},\forall
i\in G_{l},s(k-1)<\tau\leqslant sk,\;l\neq k.$

where $V=\Pi_{\Omega}(W),\;W=U-\alpha\nabla f(U).$ Hence for at most
$I=\frac{1}{4\alpha\lambda}$ iterations, $V\in\mathcal{G}_{\mathbf{m}}.$

_Proof._ By Lemma 20 we have $\forall i\in G_{l},s(k-1)<\tau\leqslant
sk,\;l\neq k,$

$[V]_{i,1}=\frac{\sqrt{K}}{\sqrt{\|W\|}}\Pi_{+}([U]_{i,1}-\alpha[\nabla
f(U)]_{i,1})\leqslant\Pi_{+}([U]_{i,1}-\alpha[\nabla f(U)]_{i,1}).$

$\blacksquare$

######  Lemma 14.

If $\|U-U^{\ast}Q\|_{F}\leqslant\epsilon_{c},$ then

$\|V-U^{\mathbf{a},\ast}\|_{F}\leqslant\eta\|U-U^{\mathbf{a},\ast}\|_{F},$

where $\eta=1+\alpha R_{1},\;R_{1}=2C_{2}(n+\sqrt{p\log
n})+2/C_{1}\Theta^{2}_{\max}+2L+12n\beta.$

_Proof._ Suppose $\|U-U^{\ast}Q\|_{F}\leqslant\epsilon_{c},$ then from Lemma
20 we have

$\displaystyle\|V-U^{\mathbf{a},\ast}\|_{F}$
$\displaystyle=\|\Pi_{\Omega}(U-\alpha\nabla
f(U))-\Pi_{\Omega}(U^{\mathbf{a},\ast}-\alpha\nabla
f(U^{\mathbf{a},\ast}))\|_{F}$
$\displaystyle=\|\Pi_{\mathcal{C}(\Omega)}(U-\alpha\nabla
f(U))-\Pi_{\mathcal{C}(\Omega)}(U^{\mathbf{a},\ast}-\alpha\nabla
f(U^{\mathbf{a},\ast}))\|_{F}$ $\displaystyle\leqslant\|(U-\alpha\nabla
f(U))-(U^{\mathbf{a},\ast}-\alpha\nabla f(U^{\mathbf{a},\ast}))\|_{F}$
$\displaystyle\leqslant(1+tR_{1})\|U-U^{\mathbf{a},\ast}\|_{F},$

where $R_{1}=2C_{2}(n+\sqrt{p\log n})+2/C_{1}\Theta^{2}_{\max}+2L+12n\beta.$
$\blacksquare$

###  B.6  Linear convergence near optimal point (Phase 2)

Here our goal is to prove Theorem 10. The main idea is to show the local
strongly convexity at optimal points, which is shown in Lemma 15 and Lemma 17.
The rest lemmas provide details for calculating the concrete form of the local
strongly convexity radius around optimal points.

######  Lemma 15 (Linear convergence for every step).

If $U\in\mathcal{G}$, define $\Delta:=UQ^{T}-U^{\ast},$ and suppose

$\langle\nabla^{2}f(U^{\ast})[\Delta],\Delta\rangle\geqslant\tilde{\beta}\|\Delta\|_{F}^{2},$

for some $\tilde{\beta}>0.$ Then $\exists\gamma\in(0,1),\epsilon>0,$ s.t.,

$\|V-U^{\ast}Q\|_{F}\leqslant\gamma\|U-U^{\ast}Q\|_{F},$

where

$Q=\text{argmin}_{\tilde{Q}\in\mathcal{F}}\|U-U^{\ast}\tilde{Q}\|_{F},$

provided that $\|U-U^{*}Q\|_{F}<\epsilon,$ for some $\epsilon>0.$

_Proof._ By Lemma 19 we have

$\displaystyle\langle\nabla f(U),U\rangle$ $\displaystyle=\langle\nabla
f(U)(Q^{a})^{T},U(Q^{a})^{T}\rangle$ $\displaystyle=\langle\nabla
f(U(Q^{a})^{T}),U(Q^{a})^{T}\rangle$ $\displaystyle<-(\lambda-L)K/2,$

provided that $\|U-U^{\mathbf{a},\ast}\|_{F}<\epsilon_{c}$. This implies that

$\|U-\alpha\nabla f(U)\|_{F}>\|U\|_{F}.$

Then it can be shown (with details provided in Lemma 20) that

$\Pi_{\Omega}(U-\alpha\nabla f(U))=\Pi_{\mathcal{C}(\Omega)}(U-\alpha\nabla
f(U)),$

where $\mathcal{C}(\Omega)$ stands for the convex hall of $\Omega.$ It is
known that

$\|\Pi_{C}(v)-\Pi_{C}(u)\|\leqslant\|v-u\|,$

for any convex set $C$. Note
$\forall\tilde{Q}\in\mathcal{F},\;U^{\ast}\tilde{Q}$ is a stationary point for
$\Pi_{\Omega}$. Then we have

$\displaystyle\|V-U^{\ast}Q\|_{F}^{2}$
$\displaystyle=\|\Pi_{\Omega}(U-\alpha\nabla
f(U))-\Pi_{\Omega}(U^{\ast}Q-\alpha\nabla f(U^{\ast}Q))\|_{F}^{2}$
$\displaystyle=\|\Pi_{\mathcal{C}(\Omega)}(U-\alpha\nabla
f(U))-\Pi_{\mathcal{C}(\Omega)}(U^{\ast}Q-\alpha\nabla
f(U^{\ast}Q))\|_{F}^{2}$ $\displaystyle\leqslant\|(U-\alpha\nabla
f(U))-(U^{\ast}Q-\alpha\nabla f(U^{\ast}Q))\|_{F}^{2}$
$\displaystyle=\|(UQ^{T}-\alpha\nabla f(UQ^{T}))-(U^{\ast}-\alpha\nabla
f(U^{\ast}))\|_{F}^{2}$
$\displaystyle=\|UQ^{T}-U^{\ast}\|_{F}^{2}+\alpha^{2}\|\nabla
f(UQ^{T}))-\nabla f(U^{\ast}))\|_{F}^{2}$ $\displaystyle-2\alpha\langle\nabla
f(UQ^{T})-\nabla f(U^{\ast}),UQ^{T}-U^{\ast}\rangle,$

since $\nabla f(\tilde{U}\tilde{Q})=\nabla
f(\tilde{U})\tilde{Q},\;\forall\tilde{Q}\in\mathcal{O}_{r\times
r},\tilde{U}\in\mathbb{R}_{n\times r}.$ Note that
$\|UQ^{T}-U^{\ast}\|_{F}^{2}=\|U-U^{\ast}Q\|_{F}^{2},\;\|\nabla
f(UQ^{T})-\nabla f(U^{\ast}))\|_{F}\leqslant R_{1}\|UQ^{T}-U^{\ast}\|_{F},$
for some constant $R_{1},$ so we only need to analyze the last term. And by
MVT we have

$\langle\nabla f(UQ^{T})-\nabla
f(U^{\ast}),UQ^{T}-U^{\ast}\rangle=\int_{0}^{1}\langle\nabla^{2}f(U^{\ast}+\tau(UQ^{T}-U^{\ast}))[UQ^{T}-U^{\ast}],UQ^{T}-U^{\ast}\rangle
d\tau.$

Suppose

$\langle\nabla^{2}f(U^{\ast})[\Delta],\Delta\rangle\geqslant\tilde{\beta}\|\Delta\|_{F}^{2},$

for some $\tilde{\beta}>0.$ By continuity we have $\exists
0<\epsilon_{1}<\epsilon_{c},$

$\langle\nabla^{2}f(U)[\Delta],\Delta\rangle\geqslant\tilde{\beta}/2\|\Delta\|_{F}^{2},$

$\forall U$ s.t. $\|U-U^{\ast}\|_{F}<\epsilon_{1}.$ Now we set
$\|UQ^{T}-U^{\ast}\|_{F}<\epsilon_{1},$ then

$\displaystyle\langle\nabla f(UQ^{T})-\nabla
f(U^{\ast}),UQ^{T}-U^{\ast}\rangle$
$\displaystyle=\int_{0}^{1}\langle\nabla^{2}f(U^{\ast}+\tau(UQ^{T}-U^{\ast}))[UQ^{T}-U^{\ast}],UQ^{T}-U^{\ast}\rangle
d\tau$ $\displaystyle\geqslant\tilde{\beta}/2\|UQ^{T}-U^{\ast}\|_{F}^{2}.$

By combining them together, we have

$\displaystyle\|V-U^{\ast}Q\|_{F}^{2}$
$\displaystyle=\|UQ^{T}-U^{\ast}\|_{F}^{2}+\alpha^{2}\|\nabla
f(UQ^{T}))-\nabla f(U^{\ast}))\|_{F}^{2}$ $\displaystyle-2\alpha\langle\nabla
f(UQ^{T})-\nabla f(U^{\ast}),UQ^{T}-U^{\ast}\rangle$
$\displaystyle\leqslant\|UQ^{T}-U^{\ast}\|_{F}^{2}+R_{1}^{2}\alpha^{2}\|UQ^{T}-U^{\ast}\|_{F}^{2}-2\alpha\tilde{\beta}/2\|UQ^{T}-U^{\ast}\|_{F}^{2}$
$\displaystyle\leqslant(1-\alpha\tilde{\beta}/2)\|UQ^{T}-U^{\ast}\|_{F}^{2}$
$\displaystyle=\gamma^{2}\|U-U^{\ast}Q\|_{F}^{2},$

by choosing
$\alpha<\tilde{\beta}/(2R_{1}^{2}),\gamma^{2}=(1-\alpha\tilde{\beta}/2)$.
$\blacksquare$

######  Lemma 16 (Bounded over iterations).

Suppose Lemma 15 holds and $U^{0}\in\mathcal{G}$. Then given fixed
$\bar{\epsilon}:=\min\left\\{\epsilon_{1},\frac{3}{4}\underline{a}\sqrt{\underline{n}}\right\\},$
we have

$\|U^{t}-U^{\mathbf{a},\ast}\|_{F}\leqslant\bar{\epsilon},\;\forall k\geqslant
1,$

where

$U^{t+1}=\Pi_{\Omega}(U^{t}-\alpha\nabla f(U^{t})),\;\forall k\geqslant 1,$

given $\|U^{0}-U^{\mathbf{a},\ast}\|_{F}\leqslant\epsilon_{0},$ for some
$\epsilon_{0}>0.$

_Proof._ Recall that we define $Q^{t}\in\mathcal{F}$ as

$Q^{t}=\text{argmin}_{Q\in\mathcal{F}}\|U^{t}-U^{\ast}Q\|_{F}.$

Define $\epsilon_{0}:=(1-\gamma)/2\bar{\epsilon}$. Note that
$U^{0}\in\mathcal{G}$, then from Lemma 15 we have

$\|U^{1}-U^{\ast}Q^{0}\|_{F}\leqslant\gamma\|U^{0}-U^{\ast}Q^{0}\|_{F}.$

Note $U^{\ast}Q^{0}$ is the projection of $U^{0}$ on to $\mathcal{F}$ that
located on a sphere, which implies that

$\|U^{\ast}Q^{0}-U^{\mathbf{a},\ast}\|_{F}<2\|U^{0}-U^{\mathbf{a},\ast}\|_{F},$

given $\|U^{0}-U^{\mathbf{a},\ast}\|_{F}<1/2\|U^{\mathbf{a},\ast}\|_{F}.$
Thus,

$\displaystyle\|U^{1}-U^{\mathbf{a},\ast}\|_{F}$
$\displaystyle\leqslant\|U^{1}-U^{\ast}Q^{0}\|_{F}+\|U^{\ast}Q^{0}-U^{\mathbf{a},\ast}\|_{F}$
$\displaystyle\leqslant\gamma\|U^{0}-U^{\ast}Q^{0}\|_{F}+2\|U^{0}-U^{\mathbf{a},\ast}\|_{F}$
$\displaystyle\leqslant\gamma\|U^{0}-U^{\mathbf{a},\ast}\|_{F}+2\|U^{0}-U^{\mathbf{a},\ast}\|_{F}$
$\displaystyle\leqslant\gamma\epsilon_{0}+2\epsilon_{0}$
$\displaystyle\leqslant\bar{\epsilon}.$

Then we will finish the proof by induction. Suppose it holds for $l\leqslant
k,$ then $U^{l+1}\in\mathcal{G},$ from Lemma 15 we have

$\|U^{l+1}-U^{\ast}Q^{l+1}\|_{F}\leqslant\|U^{l+1}-U^{\ast}Q^{l}\|_{F}\leqslant\gamma\|U^{l}-U^{\ast}Q^{l}\|_{F},\
\;\forall l\leqslant k.$

Hence

$\displaystyle\|U^{t+1}-U^{\mathbf{a},\ast}\|_{F}$
$\displaystyle\leqslant\|U^{t+1}-U^{\ast}Q^{t}\|_{F}+\|U^{\ast}Q^{t}-U^{\ast}Q^{k-1}\|_{F}+\cdots+\|U^{\ast}Q^{0}-U^{\mathbf{a},\ast}\|_{F}$
$\displaystyle\leqslant
2\sum_{l=0}^{t}\|U^{l+1}-U^{\ast}Q^{l}\|_{F}+2\|U^{0}-U^{\mathbf{a},\ast}\|_{F}$
$\displaystyle\leqslant
2\sum_{l=0}^{t+1}\gamma^{l}\|U^{0}-U^{\mathbf{a},\ast}\|_{F}$
$\displaystyle\leqslant\bar{\epsilon}.$

$\blacksquare$

######  Lemma 17 (Positive definiteness w.r.t. $\Delta$).

If $U\in\mathcal{G},$ define $\Delta:=UQ^{T}-U^{\ast},$ then

$\langle\nabla^{2}f(U^{\ast})[\Delta],\Delta\rangle\geqslant\tilde{\beta}\|\Delta\|_{F}^{2},$

for some $\tilde{\beta}>0$, provided that $L\in((\sqrt{n}+\sqrt{p}+\sqrt{2\log
n})^{2},\lambda),\;\beta>-R_{2}/R_{3},$ for some $R_{2}\in\mathbb{R},R_{3}>0.$

_Proof._ $U\in\mathcal{G},Q\in\mathcal{F}\implies U=W_{1}+W_{2},$ with

$W_{1}:=\begin{bmatrix}d_{1,1}\mathbf{1}_{n_{1}},\dots,d_{1,s}\mathbf{1}_{n_{1}}&0&\cdots&0\\\
0&d_{2,1}\mathbf{1}_{n_{2}},\dots,d_{2,s}\mathbf{1}_{n_{2}}&\cdots&0\\\
\vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&d_{K,1}\mathbf{1}_{n_{K}},\dots,d_{K,s}\mathbf{1}_{n_{K}}\end{bmatrix},$
$W_{2}:=\begin{bmatrix}\textbf{w}_{1,1},\dots,\textbf{w}_{1,s}&0&\cdots&0\\\
0&\textbf{w}_{2,1},\dots,\textbf{w}_{2,s}&\cdots&0\\\
\vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&\textbf{w}_{K,1},\dots,\textbf{w}_{K,s}\end{bmatrix},$

for some $d_{k,i}\in\mathbb{R},\textbf{w}_{k,i}\in\mathbb{R}^{n_{k}},\;$
$\langle\textbf{w}_{k,i},\mathbf{1}_{n_{k}}\rangle=0,\forall k\in[K],i\in[s].$
Note that

$U^{\ast}Q:=\begin{bmatrix}\tilde{a}_{1,1}\mathbf{1}_{n_{1}},\dots,\tilde{a}_{1,s}\mathbf{1}_{n_{1}}&0&\cdots&0\\\
0&\tilde{a}_{2,1}\mathbf{1}_{n_{2}},\dots,\tilde{a}_{2,s}\mathbf{1}_{n_{2}}&\cdots&0\\\
\vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&\tilde{a}_{K,1}\mathbf{1}_{n_{K}},\dots,\tilde{a}_{K,s}\mathbf{1}_{n_{K}}\end{bmatrix},$

for some $\tilde{a}_{k,i}\in\mathbb{R},k\in[K],i\in[s].$ Recall

$Q=\text{argmin}_{\tilde{Q}\in\mathcal{F}}\|U-U^{\ast}\tilde{Q}\|_{F},$

then there exist $c_{k}>0,k\in[K],$ s.t.

$d_{k,i}=c_{k}\tilde{a}_{k,i},\;\forall k\in[K],i\in[s].$

Hence $W_{1}Q^{T}=[W_{1}^{1}|O_{n\times(r-K)}],$ where

$W_{1}^{1}:=\begin{bmatrix}\frac{c_{1}}{\sqrt{n_{1}}}\mathbf{1}_{n_{1}}&0&\cdots&0\\\
0&\frac{c_{2}}{\sqrt{n_{2}}}\mathbf{1}_{n_{2}}&\cdots&0\\\
\vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&\frac{c_{K}}{\sqrt{n_{K}}}\mathbf{1}_{n_{K}}\end{bmatrix}.$

And $W_{2}Q^{T}=[W_{2}^{1}|W_{2}^{2}],$ where

$W_{2}^{1}:=\begin{bmatrix}\textbf{w}_{1}&0&\cdots&0\\\
0&\textbf{w}_{2}&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&\textbf{w}_{K}\end{bmatrix},$

for some $\textbf{w}_{k}\in\mathbb{R}^{n_{k}},k\in[K].$ And every column of
$W_{2}^{2}$ belongs to the space
$\Gamma_{K}:=\text{span}\\{\mathbf{1}_{G_{k}}:k\in[K]\\}^{\perp},$ which is
the orthogonal complement of the linear subspace of $\mathbb{R}^{n}$ spanned
by the vectors $\mathbf{1}_{G_{1}},\dots,\mathbf{1}_{G_{K}}.$ i.e.,

$UQ^{T}-U^{\ast}=[W_{1}^{1}+W_{2}^{1}-U^{\ast}_{1}|W_{2}^{2}].$

Denote
$\Delta_{1}=W_{1}^{1}+W_{2}^{1}-U^{\ast}_{1},\;\Delta_{2}=W_{2}^{2},\;\Delta=[\Delta_{1}|\Delta_{2}].$
By calculating the Hessian at $U^{\ast}$ we get

$\displaystyle\langle\nabla^{2}f(U^{\ast})[\Delta],\Delta\rangle$
$\displaystyle=\langle
2(L\cdot\text{Id}_{n}+A)+y\mathbf{1}_{n}^{T}+\mathbf{1}_{n}y^{T},\Delta\Delta^{T}\rangle+\beta\|\Delta(U^{\ast})^{T}\mathbf{1}_{n}+U^{\ast}\Delta^{T}\mathbf{1}_{n}\|_{F}^{2}$
$\displaystyle=\langle
2(L\cdot\text{Id}_{n}+A)+y\mathbf{1}_{n}^{T}+\mathbf{1}_{n}y^{T},\Delta_{2}\Delta_{2}^{T}\rangle$
$\displaystyle+\langle
2(L\cdot\text{Id}_{n}+A)+y\mathbf{1}_{n}^{T}+\mathbf{1}_{n}y^{T},\Delta_{1}\Delta_{1}^{T}\rangle+\beta\|\Delta_{1}(U^{\ast}_{1})^{T}\mathbf{1}_{n}+U^{\ast}_{1}\Delta_{1}^{T}\mathbf{1}_{n}\|_{F}^{2}$
$\displaystyle=2\langle(L\cdot\text{Id}_{n}+A),\Delta_{2}\Delta_{2}^{T}\rangle+\langle
2(L\cdot\text{Id}_{n}+A)+y\mathbf{1}_{n}^{T}+\mathbf{1}_{n}y^{T},\Delta_{1}\Delta_{1}^{T}\rangle$
$\displaystyle+\beta\|\Delta_{1}(U^{\ast})^{T}\mathbf{1}_{n}+U^{\ast}\Delta_{1}^{T}\mathbf{1}_{n}\|_{F}^{2}$
$\displaystyle\geqslant 2[L-(\sqrt{n}+\sqrt{p}+\sqrt{2\log
n})^{2}]\|\Delta_{2}\|_{F}^{2}+R_{2}\|\Delta_{1}\|_{F}^{2}+\beta
R_{3}\|\Delta_{1}\|_{F}^{2}$
$\displaystyle\geqslant\min\\{2[L-(\sqrt{n}+\sqrt{p}+\sqrt{2\log
n})^{2}],R_{2}+\beta R_{3}\\}\|\Delta\|_{F}^{2}$
$\displaystyle:=\tilde{\beta}\|\Delta\|_{F}^{2},$

with probability $\geqslant 1-1/n$ (by Proposition 7), where
$R_{2}\in\mathbb{R},R_{3}>0$ is certain constant ($R_{3}>0$ follows from Lemma
18 below). By choosing $L\in((\sqrt{n}+\sqrt{p}+\sqrt{2\log
n})^{2},\lambda),\;\beta>-R_{2}/R_{3},$ we will get $\tilde{\beta}>0.$
$\blacksquare$

######  Lemma 18 (Positive definiteness w.r.t. $\Delta_{1}$).

Suppose $(U^{\ast}_{1})_{ij}=0\implies(\Delta_{1})_{ij}=0.$ Then

$\|\Delta_{1}(U^{\ast}_{1})^{T}\mathbf{1}_{n}+U^{\ast}_{1}\Delta_{1}^{T}\mathbf{1}_{n}\|_{F}^{2}\geqslant
R_{3}\|\Delta_{1}\|_{F}^{2},$

for some $R_{3}>0.$

_Proof._ Recall

$U_{1}^{\ast}=\begin{bmatrix}\frac{1}{\sqrt{n_{1}}}\mathbf{1}_{n_{1}}&0&\cdots&0\\\
0&\frac{1}{\sqrt{n_{2}}}\mathbf{1}_{n_{2}}&\cdots&0\\\
\vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&\frac{1}{\sqrt{n_{K}}}\mathbf{1}_{n_{K}}\end{bmatrix}.$

Now suppose

$\Delta_{1}=\begin{bmatrix}d_{1}&0&\cdots&0\\\ 0&d_{2}&\cdots&0\\\
\vdots&\vdots&\ddots&\vdots\\\ 0&0&\cdots&\mathbf{d}_{K}\end{bmatrix},$

where $d_{k}\in\mathbb{R}^{n_{k}}.$ Then we have

$\Delta_{1}(U^{\ast}_{1})^{T}\mathbf{1}_{n}+U^{\ast}_{1}\Delta_{1}^{T}\mathbf{1}_{n}=\begin{bmatrix}\tilde{d}_{1}\\\
\tilde{d}_{2}\\\ \vdots\\\ \tilde{d}_{K}\end{bmatrix},$

where

$\tilde{d}_{k}=\sqrt{n_{k}}d_{k}+\frac{1}{\sqrt{n_{k}}}\mathbf{1}_{n_{k}}\mathbf{1}_{n_{k}}^{T}{d}_{k}.$

Note that

$\|\tilde{d}_{k}\|^{2}\geqslant\lambda^{2}_{\min}\left(\sqrt{n_{k}}\text{Id}_{n_{k}}+\frac{1}{\sqrt{n_{k}}}\mathbf{1}_{n_{k}}\mathbf{1}_{n_{k}}^{T}\right)\|d_{k}\|^{2}\geqslant
n_{k}\|d_{k}\|^{2}.$

Hence

$\displaystyle\|\Delta_{1}(U^{\ast}_{1})^{T}\mathbf{1}_{n}+U^{\ast}_{1}\Delta_{1}^{T}\mathbf{1}_{n}\|_{F}^{2}$
$\displaystyle=\sum_{k}\|\tilde{d}_{k}\|^{2}$
$\displaystyle\geqslant\min_{k}\\{n_{k}\\}\sum_{k}\|d_{k}\|^{2}$
$\displaystyle=R_{3}\|\Delta_{1}\|_{F}^{2},$

where $R_{3}=\min_{k}\\{n_{k}\\}>0.$ $\blacksquare$

######  Lemma 19.

$\exists 0<\epsilon_{c}<1,$ s.t.
$\|U-U^{*}Q\|_{F}\leqslant\epsilon_{c}\implies$

$\langle\nabla f(U),U\rangle<-(\lambda-L)K/2<0.$

_Proof._ Let $[U]_{S}$ to be the vector of entries of $U$ that located at the
position of positive entries of $U$. Hence

$\nabla f(U)=[\nabla f(U)]_{S}+[\nabla f(U)]_{S^{c}}.$

Note that

$\|[\nabla
f(U^{*}))]_{S}\|_{F}=\|-2(\lambda-L)[U^{*}]_{S}\|_{F}=2(\lambda-L)\sqrt{K},$
$\|[\nabla f(U^{*}))]_{S^{c}}\|_{F}\leqslant n\Theta^{2}_{\max}.$

Then by Lemma 22 we have

$\displaystyle|\langle\nabla f(U),U\rangle-\langle\nabla
f(U^{*}Q),U^{*}Q\rangle|$ $\displaystyle=|\langle\nabla f(U)-\nabla
f(U^{*}Q),U\rangle+\langle\nabla f(U^{*}Q),U-U^{*}Q\rangle|$
$\displaystyle\leqslant\|\nabla f(U)-\nabla f(U^{*}Q)\|_{F}\|U\|_{F}+\|\nabla
f(U^{*}Q)\|_{F}\|U-U^{*}Q\|_{F}$ $\displaystyle=\|\nabla f(U)-\nabla
f(U^{*}Q)\|_{F}\|U\|_{F}+\|\nabla f(U^{*})Q\|_{F}\|U-U^{*}Q\|_{F}$
$\displaystyle\leqslant[\sqrt{K}R_{1}+2(\lambda-L)\sqrt{K}+n\Theta^{2}_{\max}]\|U-U^{*}Q\|_{F}.$

Finally we have

$\langle\nabla f(U),U\rangle<-(\lambda-L)K/2<0,$

given
$\epsilon_{c}\leqslant\frac{(\lambda-L)\sqrt{K}}{2[\sqrt{K}R_{1}+2(\lambda-L)\sqrt{K}+n\Theta^{2}_{\max}]}.$
Here $R_{1}=2C_{2}(n+\sqrt{p\log n})+2L+2/C_{1}\Theta^{2}_{\max}+12n\beta.$
$\blacksquare$

######  Lemma 20 (Convexity).

If $U\geqslant 0$, $t>0,\langle\nabla f(U),U\rangle<0,$ then

$\Pi_{\Omega}(U-\alpha\nabla f(U))=\Pi_{\mathcal{C}(\Omega)}(U-\alpha\nabla
f(U)).$

_Proof._ Note the fact

$\Pi_{\Omega}(U)=\Pi_{\mathbb{S}^{nr-1}(\sqrt{K})}(\Pi_{+}(U)),\;\Pi_{\mathcal{C}(\Omega)}(U)=\Pi_{\mathcal{B}^{nr}(\sqrt{K})}(\Pi_{+}(U)),$

where $\Pi_{+}(U)$ stands for the projection of $U$ on to the space of
matrices with positive entries, $\mathbb{S}^{nr-1}(\sqrt{K})$ stands for the
sphere with radius $\sqrt{K}$, $\mathcal{B}^{nr}(\sqrt{K})$ stands for the
ball with radius $\sqrt{K}$. Thus we only need to show that

$\|\Pi_{+}(U-\alpha\nabla f(U))\|_{F}^{2}\geqslant K.$

Define $\tilde{U}_{i,j}=t[\nabla f(U)]_{i,j},$ if $[U-\alpha\nabla
f(U)]_{i,j}>0.$ And $\tilde{U}_{i,j}=U_{i,j}$ otherwise. Then we have

$U-\tilde{U}=\Pi_{+}(U-\tilde{U})=\Pi_{+}(U-\alpha\nabla f(U)).$

Moreover, if $[U-\alpha\nabla f(U)]_{i,j}\leqslant 0,$ then

$[\alpha\nabla f(U)]_{i,j}\geqslant U_{i,j}\geqslant 0.$

Hence

$\langle U,\tilde{U}\rangle\leqslant t\langle\nabla f(U),U\rangle<0,$

which implies that

$\|\Pi_{+}(U-\alpha\nabla
f(U))\|_{F}^{2}=\|U-\tilde{U}\|_{F}^{2}>\|U\|_{F}^{2}=K.$

$\blacksquare$

######  Lemma 21.

Let $w\in\mathbb{R}^{n},$ $S$ to be the set of non-zero positions for $w$. If
$S\subseteq G_{k},$ for some $k\in[K]$, then

$-2\lambda\|w\|_{F}^{2}\leqslant\langle
2A+y\mathbf{1}_{n}^{T}+\mathbf{1}_{n}y^{T},ww^{T}\rangle\leqslant
2C_{2}(n+\sqrt{p\log n})\|w\|_{F}^{2}.$

Hence we can choose $R_{2}=-2(\lambda-L)$.

_Proof._ From Proposition 8 we have

$W_{n}:=\lambda\text{Id}+\frac{1}{2}(y\mathbf{1}_{n}^{T}+\mathbf{1}_{n}y^{T})-B+A\succ
0.$

Hence

$\langle\frac{1}{2}(y\mathbf{1}_{n}^{T}+\mathbf{1}_{n}y^{T})+A,ww^{T}\rangle\geqslant\langle
B-\lambda\text{Id},ww^{T}\rangle.$

Note that $S\subseteq G_{k}\implies\langle B,ww^{T}\rangle=0.$ Thus,

$\langle
2A+y\mathbf{1}_{n}^{T}+\mathbf{1}_{n}y^{T},ww^{T}\rangle\geqslant-2\lambda\|w\|_{F}^{2}.$

On the other hand,

$\displaystyle\langle\frac{1}{2}(y\mathbf{1}_{n}^{T}+\mathbf{1}_{n}y^{T})+A,ww^{T}\rangle$
$\displaystyle=\langle W+B-\lambda\text{Id},ww^{T}\rangle$
$\displaystyle\leqslant\|W\|_{op}\|w\|_{F}^{2}-\lambda\|w\|_{F}^{2}$
$\displaystyle\leqslant C_{2}(n+\sqrt{p\log n})\|w\|_{F}^{2}.$

Since $\|W\|_{op}\leqslant 2C_{2}(n+\sqrt{p\log n})$ from Proposition 8.
Furthermore, if $S\not\subseteq G_{k},\;\forall k\in[K],$ then

$\displaystyle\langle\frac{1}{2}(y\mathbf{1}_{n}^{T}+\mathbf{1}_{n}y^{T})+A,ww^{T}\rangle$
$\displaystyle=\langle W+B-\lambda\text{Id},ww^{T}\rangle$
$\displaystyle\leqslant C_{2}(n+\sqrt{p\log
n})\|w\|_{F}^{2}+1/2\max_{ij}B_{ij}\|w\|_{F}^{2}.$

By the Proof of Lemma VI.1 from cut-off paper, we can derive

$B_{ij}\leqslant 2/C_{1}\Theta^{2}_{\max},\;\forall i,j\in[n].$

where universal constant $C_{1}\in(0,1),$ $\Theta_{\max}$ is the maximum of
separation. Thus,

$\langle\frac{1}{2}(y\mathbf{1}_{n}^{T}+\mathbf{1}_{n}y^{T})+A,ww^{T}\rangle\leqslant(C_{2}(n+\sqrt{p\log
n})+1/C_{1}\Theta^{2}_{\max})\|w\|_{F}^{2}.$

$\blacksquare$

######  Lemma 22.

If $\|U-U^{\mathbf{a},\ast}\|_{F}\leqslant 1,\;U\in\mathcal{G}_{\mathbf{m}}$,
then

$\|\nabla f(U)-\nabla f(U^{\mathbf{a},\ast})\|_{F}\leqslant
R_{1}\|U-U^{\mathbf{a},\ast}\|_{F},$

where $R_{1}=2C_{2}(n+\sqrt{p\log n})+2L+12n\beta.$ This also holds if we
replace $U^{\mathbf{a},\ast}$ by $U^{\ast}Q$ in Lemma 15. Furthermore, if
$U\not\in\mathcal{G}_{\mathbf{m}},$ $R_{1}=2C_{2}(n+\sqrt{p\log
n})+2/C_{1}\Theta^{2}_{\max}+2L+12n\beta,$ where $\Theta^{2}_{\max}$ is the
maximum squared separation.

_Proof._ By definition we have $\nabla f(U)-\nabla f(U^{\mathbf{a},\ast})=$

$(2A+2L\cdot\text{Id}+y\mathbf{1}_{n}^{T}+\mathbf{1}_{n}y^{T})(U-U^{\mathbf{a},\ast})+\beta[\mathbf{1}_{n}\mathbf{1}_{n}^{T}(UU^{T}-U^{\mathbf{a},\ast}(U^{\mathbf{a},\ast})^{T})+(UU^{T}-U^{\mathbf{a},\ast}(U^{\mathbf{a},\ast})^{T})\mathbf{1}_{n}\mathbf{1}_{n}^{T}]U.$

Note that

$\displaystyle\|U\|_{op}$
$\displaystyle\leqslant\|U-U^{\mathbf{a},\ast}\|_{op}+\|U^{\mathbf{a},\ast}\|_{op}$
$\displaystyle\leqslant\|U-U^{\mathbf{a},\ast}\|_{F}+\|U^{\mathbf{a},\ast}\|_{op}$
$\displaystyle\leqslant 1+1$ $\displaystyle=2.$

Hence

$\displaystyle\|(UU^{T}-U^{\mathbf{a},\ast}(U^{\mathbf{a},\ast})^{T})\mathbf{1}_{n}\mathbf{1}_{n}^{T}U\|_{F}$
$\displaystyle\leqslant\|(UU^{T}-U^{\mathbf{a},\ast}(U^{\mathbf{a},\ast})^{T})\|_{op}\|U\|_{op}\|\mathbf{1}_{n}\mathbf{1}_{n}^{T}\|_{F}$
$\displaystyle\leqslant
2n\|U(U-U^{\mathbf{a},\ast})^{T}+(U-U^{\mathbf{a},\ast})(U^{\mathbf{a},\ast})^{T}\|_{F}$
$\displaystyle\leqslant 6n\|U-U^{\mathbf{a},\ast}\|_{F}.$

Note every column of $U-U^{\mathbf{a},\ast}$ satisfies the condition of $w$ in
Lemma 21. Thus,

$\displaystyle\|\nabla f(U)-\nabla f(U^{\mathbf{a},\ast})\|_{F}$
$\displaystyle\leqslant(\|(2A+y\mathbf{1}_{n}^{T}+\mathbf{1}_{n}y^{T})\|_{F}+2L)\|(U-U^{\mathbf{a},\ast})\|_{F}+2\beta\|(UU^{T}-U^{\mathbf{a},\ast}(U^{\mathbf{a},\ast})^{T})\|_{op}\|U\|_{op}\|\mathbf{1}_{n}\mathbf{1}_{n}^{T}\|_{F}$
$\displaystyle\leqslant(2C_{2}(n+\sqrt{p\log
n})+2L)\|U-U^{\mathbf{a},\ast}\|_{F}+12n\beta\|U-U^{\mathbf{a},\ast}\|_{F}.$

$\blacksquare$

######  Lemma 23.

$\exists 0<\epsilon_{1}<1,$ s.t.

$\langle\nabla^{2}f(U)[\Delta],\Delta\rangle\geqslant\tilde{\beta}/2\|\Delta\|_{F}^{2},$

$\forall U$ s.t. $\|U-U^{\ast}\|_{F}<\epsilon_{1}$ and all the $\Delta$ of the
form in Proof of Lemma 15.

_Proof._ By calculating the Hessian of $f$ at $U^{\ast}$ and $U$ we get

$\langle\nabla^{2}f(U^{\ast})[\Delta],\Delta\rangle=\langle
2(L\cdot\text{Id}_{n}+A)+y\mathbf{1}_{n}^{T}+\mathbf{1}_{n}y^{T},\Delta\Delta^{T}\rangle+\beta\|\Delta(U^{\ast})^{T}\mathbf{1}_{n}+U^{\ast}\Delta^{T}\mathbf{1}_{n}\|_{F}^{2},$
$\displaystyle\langle\nabla^{2}f(U)[\Delta],\Delta\rangle$
$\displaystyle=\langle
2(L\cdot\text{Id}_{n}+A)+y\mathbf{1}_{n}^{T}+\mathbf{1}_{n}y^{T},\Delta\Delta^{T}\rangle+\beta\|\Delta(U)^{T}\mathbf{1}_{n}+U\Delta^{T}\mathbf{1}_{n}\|_{F}^{2}$
$\displaystyle+\beta\langle\mathbf{1}_{n}\mathbf{1}_{n}^{T}(UU^{T}-(U^{\ast})(U^{\ast})^{T})+(UU^{T}-(U^{\ast})(U^{\ast})^{T})\mathbf{1}_{n}\mathbf{1}_{n}^{T},\Delta\Delta^{T}\rangle,$

Then we have

$\displaystyle|\langle\nabla^{2}f(U^{\ast})[\Delta],\Delta\rangle-\langle\nabla^{2}f(U)[\Delta],\Delta\rangle|$
$\displaystyle\leqslant\beta(\|\Delta(U)^{T}\mathbf{1}_{n}+U\Delta^{T}\mathbf{1}_{n}\|_{F}^{2}-\|\Delta(U^{\ast})^{T}\mathbf{1}_{n}+U^{\ast}\Delta^{T}\mathbf{1}_{n}\|_{F}^{2})$
$\displaystyle+\beta|\langle\mathbf{1}_{n}\mathbf{1}_{n}^{T}(UU^{T}-(U^{\ast})(U^{\ast})^{T})+(UU^{T}-(U^{\ast})(U^{\ast})^{T})\mathbf{1}_{n}\mathbf{1}_{n}^{T},\Delta\Delta^{T}\rangle|$
$\displaystyle\leqslant\beta(12n\|U-U^{\ast}\|_{F}+6n\|U-U^{\ast}\|_{F})$
$\displaystyle=18\beta n\|U-U^{\ast}\|_{F}.$

Hence

$\langle\nabla^{2}f(U)[\Delta],\Delta\rangle\geqslant\tilde{\beta}/2\|\Delta\|_{F}^{2},$

provided that $\epsilon_{1}\leqslant\frac{\tilde{\beta}}{36\beta n}.$
$\blacksquare$