# Solving graph equipartition SDPs on an algebraic variety

Tianyun Tang , Kim-Chuan Toh Department of Mathematics, National University of
Singapore, Singapore 119076 (ttang@u.nus.edu). Department of Mathematics, and
Institute of Operations Research and Analytics, National University of
Singapore, Singapore 119076 (mattohkc@nus.edu.sg). The research of this author
is supported by the Ministry of Education, Singapore, under its Academic
Research Fund Tier 3 grant call (MOE-2019-T3-1-010).

###### Abstract

Semidefinite programs are generally challenging to solve due to their high
dimensionality. Burer and Monteiro developed a non-convex approach to solve
linear SDP problems by applying its low rank property. Their approach is fast
because they used factorization to reduce the problem size. In this paper, we
focus on solving the SDP relaxation of a graph equipartition problem, which
involves an additional semidefinite upper bound constraint over the
traditional linear SDP. By applying the factorization approach, we get a non-
convex problem with an additional non-smooth spectral inequality constraint.
We discuss when the non-convex problem is equivalent to the original SDP, and
when a second order stationary point of the non-convex problem is also a
global minimum. Our results generalize previous works on smooth non-convex
factorization approaches for linear SDP to the non-smooth case. Moreover, the
constraints of the non-convex problem involve an algebraic variety with some
conducive properties that allow us to use Riemannian optimization techniques
and non-convex augmented Lagrangian method to solve the SDP problem very
efficiently with certified global optimality.

keywords: Graph equipartition, Burer and Monteiro method, low rank SDP,
algebraic variety, Riemannian optimization, augmented Lagrangian method

##  1 Introduction

###  1.1 Literature Review

Many combinatorial optimization problems such as max-cut and maximum stable
set problems are NP-hard, which means that it is almost impossible to solve
them exactly in polynomial time. These problems are often solved approximately
via their semidefinite programming (SDP) relaxations which typically result in
linear SDP problems of the following form:

${\rm SDP:\ }\min_{X\in\mathbb{S}^{n}}\left\\{\left\langle C,X\right\rangle:\
\mathcal{A}\left(X\right)=b,\ X\succeq 0\right\\},$ (1)

where $b\in\mathbb{R}^{m}$ and $C\in\mathbb{S}^{n}$ are given data, and
$\mathcal{A}:\mathbb{S}^{n}\rightarrow\mathbb{R}^{m}$ is a given linear
mapping (see [16, 20, 31] for examples of SDP relaxation).

The above linear SDP problem (1) is a convex programming problem which can be
solved by various well developed solvers such as SDPT3 [36, 37], MOSEK [25],
and SDPNAL [42, 44]. Those solvers are efficient and accurate if $n$ is
moderate (say less than $2000$) and in the case of SDPT3 and MOSEK, $m$ is
also not too large (say less than $20000$). However, if $n$ is large, the SDP
will be challenging to solve by the aforementioned solvers because they will
incur excessive computing cost and memory usage due to the high dimensionality
in $n$. To overcome this difficulty, Burer and Monteiro [9, 10] proposed to
solve the linear SDP (1) by applying the factorization $X=RR^{\top}$ to
reformulate it as the following non-convex problem:

${\rm SDPLR:\ }\min\left\\{\big{\langle}C,\,RR^{\top}\big{\rangle}:\
\mathcal{A}\big{(}RR^{\top}\big{)}=b,\ R\in\mathbb{R}^{n\times r}\right\\}.$
(2)

Due to the low rank property (see, e.g., [26, 14, 1]) of the SDP (1), which
states that if the feasible region of (1) is compact, then it will have an
optimal solution of rank $\leq\sqrt{2m}.$ So (1) is equivalent to (2) when
$r\geq\sqrt{2m}$. Burer and Monteiro used an augmented Lagrangian method to
solve (2) and its high efficiency was verified in numerical comparison to
other algorithms such as the spectral bundle [18] and interior point methods
[3].

Since (2) is a non-convex problem, it may have spurious local minima. As far
as we know, the first foundational and thorough analysis about the global
optimality of (2) is done in a series of works [4, 6, 7, 29], by Boumal et al.
They proved that under LICQ-like regularity assumptions on the constraints of
(2), if $r\geq\sqrt{2m}$, then for almost all $C\in\mathbb{S}^{n}$ except for
a set of measure zero, every second order stationary point of (2) is a global
optimal solution. This result explains why the BM approach can almost always
find the global minimum in practice. They also did a smoothed analysis for
problem (2) and its penalty version, which was useful in choosing the stopping
criterion for the augmented Lagrangian subproblems. Moreover, with their
regularity assumption, the set $\left\\{R\in\mathbb{R}^{n\times r}:\
\mathcal{A}\left(RR^{\top}\right)=b\right\\}$ is a Riemannian manifold
embedded in $\mathbb{R}^{n\times r},$ so they can use an optimization
algorithm on manifold to solve (2). They also developed a widely used toolbox
called manopt [8] to handle optimization problems on a variety of manifolds.

One specific example where Riemannian algorithms are expected to be efficient
for solving the problem (2) is the max-cut problem for which
$\mathcal{A}(X)={\rm diag}(X)$, and the underlying manifold is the oblique
manifold:

${\rm OB}_{n,r}:=\left\\{R\in\mathbb{R}^{n\times r}:\ {\rm
diag}\big{(}RR^{\top}\big{)}=e\right\\},$

whose tangent space and retraction mapping have a simple formula. In this
case, Boumal et al. [6] used manopt to solve (2) and verified its high
efficiency as compared to Burer and Monteiro’s augmented Lagrangian method and
other convex algorithms.

The research by Boumal et al. in [4, 6, 7, 29] was a milestone in studying the
global optimality of the non-convex factorized formulation of a linear SDP
problem by combining the geometric and low-rank property of the SDP problem.
Their results have since been generalized by Cifuentes in [11] to linear SDP
problems with multiple positive semidefinite block variables and inequality
constraints. One of our goals in this paper is to generalize the theory of
Boumal et al. to a class of SDP problems with an additional semidefinite upper
bound constraint, where the feasible sets of the non-convex factorized models
are no longer smooth manifolds.

###  1.2 Our Contribution

In this paper, we consider the graph equipartition problem, that is, given a
graph $G$ of size $n$, we want to partition it into $k$ parts of equal size
$q:=n/k$, such that the total number of edges between different parts is
minimized. Its SDP relaxation is given as follows (see [31, 21]):

${\rm GEP:}\min\left\\{\left\langle L,Y\right\rangle:\ {\rm diag}(Y)=e,\
Ye=qe,\ qI\succeq Y\succeq 0\right\\},$ (3)

where $L$ is the Laplacian of the graph and $e\in\mathbb{R}^{n}$ is the vector
of all ones. By changing the variable
$X=\left(Y-\frac{1}{k}ee^{\top}\right)\frac{k}{k-1}$, problem (3) is
equivalent to the following SDP problem,

${\rm GEP1:}\min\left\\{\left\langle L,X\right\rangle:\ {\rm diag}(X)=e,\
Xe=0,\ X\succeq 0,\ \frac{n}{k-1}I\succeq X\right\\}.$ (4)

Note that for $k=2$, we can remove the redundant constraint
$\frac{n}{k-1}I\succeq X,$ and it becomes a standard linear SDP111Note that
one may replace $Xe=0$ by $\left\langle X,ee^{\top}\right\rangle=0$ since
$X\succeq 0$.. To solve (4), we consider the following more general problem:

${\rm SDP1:}\ \min\left\\{\left\langle C,X\right\rangle:\
\mathcal{A}\left(X\right)=b,\ Xe=0,\
\big{\langle}X,\,ee^{\top}\big{\rangle}=0,\ \alpha I\succeq X\succeq
0\right\\},$ (5)

where $\mathcal{A}:\ \mathbb{S}^{n}\rightarrow\mathbb{R}^{m-1}$ is a general
linear mapping and $\alpha>0.$ Note that we add the redundant constraint
$Xe=0$ to facilitate our theoretical analysis later. For later usage, we
define

$\Omega_{\alpha}^{1}:=\left\\{X\in\mathbb{S}^{n}:\ \alpha I\succeq X\succeq
0\right\\}.$ (5a)

The low-rank factorization model of (5) is as follows:

${\rm SDPLR1:\
}\min\left\\{\frac{1}{2}\big{\langle}C,\,RR^{\top}\big{\rangle}:R\in\mathbb{R}^{n\times
r},\mathcal{A}\big{(}RR^{\top}\big{)}=b,\ R^{\top}e=0\
,\|R\|_{2}\leq\sqrt{\alpha}\right\\},$ (6)

where $\|R\|_{2}$ is the spectral norm of $R$. We add the factor $\frac{1}{2}$
in the objective function for convenience. Note that the main difference
between (1),(2) and (5),(6) is that we have an extra inequality constraint
$\|R\|_{2}\leq\sqrt{\alpha}$. We will prove that (6) is equivalent to (5) if
$r$ is larger than some given bound, which is a generalization of Burer and
Monteiro’s result. We will also prove that under the constraint nondegeneracy
condition (see [34]), for almost all $C\in\mathbb{S}^{n}$ except for a set of
measure zero, any second order stationary point of (6) is also a global
optimal solution. This is a generalization of the results in [6, 7, 29].

If we choose $\mathcal{A}(\cdot)={\rm diag}(\cdot),$ $b=e$ and $C=L$, then (6)
is the low-rank factorization of (4) for $\alpha={\frac{n}{k-1}}$, and we get

$\min\left\\{\frac{1}{2}\big{\langle}C,\,RR^{\top}\big{\rangle}:R\in\mathbb{R}^{n\times
r},{\rm diag}\big{(}RR^{\top}\big{)}=e,\ R^{\top}e=0\
,\|R\|_{2}\leq\sqrt{\alpha}\right\\},$ (7)

whose feasible region is the intersection of the following two sets:

$\mathcal{B}_{n,r}:=\left\\{R\in\mathbb{R}^{n\times r}:\ {\rm
diag}\big{(}RR^{\top}\big{)}=e,\ R^{\top}e=0\right\\},$ (8)
$\Omega_{\alpha}^{2}:=\left\\{R\in\mathbb{R}^{n\times r}:\
\|R\|_{2}\leq\sqrt{\alpha}\right\\}.$ (9)

Note that if $X=RR^{\top}$, then $X\in\Omega_{\alpha}^{1}$ if and only if
$R\in\Omega_{\alpha}^{2}.$ The set $\mathcal{B}_{n,r}$ is obtained by adding
the constraint $R^{\top}e=0$ to the oblique manifold ${\rm OB}_{n,r}$. Just as
each element of ${\rm OB}_{n,r}$ can be viewed as the Cartesian product of $n$
unit vectors, which are freely chosen, each element of $\mathcal{B}_{n,r}$ can
be considered as the Cartesian product of $n$ unit vectors, but they are
constrained to have their sum equal to zero. Therefore, those vectors are no
longer independent from each other and $\mathcal{B}_{n,r}$ is not as simple as
${\rm OB}_{n,r}.$ In fact, later we will show that $\mathcal{B}_{n,r}$ is not
even a manifold, but an algebraic variety defined as the set of common zeros
of a system of polynomials (see [12] chapter 3 definition 1). However, even
though it is not a manifold, we can still derive an explicit formula to
describe the tangent cones and second order tangent sets of its singular
points. More importantly, we find that one possible retraction mapping of
$\mathcal{B}_{n,r}$ is related to the well known geometric median problem,
which means that we can compute the retraction efficiently. These findings are
crucial for our algorithmic design. We will design an algorithm that equips
the underlying Riemannian optimization method with a strategy to check the
global optimality of a singular point and escape any non-optimal singular
point efficiently. To solve (7) with the spectral upper bound, we use an
augmented Lagrangian method on the algebraic variety $\mathcal{B}_{n,r}$, that
is, we keep the primal iterations on $\mathcal{B}_{n,r}$ while penalizing
$\|R\|_{2}\leq\sqrt{\alpha}.$ Now, we summarize our contributions as follows:

  * •

We study low rank factorization of graph equipartition SDP and prove the
equivalence between the factorized problem and the original convex problem in
terms of global optimality. We also prove that under the constraint
nondegeneracy condition [34], any second order stationary point of the
factorized problem is a global optimal solution with probability 1. These
results generalize the results of Boumal et al. from [4, 6, 7, 29] since we
have an extra non-smooth inequality constraint and our feasible region is not
necessarily a smooth manifold. Our study here may serve as a prototype for
extending the non-convex factorization approach to other types of SDP problems
whose resulting feasible regions may not be smooth manifolds.

  * •

We study a special algebraic variety $\mathcal{B}_{n,r}$, which is the set of
zeros of the equality constraints of the factorized problem. We study the
local geometric properties of $\mathcal{B}_{n,r}$ on both smooth points and
singular points. We also show that retraction is equivalent to the geometric
median problem under certain conditions. These properties allow us to conduct
manifold optimization methods on the smooth part of $\mathcal{B}_{n,r}$ while
escaping from non-optimal singular points of $\mathcal{B}_{n,r}$.

  * •

We develop a gradient descent method and an augmented Lagrangian method on an
algebraic variety to solve graph equipartition SDP problems and conduct
numerical experiments to demonstrate their high efficiency. Although manifold
structures are increasingly used in algorithms for solving SDP problems (see
[6, 24, 40]), our novelty here is in developing algorithms to solve the SDP
problems on an algebraic variety that is not necessarily a smooth manifold.

###  1.3 Organization of this Paper

The paper is organised as follows, in section 2, we prove the equivalence of
the upper bounded linear SDP problem (5) and the low-rank factorized model (6)
under a certain regularity condition. We also establish a sufficient condition
for a second order stationary point of (6) to be a global optimal solution. In
section 3, we study the analyse the local geometric properties of
$\mathcal{B}_{n,r}$. We show that its projection mapping and retraction
mapping around smooth points can be computed efficiently. In section 4, we
will design algorithms to solve graph equipartition SDP problems. In section
5, we conduct numerical experiments to evaluate the performance of our
proposed approach and algorithm. In section 6, we give a brief conclusion. We
show how to deal with the singular points in $\mathcal{B}_{n,r}$ in the
appendix. Some auxiliary results and proofs are also put in the appendix.

###  1.4 Notation

Note that we often omit stating the dimension of a vector or matrix if it is
already clear from the context. We use the following notation.

  * (1)

$e\in\mathbb{R}^{n}$ is the all ones vector and $J:=I-ee^{\top}/n.$

  * (2)

$e_{i}$ is a column vector such that its $i$th entry is 1 and all other
entries are zero. We will omit the dimension of $e_{i}$ if it is already
clear.

  * (3)

${\rm S}^{n-1}:=\left\\{x\in\mathbb{R}^{n}:\ \|x\|_{2}=1\right\\}$ is the unit
sphere, and ${\rm St}(n,p):=\left\\{X\in\mathbb{R}^{n\times p}:\
X^{\top}X=I\right\\}$ is the Stiefel manifold.

  * (4)

$\mathbb{S}^{n}$ and $\mathbb{S}^{n}_{+}$ denote the space of $n\times n$
symmetric matrices and its subset of positive semidefinite matrices,
respectively. For notational simplicity, we use $X\succeq 0$ to denote
$X\in\mathbb{S}_{+}^{n}$. We denote the trace of a square matrix $X$ as ${\rm
Tr}(X).$ We use $\|X\|$ or $\|X\|_{F}$ to denote the Frobenius norm of a given
matrix $X.$

  * (5)

For $x\in\mathbb{R}^{n}$, ${\rm diag}(x)\in\mathbb{R}^{n\times n}$ is the
diagonal matrix with diagonal entries $x$. For $X\in\mathbb{R}^{n\times n}$,
${\rm diag}(X)\in\mathbb{R}^{n}$ is the diagonal vector of $X$, and ${\rm
Diag}(X)={\rm diag}({\rm diag}(X)).$

  * (6)

${\rm Proj}_{F}(X)$ is the orthogonal projection of $X$ onto $F,$ which could
be a linear space or a compact set. ${\rm Proj}_{F^{\perp}}(X)=X-{\rm
Proj}_{F}(X).$

  * (7)

For $A,B\in\mathbb{R}^{m\times n}$, $A\circ B\in\mathbb{R}^{m\times n}$ is
defined by $(A\circ B)_{ij}=a_{ij}b_{ij}.$

###  1.5 Preliminaries.

In this paper, we will frequently use the definition of tangent cone (see [5]
section 2). For any closed set $S\subset\mathbb{R}^{n}$ and $x\in S,$ the
inner tangent cone $\mathcal{T}^{i}_{S}(x)$ is defined as follows:

$\mathcal{T}^{i}_{S}(x):=\left\\{h\in\mathbb{R}^{n}:\ {\rm
dist}\left(x+th,S\right)=o(t),\ t\geq 0\right\\}.$ (10)

The contingent tangent cone is defined as follows:

$\mathcal{T}_{S}(x):=\left\\{h\in\mathbb{R}^{n}:\ \exists t_{k}\downarrow 0,\
{\rm dist}\left(x+t_{k}h,S\right)=o(t_{k})\right\\}.$ (11)

It is easy to see that $\mathcal{T}^{i}_{S}(x)\subset\mathcal{T}_{S}(x).$ Note
that there is also another tangent cone called Clarke tangent cone (see [5]
section 2), but for a closed convex set $S$ that is non-singleton, these three
cones are equal to each other (see [5] proposition 2.55-2.57). Therefore, when
we deal with a closed convex set $S$, we simply use $\mathcal{T}_{S}(x)$ to
denote the tangent cone.

In order to understand more about the local geometric property of a given
closed set $S\subset\mathbb{R}^{n}$, we need second order tangent sets which
are defined as follows (see [5] definition 3.28):

$\mathcal{T}^{i,2}_{S}(x,h):=\left\\{w\in\mathbb{R}^{n}:\ {\rm
diag}\Big{(}x+th+\frac{1}{2}t^{2}w,S\Big{)}=o(t^{2}),t\geq 0\right\\},\ {\rm
for}\ h\in\mathcal{T}^{i}_{S}(x),$ (12)
$\mathcal{T}^{2}_{S}(x,h):=\left\\{w\in\mathbb{R}^{n}:\ \exists
t_{k}\downarrow 0,\ {\rm
diag}\Big{(}x+t_{k}h+\frac{1}{2}t_{k}^{2}w,S\Big{)}=o(t_{k}^{2})\right\\},\
{\rm for}\ h\in\mathcal{T}_{S}(x),$ (13)

where the first and second sets are called inner second order tangent set and
outer second order tangent set, respectively.

##  2 Relation between SDP1 and SDPLR1

###  2.1 Low-rank property of SDP1

The following theorem states the low-rank property of a linear SDP with an
additional semidefinite upper bound constraint of the form $\alpha I\succeq
X$. This result here generalizes the results in [26, 9, 10].

######  Theorem 2.1.

Consider the following SDP222Note that we use a general affine constraint
$\mathcal{B}\left(X\right)=b$ to replace the equality constraints in (5),
excluding the redundant constraint $Xe=0$.:

$\min\left\\{\left\langle C,X\right\rangle:\ \mathcal{B}(X)=b,\ X\succeq 0,\
\alpha I\succeq X\right\\},$ (14)

where $\alpha>0,$ $b\in\mathbb{R}^{m}$ has dimension $m\geq 1$, and
$\mathcal{B}:\mathbb{S}^{n}\rightarrow\mathbb{R}^{m}$ is a linear mapping such
that its adjoint $\mathcal{B}^{*}$ is injective333One can remove this
condition and replace the rank bound
$\lceil\frac{\beta}{\alpha}\rceil+\sqrt{2m}-1$ with
$\lceil\frac{\beta}{\alpha}\rceil+\sqrt{2\cdot{\rm
rank}\left(\mathcal{B}\right)}-1,$ where ${\rm rank}\left(\mathcal{B}\right)$
is the dimension of the range space of the linear mapping $\mathcal{B}.$ The
proof is similar.. If (14) has a feasible solution and

$\beta:=\max\left\\{{\rm Tr}(X):\ \mathcal{B}\left(X\right)=b,\ \alpha
I\succeq X\succeq 0\right\\},$

then (14) has a solution $X$ such that ${\rm
rank}\left(X\right)\leq\lceil\frac{\beta}{\alpha}\rceil+\sqrt{2m}-1.$

###### Proof.

Since the feasible region is non-empty and compact, an optimal solution
exists. Suppose $X$ is an optimal solution such that it has the
maximal444Since the number of eigenvalues equal to $\alpha$ is a bounded
integer, we can always choose such an $X$ out of the optimal solution set.
number of eigenvalues equal to $\alpha$, say $\ell$. Let $r={\rm
rank}\left(X\right).$ Consider the eigenvalue decomposition $X=P\Lambda
P^{\top},$ with $\Lambda={\rm diag}(\lambda)$ and the eigenvalues $\lambda$
are arranged such that
$\lambda_{1}=\lambda_{2}=\ldots=\lambda_{\ell}=\alpha>\lambda_{\ell+1}\geq\ldots\geq\lambda_{r}>0=\lambda_{r+1}=\ldots=\lambda_{n}.$

It is easy to see that we have $\ell\alpha\leq{\rm Tr}(X)\leq\beta,$ so
$\ell\leq\frac{\beta}{\alpha}.$ If $\ell=\frac{\beta}{\alpha}$, then ${\rm
rank}\left(X\right)=\ell=\frac{\beta}{\alpha}\leq\lceil\frac{\beta}{\alpha}\rceil+\sqrt{2m}-1,$
which satisfies the bound. If $\ell<\frac{\beta}{\alpha}$, then
$\ell\leq\lceil\frac{\beta}{\alpha}\rceil-1.$ Now, suppose
$r>\lceil\frac{\beta}{\alpha}\rceil+\sqrt{2m}-1$, then $s:=r-\ell>\sqrt{2m}$.
Let $\Lambda_{1}={\rm
diag}\big{(}\left(\lambda_{\ell+1},\ldots,\lambda_{r}\right)^{\top}\big{)},$
and $P_{1}$ be its corresponding matrix of eigenvectors. Consider the
following linear system:

$\mathcal{B}\big{(}P_{1}SP_{1}^{\top}\big{)}=\vec{0},\ S\in\mathbb{S}^{s}$
(15)

Because ${\rm dim}\left(\mathbb{S}^{s}\right)=\frac{s(s+1)}{2}>m$, the linear
system has a solution $S\neq{\bf 0}.$ Consider $X_{t}=X+tP_{1}SP_{1}^{\top}.$
We have $\mathcal{B}\left(X_{t}\right)=\mathcal{B}\left(X\right)=b$ for any
$t\in\mathbb{R}.$ For $t\neq 0$ sufficiently small, $X_{t}$ is still in the
feasible region of (14). From the optimality of $X$ and varying $t$ over
slightly positive and negative values, we get $\left\langle
C,P_{1}SP_{1}^{\top}\right\rangle=0$. Thus $\left\langle
C,X_{t}\right\rangle=\left\langle C,X\right\rangle$ for any $t\in\mathbb{R}$,
and $X_{t}$ is an optimal solution for all $t$ sufficiently small. Let
$\mathcal{T}:=\left\\{t\in\mathbb{R}:\ X_{t}\succeq 0,\ \alpha I\succeq
X_{t}\right\\}.$ Since $S\neq 0$, $\mathcal{T}$ is a bounded closed interval
$[t_{1},t_{2}]$. We consider $X_{t_{1}}.$ It has eigenvalues
$\lambda_{1}^{\prime}=\lambda_{2}^{\prime}=\ldots=\lambda_{\ell}^{\prime}=\alpha\geq\lambda_{\ell+1}^{\prime}\geq\ldots\geq\lambda_{r}^{\prime}\geq
0=\lambda_{r+1}^{\prime}=\ldots=\lambda_{n}^{\prime}=0.$ Since $t_{1}$ is on
the boundary of $\mathcal{T}$, either $\lambda_{\ell+1}^{\prime}=\alpha$ or
$\lambda_{r}^{\prime}=0.$ However, since $X$ has the maximum number of
eigenvalues equal to $\alpha$, $\lambda_{\ell+1}^{\prime}<\alpha$ and we have
$\lambda_{r}^{\prime}=0.$ Therefore, we have found another optimal solution
$X_{t_{1}}$ whose number of eigenvalues equal to $\alpha$ is still maximal but
has rank less than ${\rm rank}\left(X\right).$ We can continue this operation
from $X_{t_{1}}$ until we find an optimal solution with rank
$\leq\lceil\frac{\beta}{\alpha}\rceil+\sqrt{2m}-1.$

If we do not have the constraint $\alpha I\succeq X$ in Theorem 2.1 but
$\left\\{X\in\mathbb{S}^{n}:\ X\succeq 0,\ \mathcal{B}(X)=b\right\\}$ is
compact, then we can choose a sufficiently large $\alpha$ such that any point
$X$ in the compact feasible region satisfies ${\rm Tr}(X)\leq\alpha$ (which
also implies that $X\preceq\alpha I$). In this case, the constant $\beta$ in
Theorem 2.1 is at most $\alpha$ and then the rank-bound is $\sqrt{2m},$ which
is almost the same as the classical bound of
$\frac{1}{2}\left(\sqrt{8m+1}-1\right).$

Let us now consider (5), we can remove $Xe=0$ since it is redundant. Then
there are $m$ affine constraints. From Theorem 2.1, there exists an optimal
solution of rank $\leq\lceil\frac{\beta}{\alpha}\rceil+\sqrt{2m}-1.$ This
means that (5) is equivalent to (6) as long as
$r\geq\lceil\frac{\beta}{\alpha}\rceil+\sqrt{2m}-1.$ For (4), since ${\rm
diag}(X)=e$, we may choose $\beta=n$ and the rank bound is
$k+\sqrt{2(n+1)}-2.$

###  2.2 A sufficient condition for a second order stationary point of SDPLR1
to be globally optimal

Here we establish a sufficient condition for a second order stationary point
of the non-convex problem SDPLR1 to be globally optimal. It turns out that the
condition is closely related to the constraint nondegeneracy condition of the
problem SDPLR1 (6). Thus in the first part of this subsection, we elaborate on
the notation for matrix decomposition and several tangent cones and normal
cones. In the second part, we state the regularity condition and its
equivalence to the constraint nondegeneracy condition of (6). In the third
part, we prove that under the regularity condition, a rank-deficient second
order stationary point of (6) is also a global optimal solution. In the last
part, we show that the rank-deficient condition usually holds provided the
parameter $r$ in (6) is large enough.

Let $R\in\mathbb{R}^{n\times r}$ be a feasible solution of (6). Then
$RR^{\top}$ is a feasible solution of (5). From now on, we suppose $R$ has the
following singular value decomposition (SVD)555These notation for the SVD will
be used frequently later.

$\displaystyle
R=\begin{pmatrix}U_{1}&U_{2}&U_{3}\end{pmatrix}\begin{pmatrix}\sqrt{\Lambda_{1}}&&\\\
&\sqrt{\Lambda_{2}}&\\\ &&{\bf 0}\end{pmatrix}\ \begin{pmatrix}V_{1}^{\top}\\\
V_{2}^{\top}\\\ V_{3}^{\top}\end{pmatrix},$ (16)

where $U_{1}\in\mathbb{R}^{n\times r_{1}}$, $U_{2}\in\mathbb{R}^{n\times
r_{2}}$, $U_{3}\in\mathbb{R}^{n\times(n-r_{1}-r_{2})}$ and
$\begin{pmatrix}U_{1}&U_{2}&U_{3}\end{pmatrix}\in{\rm St}(n,n)$;
$V_{1}\in\mathbb{R}^{r\times r_{1}},$ $V_{2}\in\mathbb{R}^{r\times r_{2}},$
$V_{3}\in\mathbb{R}^{r\times(r-r_{1}-r_{2})}$ and
$\begin{pmatrix}V_{1}&V_{2}&V_{3}\end{pmatrix}\in{\rm St}(r,r)$;
$\sqrt{\Lambda_{1}}$ corresponds to singular values equal to $\sqrt{\alpha},$
$\sqrt{\Lambda_{2}}$ corresponds to singular values in the interval
$(0,\sqrt{\alpha})$; ${\bf 0}$ is the zero matrix in
$\mathbb{R}^{(n-r_{1}-r_{2})\times(r-r_{1}-r_{2})}.$ From $R^{\top}e=0$ and
hence $U_{1}^{\top}e=0$, $U_{2}^{\top}e=0$, we have

$J\left(U_{1},U_{2},R\right)=\left(U_{1},U_{2},R\right).$

Based on the SVD of $R$, we know that $X=RR^{\top}$ has the following
eigenvalue decomposition,

$X=\begin{pmatrix}U_{1}&U_{2}&U_{3}\end{pmatrix}\begin{pmatrix}\Lambda_{1}&&\\\
&\Lambda_{2}&\\\ &&\hat{\bf 0}\end{pmatrix}\ \begin{pmatrix}U_{1}^{\top}\\\
U_{2}^{\top}\\\ U_{3}^{\top}\end{pmatrix},$

where $\hat{\bf 0}\in\mathbb{R}^{(n-r_{1}-r_{2})\times(n-r_{1}-r_{2})}.$ From
these decompositions, we have ${\rm rank}(R)={\rm rank}(X)=r_{1}+r_{2}.$

We let $\mathcal{T}_{\Omega_{\alpha}^{1}}(X)$,
$\mathcal{T}_{\Omega_{\alpha}^{2}}(R)$,
$\mathcal{N}_{\Omega_{\alpha}^{1}}(X)$, $\mathcal{N}_{\Omega_{\alpha}^{2}}(R)$
be their tangent cones and normal cones666Note that the standard notation for
normal cone is
$\mathcal{N}_{\mathcal{T}_{\Omega_{\alpha}}\left(X\right)}\left(Y\right).$ As
we always consider $Y=0,$ so we choose a more convenient notation. . From
proposition 2.5 and proposition 2.9 in [13] (also see [43]), if
$\lambda_{\max}\left(X\right)=\alpha$ and $\lambda_{\min}\left(X\right)=0,$
then for any $H\in\mathbb{S}^{n}$, $Y\in\mathbb{R}^{n\times r},$ we have the
following formulas for various directional derivatives:

$\lambda_{\max}^{\prime}\left(X;H\right)=\lambda_{\max}\big{(}U_{1}^{\top}HU_{1}\big{)},\
\lambda_{\min}^{\prime}\left(X;H\right)=\lambda_{\min}\big{(}U_{3}^{\top}HU_{3}\big{)}$

and

$\sigma_{\max}^{\prime}(R;Y)=\lambda_{\max}\left(\frac{U_{1}^{\top}YV_{1}+V_{1}^{\top}Y^{\top}U_{1}}{2}\right).$

Thus, we have777Those cones are also correct if
$\lambda_{\max}\left(X\right)<\alpha$ or $\lambda_{\min}\left(X\right)>0.$

$\displaystyle\mathcal{T}_{\Omega_{\alpha}^{1}}\left(X\right)$
$\displaystyle=$ $\displaystyle\left\\{H\in\mathbb{S}^{n}:\
U_{1}^{\top}HU_{1}\preceq 0,\ U_{3}^{\top}HU_{3}\succeq 0\right\\},$
$\displaystyle\mathcal{T}_{\Omega_{\alpha}^{2}}\left(R\right)$
$\displaystyle=$ $\displaystyle\left\\{Y\in\mathbb{R}^{n\times r}:\
U_{1}^{\top}YV_{1}+V_{1}^{\top}Y^{\top}U_{1}\preceq 0\right\\},$
$\displaystyle\mathcal{N}_{\Omega_{\alpha}^{1}}\left(X\right)$
$\displaystyle=$ $\displaystyle\left\\{U_{1}SU_{1}^{\top}-U_{3}WU_{3}^{\top}:\
S\in\mathbb{S}^{r_{1}}_{+},\ W\in\mathbb{S}^{n-r_{1}-r_{2}}_{+}\right\\},$
$\displaystyle\mathcal{N}_{\Omega_{\alpha}^{2}}\left(R\right)$
$\displaystyle=$ $\displaystyle\left\\{U_{1}ZV_{1}^{\top}:\
Z\in\mathbb{S}^{r_{1}}_{+}\right\\}.$

Since (5) is a convex problem, the following condition is a sufficient
optimality condition for SDP1.

######  Proposition 2.2.

If $\hat{X}$ is feasible for (5) and there exists $y\in\mathbb{R}^{m-1}$ such
that
$J\left(C-\mathcal{A}^{*}(y)\right)J\in-\mathcal{N}_{\Omega_{\alpha}^{1}}(\hat{X}),$
then $\hat{X}$ is an optimal solution of (5).

###### Proof.

Define
$\tilde{\mathcal{A}}:\mathbb{S}^{n}\rightarrow\mathbb{R}^{m-1}\times\mathbb{R}^{n}\times\mathbb{R}$
such that
$\tilde{\mathcal{A}}\left(X\right)=\left(\mathcal{A}(X),Xe,\left\langle
X,ee^{\top}\right\rangle\right).$ Then
$\tilde{\mathcal{A}}^{*}:\mathbb{R}^{m-1}\times\mathbb{R}^{n}\times\mathbb{R}\rightarrow\mathbb{S}^{n}$
is given by
$\tilde{\mathcal{A}}^{*}(y,z,\beta)=\mathcal{A}^{*}y+\frac{1}{2}(ez^{\top}+ze^{\top})+\beta
ee^{\top}.$ Define $\tilde{b}:=\left(b,\vec{0},0\right).$ Then (5) can be
simplified as:

$\min\left\\{\left\langle C,X\right\rangle:\
\tilde{\mathcal{A}}(X)=\tilde{b},\ X\in\Omega_{\alpha}^{1}\right\\}.$

Note that
$J\left(C-\mathcal{A}^{*}(y)\right)J\in-\mathcal{N}_{\Omega^{1}_{\alpha}}(\hat{X})$
implies that there exists
$\tilde{y}\in\mathbb{R}^{m-1}\times\mathbb{R}^{n}\times\mathbb{R}$ such that
$C-\tilde{\mathcal{A}}^{*}(\tilde{y})\in-\mathcal{N}_{\Omega_{\alpha}^{1}}(\hat{X}).$
Now, suppose $\hat{X}$ is not an optimal solution, then there exists
$\bar{X}\in\mathbb{S}^{n}$ such that $\bar{X}\in\Omega_{\alpha}^{1},$
$\tilde{\mathcal{A}}(\bar{X})=\tilde{b}$ and $\left\langle
C,\bar{X}\right\rangle<\big{\langle}C,\,\hat{X}\big{\rangle}.$ Then we have
$\tilde{\mathcal{A}}(\hat{X}-\bar{X})=0$ which implies that
$\big{\langle}\tilde{\mathcal{A}}^{*}(\tilde{y}),\,\hat{X}-\bar{X}\big{\rangle}=0.$
Hence
$\big{\langle}C-\tilde{\mathcal{A}}^{*}(\tilde{y}),\,\hat{X}-\bar{X}\big{\rangle}>0.$
Since
$C-\tilde{\mathcal{A}}^{*}(\tilde{y})\in-\mathcal{N}_{\Omega_{\alpha}^{1}}(\hat{X})$
and $\bar{X}\in{\Omega_{\alpha}^{1}}(\hat{X})$, we get by the definition of
the normal cone that
$\big{\langle}C-\tilde{\mathcal{A}}^{*}(\tilde{y}),\,\hat{X}-\bar{X}\big{\rangle}\leq
0,$ which is a contradiction. Thus $\hat{X}$ must be an optimal solution of
(5).

Now, we define the mapping $G:\ \mathbb{R}^{n\times
r}\rightarrow\mathbb{R}^{m-1}\times\mathbb{R}^{r}\times\mathbb{R}^{n\times r}\
{\rm such\ that}$

$G(R):=\Big{(}\frac{\mathcal{A}\left(RR^{\top}\right)-b}{2},R^{\top}e,R\Big{)}.$

Let $\mathcal{K}:=\\{0\\}^{m-1}\times\\{0\\}^{r}\times\Omega_{\alpha}^{2}.$
Then the problem SDPLR1 (6) can be simplified as

$\min\left\\{\frac{1}{2}\big{\langle}C,\,RR^{\top}\big{\rangle}:\
G(R)\in\mathcal{K}\right\\}.{\rm}$ (6 ′)

Robinson’s constraints qualification (CQ) [32] is said to hold at $R$ for (6
′) if

$0\in{\rm int}\left\\{G(R)+G^{\prime}(R)\mathbb{R}^{n\times
r}-\mathcal{K}\right\\}.$ (17)

Note that by Corollary 2.98 in [5], the above condition is equivalent to

$\displaystyle G^{\prime}(R)\mathbb{R}^{n\times
r}-\mathcal{T}_{\mathcal{K}}(G(R))=\mathbb{R}^{m-1}\times\mathbb{R}^{r}\times\mathbb{R}^{n\times
r}$ (18)

where $\mathcal{T}_{\mathcal{K}}(G(R))$ is the tangent cone to $\mathcal{K}$
at $G(R).$

For any $R\in\mathbb{R}^{n\times r}$, define the linear mapping
$g_{\mbox{\tiny$R$}}:\mathbb{R}^{n\times
r}\rightarrow\mathbb{R}^{m-1}\times\mathbb{R}^{r}$ by

$g_{\mbox{\tiny$R$}}(H)=\Big{(}\frac{1}{2}{\mathcal{A}\big{(}RH^{\top}+HR^{\top}\big{)}},H^{\top}e\Big{)}.$

Then $g_{\mbox{\tiny$R$}}^{*}(y,z)=\mathcal{A}^{*}(y)R+ez^{\top}.$ Note that
$g_{\mbox{\tiny$R$}}$ is the differential of the equality constraints of (6)
and $G^{\prime}(R)[H]=\left(g_{\mbox{\tiny$R$}}(H),H\right)$. Let

$\Gamma_{R}:=\left\\{H\in\mathbb{R}^{n\times r}:\
U_{1}^{\top}HV_{1}+V_{1}^{\top}H^{\top}U_{1}=0\right\\}.$

The following definition is important for our later analysis.

######  Definition 2.3.

$R\in\mathbb{R}^{n\times r}$ such that $G(R)\in\mathcal{K}$ is called a
regular point for (6 ′) if
$g_{\mbox{\tiny$R$}}\left(\Gamma_{R}\right)=\mathbb{R}^{m-1}\times\mathbb{R}^{r}.$

The above definition of a regular point is equivalent to the constraint
nondegeneracy condition in Definition 2.1 of [34], which states that

$\displaystyle G^{\prime}(R){\cal X}+{\rm
lin}(\mathcal{T}_{\mathcal{K}}(G(R))={\cal Y}$ (19)

where ${\cal X}=\mathbb{R}^{n\times r}$, ${\cal
Y}=\mathbb{R}^{m-1}\times\mathbb{R}^{r}\times\mathbb{R}^{n\times r}$,
$\mathcal{K}=\\{0\\}^{m-1}\times\\{0\\}^{r}\times\Omega_{\alpha}^{2}$, and
${\rm lin}(\mathcal{T}_{\mathcal{K}}(G(R))$ is the largest linear subspace
contained in $\mathcal{T}_{\mathcal{K}}(G(R)).$ In particular

${\rm
lin}(\mathcal{T}_{\mathcal{K}}(G(R))=\\{0\\}^{m-1}\times\\{0\\}^{r}\times{\rm
lin}(\mathcal{T}_{\Omega^{2}_{\alpha}}(R)).$

It is easy to show that the constraint nondegeneracy condition is equivalent
to

$g_{\mbox{\tiny$R$}}({\rm
lin}(\mathcal{T}_{\Omega^{2}_{\alpha}}(R)))=\mathbb{R}^{m-1}\times\mathbb{R}^{r}.$

Since

${\rm lin}(\mathcal{T}_{\Omega^{2}_{\alpha}}(R))=\\{Y\in\mathbb{R}^{n\times
r}\mid U_{1}^{\top}YV_{1}+V_{1}^{\top}Y^{\top}U_{1}=0\\}=\Gamma_{R},$

we have shown the equivalence between a regular point and the constraint
nondegeneracy condition.

######  Remark 2.4.

Constraint nondegeneracy condition is also defined in [5] Remark 4.72. They
additionally require that the set $\mathcal{K}$ is $\mathcal{C}^{1}$-cone
reducible (see Definition 3.135 of [5]). One can prove that
$\mathcal{K}=\\{0\\}^{m-1}\times\\{0\\}^{r}\times\Omega_{\alpha}^{2}$ is
$\mathcal{C}^{\infty}$-cone reducible using a similar approach as in Example
3.140 of [5]. We omit the proof here since we use the more convenient
definition from [34].

From (18) and (19), we have the following lemma.

######  Lemma 2.5.

If $R\in\mathbb{R}^{n\times r}$ such that $G(R)\in\mathcal{K}$ is regular,
then it satisfies Robinson’s constraints qualification for (6 ′).

The following lemma gives a sufficient condition for a second order stationary
point of (6) to be a global optimal solution.

######  Lemma 2.6.

Suppose $R\in\mathbb{R}^{n\times r}$ such that $G(R)\in\mathcal{K}$ is a
regular rank deficient second order stationary point of (6). Then $RR^{\top}$
is a global optimal solution of (5) and so $R$ is a global optimal solution of
(6).

###### Proof.

First, from definition of $G$, we have

$G^{\prime}(R)[H]=\left(g_{\mbox{\tiny$R$}}(H),H\right)=\Big{(}\mathcal{A}\Big{(}\frac{RH^{\top}+HR^{\top}}{2}\Big{)},H^{\top}e,H\Big{)},\
G^{\prime\prime}(R)[H,H]=\big{(}\mathcal{A}(HH^{\top}),0,0\big{)}.$

Since $R$ is a regular second order stationary point of (6), by Lemma 2.5, it
satisfies Robinson’s constraints qualification. This means that it satisfies
first and second order necessary KKT conditions (See [5] (3.16) and (3.99)),
i.e. there exists $(y,z)\in\mathbb{R}^{m-1}\times\mathbb{R}^{r}$ such that

$G(R)\in\mathcal{K},\
\left(C-\mathcal{A}^{*}(y)\right)R-ez^{\top}\in-\mathcal{N}_{\Omega_{\alpha}^{2}}(R);$
(20) $\forall\ d,w\in\mathbb{R}^{n\times r}\ \mbox{s.t}\
G^{\prime}(R)d=\left(g_{\mbox{\tiny$R$}}(d),d\right)\in\mathcal{T}_{\mathcal{K}}\left(G(R)\right)\
{\rm and}$
$G^{\prime}(R)w+G^{\prime\prime}(R)\left(d,d\right)\in\mathcal{T}^{i,2}_{\mathcal{K}}\left(G(R),G^{\prime}(R)d\right),$
${\rm if}\ \left\langle CR,d\right\rangle=0,\ {\rm then}\ \left\langle
CR,w\right\rangle+\left\langle Cd,d\right\rangle\geq 0.$ (21)

Here $\mathcal{T}_{\mathcal{K}}^{i,2}(G(R),G^{\prime}(R)d)$ is the inner
second order tangent set defined in (12) as follows:

$\mathcal{T}_{\mathcal{K}}^{i,2}(G(R),G^{\prime}(R)d)=\left\\{v\in\mathbb{R}^{n\times
r}:\ {\rm
dist}\Big{(}G(R)+tG^{\prime}(R)d+\frac{t^{2}}{2}v,\mathcal{K}\Big{)}=o(t^{2}),\
t\geq 0\right\\}.$

From (20), we know that $X=RR^{\top}$ is feasible for (5) and there exists
$Z\in\mathbb{S}^{r_{1}}$ such that $Z\preceq 0$ and

$\left(C-\mathcal{A}^{*}(y)\right)R-ez^{\top}=U_{1}ZV_{1}^{\top}.$ (22)

Our goal next is to show that

$J\left(C-\mathcal{A}^{*}y\right)J\in-\mathcal{N}_{\Omega_{\alpha}^{1}}\left(X\right).$

Multiplying $J$ to the left-hand side of (22) and noting that
$J(U_{1},U_{2},R)=(U_{1},U_{2},R),$ we get

$J\left(C-\mathcal{A}^{*}(y)\right)JR=U_{1}ZV_{1}^{\top}.$ (23)

Thus we have

$\displaystyle\begin{pmatrix}U_{1}^{\top}\\\ U_{2}^{\top}\\\
U_{3}^{\top}\end{pmatrix}J\left(C-\mathcal{A}^{*}(y)\right)J\begin{pmatrix}U_{1}&U_{2}&U_{3}\end{pmatrix}\begin{pmatrix}\sqrt{\Lambda_{1}}&&\\\
&\sqrt{\Lambda_{2}}&\\\ &&{\bf 0}\end{pmatrix}=\begin{pmatrix}Z&0&0\\\
0&0&0\\\ 0&0&0\end{pmatrix},$

which implies that

$J\left(C-\mathcal{A}^{*}(y)\right)J=\begin{pmatrix}U_{1}&U_{2}&U_{3}\end{pmatrix}\begin{pmatrix}Z\sqrt{\alpha}^{-1}&0&0\\\
0&0&0\\\ 0&0&W\end{pmatrix}\begin{pmatrix}U_{1}^{\top}\\\ U_{2}^{\top}\\\
U_{3}^{\top}\end{pmatrix},$ (24)

for some $W\in\mathbb{S}^{n-r_{1}-r_{2}}.$ For any
$S\in\mathbb{R}^{(n-r_{1}-r_{2})\times(r-r_{1}-r_{2})}$, define
$d_{S}=JU_{3}SV_{3}^{\top}\in\mathbb{R}^{n\times r}.$ Note that we can define
$d_{S}$ because $R$ is rank deficient and so $r-r_{1}-r_{2}>0$. Then
$d_{S}^{\top}e=0$, and $Rd_{S}^{\top}=RV_{3}S^{\top}U_{3}^{\top}J=0$ because
$RV_{3}=0.$ So $g_{\mbox{\tiny$R$}}(d_{S})=0.$ We also have

$R+td_{S}=J\begin{pmatrix}U_{1}&U_{2}&U_{3}\end{pmatrix}\begin{pmatrix}\sqrt{\Lambda_{1}}&&\\\
&\sqrt{\Lambda_{2}}&\\\ &&tS\end{pmatrix}\begin{pmatrix}V_{1}^{\top}\\\
V_{2}^{\top}\\\ V_{3}^{\top}\end{pmatrix},$

which means that for $t$ sufficiently small
$\|R+td_{S}\|_{2}\leq\sqrt{\alpha}$ i.e. $R+td_{S}\in\Omega_{\alpha}^{2}.$
Therefore,

$G^{\prime}(R)d_{S}\in\mathcal{T}_{\mathcal{K}}\left(G(R)\right).$

From
$g_{\mbox{\tiny$R$}}\left(\Gamma_{R}\right)=\mathbb{R}^{m-1}\times\mathbb{R}^{r},$
we may choose $w_{S}\in\Gamma_{R}$ such that

$g_{\mbox{\tiny$R$}}\left(w_{S}\right)=\left(-\mathcal{A}\big{(}d_{S}d_{S}^{\top}\big{)},0\right).$

Thus
$\mathcal{A}\left(\frac{Rw_{S}^{\top}+w_{S}R^{\top}}{2}+d_{S}d_{S}^{\top}\right)=0,$
and $w_{S}^{\top}e=0$.

Next we show that

$\|R+td_{S}+\mbox{$\frac{1}{2}$}t^{2}w_{S}\|_{2}\leq\sqrt{\alpha}+o(t^{2}).$
(25)

For the case where $\|R\|_{2}<\sqrt{\alpha},$ it is easy to see that (25)
holds true. Now suppose $\|R\|_{2}=\sqrt{\alpha}.$ Then we have

$\displaystyle\lambda_{\max}\Big{(}\Big{(}R+td_{S}+\frac{t^{2}}{2}w_{S}\Big{)}\Big{(}R+td_{S}+\frac{t^{2}}{2}w_{S}\Big{)}^{\top}\Big{)}$
$\displaystyle=\lambda_{\max}\Big{(}RR^{\top}+tRd_{S}^{\top}+td_{S}R^{\top}+t^{2}\Big{(}d_{S}d_{S}^{\top}+\frac{w_{S}R^{\top}+Rw_{S}^{\top}}{2}\Big{)}\Big{)}+O(t^{3})$
$\displaystyle=\lambda_{\max}\Big{(}RR^{\top}+t^{2}\Big{(}d_{S}d_{S}^{\top}+\frac{w_{S}R^{\top}+Rw_{S}^{\top}}{2}\Big{)}\Big{)}+O(t^{3})$
$\displaystyle=\lambda_{\max}\big{(}RR^{\top}\big{)}+t^{2}\lambda_{\max}\Big{(}U_{1}^{\top}\Big{(}d_{S}d_{S}^{\top}+\frac{w_{S}R^{\top}+Rw_{S}^{\top}}{2}\Big{)}U_{1}\Big{)}+o(t^{2})=\alpha+o(t^{2}),$

where the first equality comes from the Lipschitz continuity of
$\lambda_{\max},$ the second equality comes from $Rd_{S}^{\top}=0$, the third
equality comes from the directional derivative of $\lambda_{\max},$ the fourth
equality comes from $U_{1}^{\top}d_{S}=0$ and
$U_{1}^{\top}w_{S}R^{\top}U_{1}+U_{1}^{\top}Rw_{S}^{\top}U_{1}=0,$ where the
latter holds because $w_{S}\in\Gamma_{R}.$ This completes the proof of (25).

Now for $t>0$, we have

$\displaystyle{\rm
dist}\left(R+td_{S}+\frac{1}{2}t^{2}w_{S},\Omega_{\alpha}^{2}\right)=o(t^{2})$
$\displaystyle\left(R+td_{S}+\frac{1}{2}t^{2}w_{S}\right)^{\top}e=0$
$\displaystyle\mathcal{A}\left(\Big{(}R+td_{S}+\frac{1}{2}t^{2}w_{S}\Big{)}\Big{(}R+td_{S}+\frac{1}{2}t^{2}w_{S}\Big{)}^{\top}\right)-b$
$\displaystyle=t^{2}\mathcal{A}\left(d_{S}d_{S}^{\top}+\left(Rw_{S}^{\top}+w_{S}R^{\top}\right)/2\right)+O(t^{3})=O(t^{3}),$

where the first inequality comes from (25) and the rest come from
$g_{\mbox{\tiny$R$}}(w_{S})=\left(-\mathcal{A}\left(d_{S}d_{S}^{\top}\right),0\right).$
The above results together imply that

$G^{\prime}(R)w_{S}+G^{\prime\prime}(R)(d_{S},d_{S})\in\mathcal{T}^{i,2}_{\mathcal{K}}\left(G(R),G^{\prime}(R)d_{S}\right).$

Since $\left\langle CR,d_{S}\right\rangle=\left\langle
C,Rd_{S}^{\top}\right\rangle=0$, from (21), we have

$\left\langle CR,w_{S}\right\rangle+\left\langle Cd_{S},d_{S}\right\rangle\geq
0.$

From
$\mathcal{A}\left(d_{S}d_{S}^{\top}+\left(Rw_{S}^{\top}+w_{S}R^{\top}\right)/2\right)=0$,
we have
$\left\langle\mathcal{A}^{*}(y),d_{S}d_{S}^{\top}\right\rangle+\left\langle\mathcal{A}^{*}(y)R,w_{S}^{\top}\right\rangle=0.$
Thus we have

$\left\langle\left(C-\mathcal{A}^{*}(y)\right)R,w_{S}\right\rangle+\big{\langle}\left(C-\mathcal{A}^{*}(y)\right),\,d_{S}d_{S}^{\top}\big{\rangle}\geq
0.$

Since $Jw_{S}=w_{S}$, $JR=R$, $Jd_{S}=d_{S}$, it follows from the above
inequality that

$\left\langle
J\left(C-\mathcal{A}^{*}(y)\right)JR,w_{S}\right\rangle+\big{\langle}J\left(C-\mathcal{A}^{*}(y)\right)J,\,d_{S}d_{S}^{\top}\big{\rangle}\geq
0.$

From (23) and $w_{S}\in\Gamma_{R}$, we have $\left\langle
J\left(C-\mathcal{A}^{*}(y)\right)JR,w_{S}\right\rangle=\left\langle
U_{1}ZV_{1}^{\top},w_{S}\right\rangle=\big{\langle}Z,\,U_{1}^{\top}w_{S}V_{1}\big{\rangle}=0.$
So

$\left\langle
J\left(C-\mathcal{A}^{*}(y)\right)J,d_{S}d_{S}^{\top}\right\rangle=\big{\langle}J\left(C-\mathcal{A}^{*}(y)\right)J,\,U_{3}SS^{\top}U_{3}^{\top}\big{\rangle}\geq
0,$

for any $S\in\mathbb{R}^{(n-r_{1}-r_{2})\times(r-r_{1}-r_{2})}.$ This means
that $W=U_{3}^{\top}J\left(C-\mathcal{A}^{*}(y)\right)JU_{3}\succeq 0.$ Thus,
from (24),

$J\left(C-\mathcal{A}^{*}(y)\right)J\in-\mathcal{N}_{\Omega_{\alpha}^{1}}\left(X\right).$

From Proposition 2.2, $X=RR^{\top}$ is an optimal solution for (5).

The following proposition is inspired by Lemma 3 of [4] and Theorem 2 in [6].
It tells that the rank deficient condition usually holds.

######  Lemma 2.7.

Suppose $\beta:=\max\left\\{{\rm Tr(X)}:\ \mathcal{A}\left(X\right)=b,\ Xe=0,\
X\succeq 0\right\\}<\infty$ and
$r\geq\lceil\frac{\beta}{\alpha}\rceil+\sqrt{2m}.$ Then for all
$C\in\mathbb{S}^{n}$ except for a zero measure set, any first order stationary
point $R$ of (6) that satisfies Robinson’s constraints qualification is rank
deficient.

###### Proof.

Let $C\in\mathbb{S}^{n}$ be such that (6) has a full rank first order
stationary point $R\in\mathbb{R}^{n\times r}$ satisfying Robinson’s constraint
qualification. It necessarily satisfies the first order KKT condition (20),
which means that

$\left(C-\mathcal{A}^{*}(y)\right)R\in ez^{\top}-\left\\{U_{1}ZV_{1}^{\top}:\
Z\succeq 0\right\\}\subset\mathcal{M}_{\leq 1+r_{1}},$

where $\mathcal{M}_{\leq 1+r_{1}}=\left\\{A\in\mathbb{S}^{n}:\ {\rm
rank}\left(A\right)\leq 1+r_{1}\right\\}.$ Thus we have

${\rm rank}\left(R\right)+{\rm
rank}\left(C-\mathcal{A}^{*}(y)\right)-n\leq{\rm
rank}\left(\left(C-\mathcal{A}^{*}(y)\right)R\right)\leq r_{1}+1,$

where the first inequality follows from Sylvester’s rank inequality (see 10.5
of [28]). Note that $r_{1}\leq\lceil\frac{\beta}{\alpha}\rceil-1$ (for
otherwise because $r>\lceil\frac{\beta}{\alpha}\rceil$ and $R$ is full rank,
we would get the contradiction that ${\rm
Tr}\left(RR^{\top}\right)>\lceil\frac{\beta}{\alpha}\rceil\alpha\geq\beta$).

Since $R$ has full rank, we have

${\rm rank}\left(C-\mathcal{A}^{*}(y)\right)\leq n+r_{1}+1-r\leq n-\sqrt{2m}.$

Then $C\in\mathcal{M}_{\leq
n-\sqrt{2m}}+\mathcal{A}^{*}\left(\mathbb{R}^{m-1}\right)$ whose dimension is
bounded by
$\frac{n(n+1)}{2}-\frac{\sqrt{2m}\left(\sqrt{2m}+1\right)}{2}+m-1<\frac{n(n+1)}{2}.$
This means that $C\in\mathbb{S}^{n}$ is contained in a set of measure zero.

The following theorem directly follows from Lemma 2.6 and Lemma 2.7.

######  Theorem 2.8.

Suppose $\beta:=\max\left\\{{\rm Tr}\left(X\right):\
\mathcal{A}\left(X\right)=b,\ Xe=0,\ X\succeq 0\right\\}<\infty$ and
$r\geq\lceil\frac{\beta}{\alpha}\rceil+\sqrt{2m}.$ Then for all
$C\in\mathbb{S}^{n}$ except for a zero measure set, any regular second order
stationary point $R$ of (6) is a global minimizer of $(\ref{SDPLR}).$
Moreover, if we can get $\Theta\in\mathcal{N}_{\Omega_{2}^{\alpha}}(R)$ in the
KKT system (20), we have $g_{\mbox{\tiny$R$}}^{*}(y,z)=CR+\Theta.$ Because
$g_{\mbox{\tiny$R$}}$ is an onto linear mapping,
$(y,z)=\left(g_{\mbox{\tiny$R$}}\,g_{\mbox{\tiny$R$}}^{*}\right)^{-1}g_{\mbox{\tiny$R$}}\left(CR+\Theta\right).$
Thus, we can recover the dual variable for problem (5).

######  Remark 2.9.

We note that the construction of the dual variables $(y,Z)$ from the primal
variable $R$ in Theorem 2.8 is necessary since in our Riemannian based
algorithm, the dual variables are not constructed explicitly. We also note
that Theorem 2.8 has the prerequisite that the output of an optimization
algorithm is a regular point. Such a regularity assumption is commonly used in
the optimization literature because without it, the solution obtained by an
algorithm may not even be a KKT solution.

###  2.3 Comparison between Theorem 2.8 and Theorem 4 in [11]

We should mention that the SDP upper bound of (5) can also be handled by
introducing a slack variable $Y$ and writing problem (5) equivalently as
follows:

$\min\left\\{\left\langle C,X\right\rangle:\ \mathcal{A}(X)=b,\ \left\langle
X,ee^{\top}\right\rangle=0,\ X+Y=\alpha I,\
X,Y\in\mathbb{S}^{n}_{+}\right\\}.$ (26)

In this case, problem (26) is a multi-block SDP problem. By using the partial
BM factorization in section 4 of [11], we get the following problem:

$\min_{R\in\mathbb{R}^{n\times r}}\left\\{\left\langle
C,RR^{\top}\right\rangle:\ \mathcal{A}(RR^{\top})=b,\ \left\langle
RR^{\top},ee^{\top}\right\rangle=0,\ RR^{\top}+Y=\alpha I,\
Y\in\mathbb{S}^{n}_{+}\right\\}.$ (27)

Theorem 4 in [11] says that if

$r(r+1)/2>m+n(n+1)/2-r(Y)\left(r(Y)+1\right)/2,$ (28)

then for a generic $C\in\mathbb{S}^{n}$, any 2-critical point
$\left(R,Y\right)$ of (27) satisfies that $\left(RR^{\top},Y\right)$ is global
optimal for (26) (see section 2 of [11] for the definition of a 2-critical
point). Here $r(Y)$ is a bound on the smallest rank of $Y$ for some feasible
solution $\left(X,Y\right)$ of (26). From the definition of $\beta$ and
$\alpha I\succeq X,$ we have that $r(Y)\geq n-\beta/\alpha.$ This implies the
following rank condition of $r$:

$r(r+1)/2>m+n(n+1)/2-\left(n-\beta/\alpha\right)\left(n-\beta/\alpha+1\right)/2.$
(29)

We assume that $\beta/\alpha\geq 1$; for otherwise the constraint $\alpha
I\succeq X$ can be removed. If we have $1\leq\beta/\alpha=o(n)$999This is true
for SDP relaxations of graph equipartition and $k$-means clustering in (see
[27] (20)). Actually, for an SDP problem where the trace of $X$ is fixed as a
constant $\beta$, the SDP upper bound $\alpha I\succeq X$ implies that ${\rm
rank}\left(X\right)\geq\beta/\alpha.$ Thus, we have to assume
$\beta/\alpha=o(n)$ if we want the BM approach to be efficient. and
$m=o(n^{2})$101010This is often the case when we apply the BM approach.
Otherwise there are too many constraints, and the SDP relaxation is not
guaranteed to be low rank anymore., then (29) gives the following rank bound

$r>\sqrt{2m+(2+o(1))n\beta/\alpha}.$ (30)

In comparison, our rank condition is

$r\geq\lceil\frac{\beta}{\alpha}\rceil+\sqrt{2m}=\sqrt{2m+o(1)n\beta/\alpha},$
(31)

which is strictly better than the condition in (30), especially when $n$ is
sufficiently large. In particular, in graph equipartition with $k$ parts, (30)
implies that $r>\sqrt{(1+o(1))2kn},$ while (31) implies that
$r\geq(1+o(1))\sqrt{2n}.$ The latter is clearly better than the former bound
when $k$ is large. Moreover, if one considers the SDP relaxation of a
$k$-means clustering problem (see [27] (20)), then (30) implies that
$r>\sqrt{(1+o(1))(k-1)n},$ while our condition gives $r\geq\sqrt{2}+k-1,$
which is a constant that is independent of $n$. In contrast, the former lower
bound derived from (30) can be large when $n$ is much larger than $k$, which
is usually the case for a clustering problem. Apart from the rank condition,
our BM factorization (6) only involves one low-rank matrix variable $R,$ while
(27) also involves a large high rank positive semidefinite matrix variable
$Y$111111Note that $RR^{\top}+Y=\alpha I.$ If $r$ is small, then $Y$ must by
high rank., which would lead to expensive computation when one has to handle
such a high-dimensional matrix variable in any algorithm for solving (27).
Moreover, we are not aware of an existing algorithm that can handle the
problem (27) efficiently.

##  3 An algebraic variety

In this section, we consider the case where the mapping
$\mathcal{A}(\cdot)={\rm diag}(\cdot)$. The constraints of (6) come from two
parts, one of them is the following set,

$\mathcal{B}_{n,r}=\left\\{R\in\mathbb{R}^{n\times r}:\ {\rm
diag}\big{(}RR^{\top}\big{)}=e,\ R^{\top}e=0\right\\}.$ (32)

We remind the reader that ${\rm diag}(RR^{\top})=e$ is equivalent to
$\|e_{i}^{\top}R\|_{2}^{2}=1$ for all $i\in[n].$ We call a point
$R\in\mathcal{B}_{n,r}$ smooth (nonsingular) point if the linear independent
constraint qualification (LICQ) holds at $R$. Otherwise we say that $R$ is
nonsmooth (singular). Strictly speaking, LICQ is a sufficient but not
necessary condition for a point to be smooth. However, we identify smoothness
and LICQ because LICQ is convenient for us to use in optimization.

Observe that $\mathcal{B}_{n,r}$ is an algebraic variety defined as the common
zeros of $n+r$ polynomials (see [12, 17]). Its tangent space at a smooth point
$R\in\mathcal{B}_{n,r}$ is given by

$T_{R}\mathcal{B}_{n,r}:=\left\\{H\in\mathbb{R}^{n\times r}:\ H^{\top}e=0,\
{\rm diag}(HR^{\top})=0\right\\}.$ (33)

###  3.1 Local geometric properties of $\mathcal{B}_{n,r}$

From now on, we let $E^{n}:=\left\\{v\in\\{-1,1\\}^{n}:\
e^{\top}v=0\right\\}.$ Note that $\mathcal{B}_{n,r}$ is not always a smooth
manifold and we have the following decomposition result.

######  Proposition 3.1.

$\mathcal{B}_{n,r}=\mathcal{B}_{n,r}^{0}\sqcup\mathcal{B}_{n,r}^{1}$, where
$\sqcup$ denotes the disjoint union of the sets
$\mathcal{B}_{n,r}^{0}:=\left\\{vb^{\top}:\ v\in E^{n},\ b\in{\rm
S}^{r-1}\right\\}$ and $\mathcal{B}^{1}_{n,r}:=\left\\{R\in\mathcal{B}_{n,r}:\
{\rm rank}\left(R\right)\geq 2\right\\}.$ Moreover, $\mathcal{B}_{n,r}^{0}$
contains the singular points of $\mathcal{B}_{n,r}$ and
$\mathcal{B}_{n,r}^{1}$ contains the smooth points of $\mathcal{B}_{n,r}.$
Note that when $n$ is odd, $\mathcal{B}^{0}_{n,r}=\emptyset.$

###### Proof.

The decomposition comes from the fact that any $R\in\mathcal{B}_{n,r}$ with
rank equal to 1 has the form $R=vb^{\top},$ for some $v\in E^{n}$ and
$b\in{\rm S}^{r-1}$. For any $R\in\mathcal{B}_{n,r}$, define the linear
mapping $\mathcal{A}_{R}:\mathbb{R}^{n\times r}:\
\rightarrow\mathbb{R}^{r}\times\mathbb{R}^{n},$ such that for any
$H\in\mathbb{R}^{n\times r}$, $\mathcal{A}_{R}(H)=\left(H^{\top}e,{\rm
diag}(HR^{\top})\right)$. Then for any
$(\lambda,\mu)\in\mathbb{R}^{r}\times\mathbb{R}^{n}$,
$\mathcal{A}_{R}^{*}\left((\lambda,\mu)\right)=e\lambda^{\top}+{\rm
diag}(\mu)R.$ It is easy to see that LICQ is equivalent to that
$\mathcal{A}_{R}^{*}$ is an injection. For
$R=vb^{\top}\in\mathcal{B}_{n,r}^{0},$ we have that
$\mathcal{A}_{R}^{*}\left((-b,v)\right)=0.$ So, $R$ is singular. For any
$R\in\mathcal{B}_{n,r}^{1}$, we only have to prove that $\mathcal{A}_{R}^{*}$
is an injection. Suppose
$\mathcal{A}_{R}^{*}\left((\lambda,\mu)\right)=e\lambda^{\top}+{\rm
diag}(\mu)R=0$. If there exists $i\in[n]$ such that $\mu_{i}=0$, then
$\lambda=0$ and ${\rm diag}(\mu)R=0.$ Because ${\rm diag}(RR^{\top})=e$, we
have $\mu=0.$ If for all $i\in[n]$, $\mu_{i}\neq 0$, then $R=-{\rm
diag}(\mu)^{-1}ea^{\top}$, which is rank 1, and we get a contradiction.
Therefore, $(\lambda,\mu)=0$ and $\mathcal{A}_{R}^{*}$ is injective. We have
that $R$ is a smooth point of $\mathcal{B}_{n,r}$.

Proposition 3.1 shows that the only singular points in $\mathcal{B}_{n,r}$ are
rank-1 points. We should note that when using the Burer and Monteiro
factorization, $R\in\mathbb{R}^{n\times r}$ is equivalent to $RQ$ for any
orthogonal $Q\in{\rm St}(r,r).$ This is because
$RR^{\top}=RQQ^{\top}R^{\top}.$ Thus, for any singular point $R=vb^{\top}$
such that $v\in E^{n}$ and $b\in{\rm S}^{r-1},$ there exists an orthogonal
matrix $Q\in{\rm St}(r,r)$ such that $RQ=ve_{1}^{\top}.$ This implies that
there are only finitely many classes of singular points. Moreover, the number
of classes is $\frac{1}{2}\binom{n}{\frac{n}{2}}$ if $n$ is even and $0$ if
$n$ is odd. For convenience, we only consider singular points of the form
$R=ae_{1}^{\top}$ for some $a\in E^{n}.$

Now we move on to characterize the tangent cone and second order tangent set
at a singular point of $\mathcal{B}_{n,r}$ in the following two propositions.
We defer their proofs to Appendix C.

######  Proposition 3.2.

Suppose $r>1$ and $R=ae_{1}^{\top}\in\mathcal{B}_{n,r}$ is a singular point
for some $a\in E^{n}$, then its tangent cone is given by:

$\mathcal{T}^{i}_{\mathcal{B}_{n,r}}(R)=\mathcal{T}_{\mathcal{B}_{n,r}}(R)=\widetilde{\mathcal{T}}:=\big{\\{}[0,H_{1}]:H_{1}\in\mathbb{R}^{n\times(r-1)},e^{\top}H_{1}=0,a^{\top}{\rm
diag}\big{(}H_{1}H_{1}^{\top}\big{)}=0\big{\\}}.$ (34)

######  Proposition 3.3.

Suppose $r>1.$ Let $R=ae_{1}^{\top}$ be a singular point for some $a\in
E^{n}.$ Then for any
$H=[0,H_{1}]\in\mathcal{T}^{i}_{\mathcal{B}_{n,r}}(R)=\mathcal{T}_{\mathcal{B}_{n,r}}(R).$

  * (i)

If $H_{1}\notin a\mathbb{R}^{1\times(r-1)},$ then

$\displaystyle\mathcal{T}^{i,2}_{\mathcal{B}_{n,r}}(R,H)=\mathcal{T}^{2}_{\mathcal{B}_{n,r}}(R,H)=\widetilde{\mathcal{T}}^{2}(R,H)$
$\displaystyle:=\Big{\\{}\left[-a\circ{\rm
diag}\big{(}H_{1}H_{1}^{\top}\big{)},W_{1}\right]:\ a^{\top}{\rm
diag}\big{(}H_{1}W^{\top}_{1}\big{)}=0,\,e^{\top}W_{1}=0,\,W_{1}\in\mathbb{R}^{n\times(r-1)}\Big{\\}}.$

  * (ii)

If $H_{1}=a\lambda^{\top}$ for some $\lambda\in\mathbb{R}^{r-1},$ then

$\displaystyle\mathcal{T}^{i,2}_{\mathcal{B}_{n,r}}(R,H)=\mathcal{T}^{2}_{\mathcal{B}_{n,r}}(R,H)=\widetilde{\mathcal{T}}^{2}(R,H)$
$\displaystyle:=\Big{\\{}\left[-a\circ{\rm
diag}\big{(}H_{1}H_{1}^{\top}\big{)},W_{1}\right]:\ a^{\top}{\rm
diag}\big{(}W_{1}W^{\top}_{1}\big{)}=0,\,e^{\top}W_{1}=0,\,W_{1}\in\mathbb{R}^{n\times(r-1)}\Big{\\}}.$

Proposition 3.2 and 3.3 together imply that $\mathcal{B}_{n,r}$ is second
order directionally differentiable at every singular point (see definition
3.32 of [5]). With these geometric properties, we can move on to study
operations on $\mathcal{B}_{n,r},$ which are crucial for designing algorithms.
We will first discuss some important operations on the smooth part of
$\mathcal{B}_{n,r}$ such as projection and retraction in the next two
sections, while the discussion on how to handle the singular points will be
given later in Appendix B since we rarely encounter singular points in
practical computations.

###  3.2 Projection onto tangent spaces of smooth points on
$\mathcal{B}_{n,r}$

For notational convenience, let ${\rm Proj}_{R}(C):={\rm
Proj}_{T_{R}\mathcal{B}_{n,r}}(C).$

######  Theorem 3.4.

For any smooth point $\bar{R}\in\mathcal{B}_{n,r}$ and
$C\in\mathbb{R}^{n\times r},$

${\rm Proj}_{\bar{R}}(C)=J\left(C-{\rm
diag}\left(\lambda_{C}^{\bar{R}}\right)\bar{R}\right),$ (35)

where

$\lambda_{C}^{\bar{R}}:=\Big{(}I+\bar{R}\left(I-\bar{R}^{\top}\bar{R}/n\right)^{-1}\bar{R}^{\top}/n\Big{)}{\rm
diag}\big{(}JC\bar{R}^{\top}\big{)}.$ (36)

###### Proof.

Suppose $\bar{R}\in\mathcal{B}_{n,r}^{1}$, we have

${\rm Proj}_{\bar{R}}(C):=\arg\min_{H\in\mathbb{R}^{n\times
r}}\left\\{\frac{1}{2}\|H-C\|^{2}:\ H^{\top}e=0,\ {\rm
diag}(H\bar{R}^{\top})=0\right\\}.$ (37)

Since, the feasible region in (37) is convex and the objective function is
strongly convex, it has a unique optimal solution, say $H.$ From LICQ property
of $R$, $H$ also satisfies LICQ. Then $H$ must be a KKT solution, that is,
there exists $\lambda_{1}\in\mathbb{R}^{r}$, $\lambda_{2}\in\mathbb{R}^{n}$
such that

$\displaystyle H-C+e\lambda_{1}^{\top}+{\rm diag}(\lambda_{2})\bar{R}=0$ (38)
$\displaystyle H^{\top}e=0$ (39) $\displaystyle{\rm diag}(H\bar{R}^{\top})=0$
(40)

Multiplying $J$ to the left hand side of (38), and noting that $e^{\top}H=0$,
we have

$H=J(C-{\rm diag}(\lambda_{2})\bar{R}).$

Plug this into (40), we have

${\rm diag}\big{(}J{\rm diag}(\lambda_{2})\bar{R}\bar{R}^{\top}\big{)}={\rm
diag}\big{(}JC\bar{R}^{\top}\big{)}.$ (41)

Apply the following formula to (41)

${\rm diag}\left(A{\rm diag}(\lambda)B\right)=\big{(}A\circ
B^{\top}\big{)}\lambda,$ (42)

we have

$\left(I-\bar{R}\bar{R}^{\top}/n\right)\lambda_{2}={\rm
diag}\big{(}JC\bar{R}^{\top}\big{)}$ (43)

Since $\bar{R}\bar{R}^{\top}/n\succeq 0$ and ${\rm
Tr}(\bar{R}\bar{R}^{\top}/n)={\rm Tr}(I/n)=1,$ we have
$\|\bar{R}\bar{R}^{\top}/n\|_{2}\leq 1$ and the equality holds if and only if
$\bar{R}$ is rank 1. Because ${\rm rank}(\bar{R})\geq 2$, we must have
$I-\bar{R}\bar{R}^{\top}/n\succ 0.$ From the Sherman-Morrison-Woodbury
formula, we have

$\lambda_{2}=\Big{(}I+\bar{R}\left(I-\bar{R}^{\top}\bar{R}/n\right)^{-1}\bar{R}^{\top}/n\Big{)}{\rm
diag}\big{(}JC\bar{R}^{\top}\big{)},$ (44)

which is exactly $\lambda_{C}^{\bar{R}}.$

###  3.3 Retraction around smooth points of $\mathcal{B}_{n,r}$

For $\bar{R}\in\mathcal{B}_{n,r}$ and $H\in T_{\bar{R}}\mathcal{B}_{n,r}$,
define

$\mathcal{R}_{\bar{R}}(H):={\rm
Proj}_{\mathcal{B}_{n,r}}(\bar{R}+H)\in\arg\min\left\\{\|R-(\bar{R}+H)\|^{2}/2:\
R\in\mathcal{B}_{n,r}\right\\}.$ (45)

Since $\mathcal{B}_{n,r}$ is compact, the optimization problem (45) always has
an optimal solution, but it may not be unique as $\mathcal{B}_{n,r}$ is non-
convex. Now, we focus on the following optimization problem. Given
$Y\in\mathbb{R}^{n\times r},$

$\min\Big{\\{}\frac{1}{2}\|R-Y\|^{2}:\ R\in\mathcal{B}_{n,r}\Big{\\}}.$ (46)

For $Y_{1},Y_{2}\in\mathbb{R}^{n\times r}$, we say $Y_{1}\sim Y_{2}$ if there
exists $b\in\mathbb{R}^{r}$ such that $Y_{1}=Y_{2}+eb^{\top}.$ It is easy to
check that $\sim$ is an equivalence relation, which partitions
$\mathbb{R}^{n\times r}$ into equivalence classes
$[Y]:=\left\\{Y_{1}\in\mathbb{R}^{n\times r}:\ Y_{1}\sim Y\right\\}.$

######  Lemma 3.5.

For $Y\in\mathbb{R}^{n\times r}$, and $Y_{1}\in[Y],$

$\arg\min\left\\{\|R-Y\|^{2}/2:\
R\in\mathcal{B}_{n,r}\right\\}=\arg\min\left\\{\|R-Y_{1}\|^{2}/2:\
R\in\mathcal{B}_{n,r}\right\\}.$ (47)

###### Proof.

Since $Y_{1}\in[Y]$, $Y_{1}=Y+eb^{\top}$ for some $b\in\mathbb{R}^{r}$.
Because for any $R\in\mathcal{B}_{n,r}$, $R^{\top}e=0$, the two objective
functions in (47) only differs by a constant. Thus, the optimal solutions set
are the same.

Now, we move on to show that problem (46) is related to the well known
geometric median problem. This is the key point for solving (46) efficiently.
For $Y\in\mathbb{R}^{n\times r}$ such that every row is non-zero, define the
mapping $\mathcal{O}:\mathbb{R}^{n\times r}\rightarrow\mathbb{R}^{n\times r}$
such that $\mathcal{O}(Y):=Y\left({\rm
Diag}(YY^{\top})\right)^{-\frac{1}{2}},$ which is derived from $Y$ by
normalizing each row of $Y$ to have unit $2$-norm. For any
$Y\in\mathbb{R}^{n\times r},$ define

$\mathcal{F}(Y):=\\{\mathcal{O}(Y_{1}):\,Y_{1}\in[Y]\\}\,\cap\,\mathcal{B}_{n,r}.$

######  Lemma 3.6.

For any $Y\in\mathbb{R}^{n\times r}$, either $\mathcal{F}(Y)=\emptyset$ or
$\mathcal{F}(Y)=\arg\min\left\\{\frac{1}{2}\|R-Y\|^{2}:\,R\in\mathcal{B}_{n,r}\right\\}$
and $|\mathcal{F}(Y)|=1.$

###### Proof.

Suppose $\mathcal{F}(Y)\neq\emptyset$. Then there exists $Y_{1}\in[Y]$ such
that $\mathcal{O}(Y_{1})\in\mathcal{B}_{n,r}.$ Because
$\mathcal{O}(Y_{1})={\rm Proj}_{{\rm OB}_{n,r}}(Y_{1})$ and
$\mathcal{B}_{n,r}\subset{\rm OB}_{n,r},$ we have
$\\{\mathcal{O}(Y_{1})\\}=\arg\min\left\\{\frac{1}{2}\|R-Y_{1}\|^{2}:\,R\in\mathcal{B}_{n,r}\right\\}.$
From Lemma 3.5,
$\\{\mathcal{O}(Y_{1})\\}=\arg\min\left\\{\frac{1}{2}\|R-Y\|^{2}:\,R\in\mathcal{B}_{n,r}\right\\}.$
From the arbitrariness of $Y_{1}$, we have
$\mathcal{F}(Y)\subset\arg\min\left\\{\frac{1}{2}\|R-Y\|^{2}:\,R\in\mathcal{B}_{n,r}\right\\}$.
Thus,

$\mathcal{O}(Y_{1})\in\mathcal{F}(Y)\subset\arg\min\Big{\\{}\frac{1}{2}\|R-Y\|^{2}:\,R\in\mathcal{B}_{n,r}\Big{\\}}=\\{\mathcal{O}(Y_{1})\\},$

and the proof is completed.

From Lemma 3.6, we know that if $\mathcal{F}(Y)\neq\emptyset$, then
$\mathcal{F}(Y)={\rm Proj}_{\mathcal{B}_{n,r}}(Y).$ The next proposition gives
a sufficient and necessary condition for $\mathcal{F}(Y)\neq\emptyset.$

######  Proposition 3.7.

For any $Y\in\mathbb{R}^{n\times r}$, $\mathcal{F}(Y)\neq\emptyset$ if and
only if the following geometric median problem,

$\min\Big{\\{}f(b):=\sum_{i=1}^{n}\|y_{i}-b\|_{2}:\
b\in\mathbb{R}^{r}\Big{\\}},$ (48)

has an optimal solution $b_{0}\notin\\{y_{i}:i\in[n]\\}.$ Here, $y_{i}^{\top}$
is the $i$th row of $Y.$

###### Proof.

Note that (48) is a convex optimization problem, it has an optimal solution
$b_{0}\notin\\{y_{i}:i\in[n]\\}$ if and only if

$\nabla
f(b)\Big{|}_{b=b_{0}}=\sum_{i=1}^{n}\frac{b_{0}-y_{i}}{\|b_{0}-y_{i}\|}=0.$

The latter condition is equivalent to $\mathcal{F}(Y)\neq\emptyset.$

From Proposition 3.7 and Lemma 3.6, we may solve (46) by solving (48). Thus, a
nonconvex problem is reduced to a well studied convex problem. There are many
algorithms to solve (48) with guaranteed convergence, one simple method is the
Weiszfeld algorithm ([41]), which performs the following iteration:

$b_{i+1}=\left(\sum_{j=1}^{n}\frac{y_{j}}{\|b_{i}-y_{j}\|}\right)\Big{/}\left(\sum_{j=1}^{n}\frac{1}{\|b_{i}-y_{j}\|}\right).$
(49)

However, the condition $b_{0}\notin\\{y_{i}:i\in[n]\\}$ in Proposition 3.7 may
not always hold. Consider $Y\in\mathbb{R}^{3\times 2}$, whose rows are not
collinear (not on the same line), it is well known that the geometric median
of $y_{1},y_{2},y_{3}$ belongs to $\\{y_{i}:\ i\in[3]\\}$ if and only if the
triangle $y_{1}y_{2}y_{3}$ has an angle $\geq 120^{\circ}.$ This means that
the condition $b_{0}\notin\\{y_{i}:i\in[n]\\}$ may not hold with probability
1. Moreover, the following example show that one cannot find a neighbourhood
of $\mathcal{B}_{n,r}$ such that every $Y$ in it satisfies
$\mathcal{F}(Y)\neq\emptyset.$

######  Example 3.8.

Consider $\mathbb{R}^{4\times 2}.$ For any $\epsilon>0$, define
$Y_{\epsilon}\in\mathbb{R}^{4\times 2}$ such that

$Y_{\epsilon}:=\begin{pmatrix}0&0&0&-\frac{\epsilon}{2}\\\
1&1-\frac{\epsilon}{2}&-1&-1+\frac{\epsilon}{2}\end{pmatrix}^{\top}.$

We have ${\rm dist}\left(Y_{\epsilon},\mathcal{B}_{n,r}\right)<\epsilon$.
Because the first three rows of $Y_{\epsilon}$ are on the same line, the only
geometric median of its rows vector is $(0,1-\frac{\epsilon}{2})^{\top}$, so
$\mathcal{F}\left(Y_{\epsilon}\right)=\emptyset.$ However, as
$\epsilon\rightarrow 0$, $Y_{\epsilon}$ can be arbitrarily close to
$\mathcal{B}_{4,2}.$

The following theorem is a sufficient condition for the geometric median
method to be useful.

######  Theorem 3.9.

Suppose $Y\in\mathbb{R}^{n,r}$ satisfies $\|Y\|_{2}<(1-\beta)\sqrt{n}$ for
some $0<\beta<1,$ and ${\rm
dist}\left(Y,\mathcal{B}_{n,r}\right)<\frac{\beta\sqrt{n}}{\sqrt{n}+1}<1.$
Then $\mathcal{F}(Y)\neq\emptyset.$

###### Proof.

Let $R$ be an optimal solution of (46). Then $\|R-Y\|_{F}={\rm
dist}\left(Y,\mathcal{B}_{n,r}\right).$ Let $\|R-Y\|_{F}=\alpha\sqrt{n}.$ Then
we have $\alpha<\frac{\beta}{\sqrt{n}+1}$ and hence

$\beta-\alpha>\alpha\sqrt{n}.$

Moreover
$\|R\|_{2}\leq\|Y\|_{2}+\|Y-R\|_{F}<\left(1-\beta+\alpha\right)\sqrt{n}.$
Since any rank-one solution of (46) must have spectral norm equal to
$\sqrt{n}$, the last inequality implies that ${\rm rank}\left(R\right)\geq 2$
and

$\Big{\|}\left(I-RR^{\top}/n\right)^{-1}\Big{\|}_{2}\leq\frac{1}{1-\|RR^{\top}/n\|_{2}}<\frac{1}{1-\left(1-\beta+\alpha\right)^{2}}.$

By Proposition 3.1, $R$ is smooth. So it is a KKT solution of (46). There
exist $\mu\in\mathbb{R}^{r}$, $\lambda\in\mathbb{R}^{n}$ such that the
following system holds:

$\displaystyle R-Y+e\mu^{\top}+{\rm diag}(\lambda)R=0$ (50) $\displaystyle
R^{\top}e=0$ (51) $\displaystyle{\rm diag}(RR^{\top})=e$ (52)

Multiplying $J$ to the left hand side of (50) and using (51), we have

$R=J\left(Y-{\rm diag}(\lambda)R\right).$

Then

${\rm diag}(RR^{\top})={\rm diag}\big{(}J\left(Y-{\rm
diag}(\lambda)R\right)R^{\top}\big{)}.$

Using (52) and (42), we have

$\left(I-RR^{\top}/n\right)\lambda={\rm diag}\big{(}JYR^{\top}\big{)}-e.$

Thus we have that

$\lambda=\left(I-RR^{\top}/n\right)^{-1}\big{(}{\rm
diag}(JYR^{\top})-e\big{)}=\left(I-RR^{\top}/n\right)^{-1}{\rm
diag}\big{(}J(Y-R)R^{\top}\big{)}.$

Hence

$\displaystyle\|\lambda\|_{2}\leq\Big{\|}\left(I-RR^{\top}/n\right)^{-1}\Big{\|}_{2}\,\big{\|}{\rm
diag}\big{(}J(Y-R)R^{\top}\big{)}\big{\|}_{2}\leq\frac{\big{\|}{\rm
diag}\left(J(Y-R)R^{\top}\right)\big{\|}_{2}}{1-\left(1-\beta+\alpha\right)^{2}}.$
(53)

Note that we also have

$\displaystyle\Big{\|}{\rm
diag}\big{(}J(Y-R)R^{\top}\big{)}\Big{\|}_{2}\leq\big{\|}{\rm
diag}\big{(}(Y-R)R^{\top}\big{)}\big{\|}_{2}+\big{\|}R(Y-R)^{\top}e/n\big{\|}_{2}$
$\displaystyle\leq\|Y-R\|_{F}+\|R\|_{2}\|Y-R\|_{2}\|e\|_{2}/n<\left(2-\beta+\alpha\right)\alpha\sqrt{n},$
(54)

where we used Cauchy-Schwarz inequality and ${\rm
diag}\left(RR^{\top}\right)=e$ in the second inequality of (54). Plug (54)
into (53), we have

$\|\lambda\|_{2}\leq\frac{\alpha\sqrt{n}}{\beta-\alpha}<1.$

Now, from (50), we have $\left(I+{\rm
diag}\left(\lambda\right)\right)R=Y-e\mu^{\top}.$ Since $\|\lambda\|_{2}<1$,
we have $\left(I+{\rm diag}\left(\lambda\right)\right)\succ 0$. Hence
$R=\mathcal{O}(Y-e\mu^{\top})$ and so $\mathcal{F}(Y)\neq\emptyset.$

######  Remark 3.10.

One may get the impression that the result in Theorem 3.9 contradicts Example
3.8. Intuitively speaking, Theorem 3.9 requires that a point’s distance to
$\mathcal{B}_{n,r}$ is smaller than its distance to a rank-1 matrix up to some
constant, which is not equivalent to saying that the matrix is sufficiently
close to $\mathcal{B}_{n,r}.$ Theorem 3.9 is useful for designing algorithms
for solving the $k-$equipartition problem GEP1. Consider the problem in (45).
In practice, $R$ is usually far from being rank-1. Moreover, for $k\geq 3$,
the $k-$equipartition problem has the constraints
$\|R\|_{2}\leq\sqrt{\frac{n}{k-1}}.$ If we choose
$\beta\in\big{(}0,\frac{\sqrt{n}+1}{\sqrt{n}+2}\big{(}1-\frac{1}{\sqrt{k-1}}\big{)}\big{)}$
in Theorem 3.9 and $H$ to be sufficiently small such that ${\rm
dist}\left(R+H,\mathcal{B}_{n,r}\right)\leq\|H\|_{F}<\frac{\beta\sqrt{n}}{\sqrt{n}+1},$
then we have

$\|Y\|_{2}=\|R+H\|_{2}\leq\|R\|_{2}+\|H\|_{2}=\sqrt{\frac{n}{k-1}}+\frac{\beta\sqrt{n}}{\sqrt{n}+1}<(1-\beta)\sqrt{n}.$

This means that, if ${\rm rank}\left(R\right)\geq 2$, we can first compute the
retraction $\mathcal{R}_{R}(H)$ by solving (48), with $Y=R+H$. If
$\mathcal{F}(Y)\neq\emptyset$, then we have computed the retraction
successfully, otherwise, we can perform a backtracking step: $H\leftarrow
H/\sigma$ for some $\sigma>1$. Theorem 3.9 tells us that after a finite number
of backtracking steps, $\mathcal{F}(Y)\neq\emptyset$ and the retraction
computation will be successful. Note that the condition ${\rm
rank}\left(R\right)\geq 2$ is usually satisfied in our experiments. Moreover,
for $k\geq 3$, any feasible $R$ must satisfy ${\rm rank}(R)\geq 2$ since we
have $\|R\|_{2}\leq\sqrt{\frac{n}{k-1}}<\sqrt{n}$. The last inequality rules
out the possibility for $R$ to have rank one since in that case its spectral
norm must be equal to $\sqrt{n}$.

Up to now, we are able to compute the retraction of a point only if it is
“smooth”. It is an interesting mathematical problem to compute the retraction
in the general case. Suppose $Y\in\mathbb{R}^{n\times r}=vb^{\top}$ such that
$b\in{\rm S}^{r-1}$ is a rank-1 matrix. Without loss of generality, assume
that the entries of $v$ are arranged in an ascending order: $v_{1}\leq
v_{2}\leq\ldots,\leq v_{n}$. If $n$ is an even number, then ${\rm
Proj}_{\mathcal{B}_{n,r}}\left(Y\right)=wb^{\top},$ where
$w_{1}=\ldots=w_{n/2}=-1$, $w_{n/2+1}=\ldots=w_{n}=1.$ This can be proved by
first assuming $v_{1}<v_{2}<\ldots<v_{n}$, in which case we can apply
Proposition 3.7. Then use the continuity of ${\rm
Proj}_{\mathcal{B}_{n,r}}\left(Y\right)$ to prove the case for $v_{1}\leq
v_{2}\leq\ldots\leq v_{n}.$

##  4 Algorithms

In this section, we consider algorithms to solve (56), with
$\mathcal{A}(\cdot)={\rm diag}(\cdot)$ and $\alpha=\frac{n}{k-1}.$ We will
discuss the algorithms for minimum bisection and graph multi-equipartition
seperately.

###  4.1 Minimum bisection

In this subsection, we consider the case $k=2$ and hence $\alpha=n$. In this
case, since ${\rm diag}(RR^{\top})=e$ implies that $\|R\|_{F}=\sqrt{n}$, the
constraint $\|R\|_{2}\leq\sqrt{\alpha}$ is redundant and it can be removed. We
rewrite the minimum bisection problems as follows:

$\min\Big{\\{}\left\langle C,X\right\rangle:\ {\rm diag}\left(X\right)=e,\
Xe=0,\ \left\langle X,ee^{\top}\right\rangle=0,\
X\in\mathbb{S}^{n}_{+}\Big{\\}},$ (55)
$\min\Big{\\{}f(R):=\frac{1}{2}\left\langle C,RR^{\top}\right\rangle:\
R\in\mathcal{B}_{n,r}\Big{\\}}.$ (56)

From section 3, we know that $\mathcal{B}_{n,r}$ can be viewed as an affine
variety with finitely many classes of singular points. In particular, for any
$\delta\in(0,1),$ $\mathcal{B}_{n,r}$ can be separated as
$\mathcal{B}_{n,r}^{\delta-}$ and $\mathcal{B}_{n,r}^{\delta+}$ such that

$\mathcal{B}_{n,r}^{\delta-}:=\left\\{R\in\mathcal{B}_{n,r}:\
\|R\|_{2}<\sqrt{1-\delta}\sqrt{n}\right\\},$ (57)
$\mathcal{B}_{n,r}^{\delta+}:=\Big{\\{}R\in\mathcal{B}_{n,r}:\
\|R\|_{2}\geq\sqrt{1-\delta}\sqrt{n}\Big{\\}}.$ (58)

For any $\delta\in(0,1),$ all singular points are contained in
$\mathcal{B}_{n,r}^{\delta+}.$ If there exists some constant $\delta\in(0,1)$
such that all iterations are in $\mathcal{B}_{n,r}^{\delta-},$ then we can use
any kind of Riemannian optimization algorithm to solve (56). Actually, this is
often the case in practice. Only in very rare situation will one encounter a
singular point. For example, if the optimal solution of (56) is rank-1, then
the iterations will approach a singular point. Note that our escaping strategy
in Appendix B only works for an exact singular point. Thus, we will use a
rounding procedure that will round any “nearly singular point” to a singular
point and then apply Theorem B.3 to handle it. For $\delta<\frac{1}{2},$
define the function
$\texttt{Round}:\mathcal{B}_{n,r}^{\delta+}\rightarrow\mathbb{R}^{n\times r}$
such that for any $R\in\mathcal{B}_{n,r}^{\delta+},$
$\texttt{Round}(R):=\texttt{sgn}(u)e_{1}^{\top},$ where $u\in\mathbb{R}^{n}$
is the singular vector (with its first nonzero component being positive) that
corresponds to the largest singular value of $R$ and sgn is the sign function.
Note that since $\|R\|_{F}^{2}=n$ for $R\in\mathcal{B}_{n,r},$ so when
$\delta<\frac{1}{2},$ the leading singular value has multiplicity one and $u$
must be unique. Hence $\texttt{Round}(.)$ is well-defined. We now state our
algorithm.

Algorithm 1:

Choose $R_{0}\in\mathcal{B}_{n,r}$ and $\delta_{0}\in(0,\frac{1}{2}).$ Set
$k=0.$

Step 1.

Set $R_{k}$ as the initial point, $f_{k}:=f(R_{k}).$ Use any Riemannian
optimization method to solve (56) such that every iteration
$R\in\mathcal{B}_{n,r}$ satisfies $f(R)\leq f_{k}$. If for some iteration, we
obtain a point $\widehat{R}_{k}\in\mathcal{B}_{n,r}^{\delta_{k}+},$ goto Step
2.

Step 2.

Case 1. $\texttt{Round}\big{(}\widehat{R}_{k}\big{)}\notin\mathcal{B}_{n,r},$
set $R=\widehat{R}_{k},$ goto Step 3.
Case 2. $\texttt{Round}\big{(}\widehat{R}_{k}\big{)}\in\mathcal{B}_{n,r}$ and
is optimal, stop.
Case 3. $\texttt{Round}\big{(}\widehat{R}_{k}\big{)}\in\mathcal{B}_{n,r}$ and
is non-optimal: Use the escaping strategy in Appendix B to find another smooth
point $R^{+}\in\mathcal{B}_{n,r}$ such that
$f(R^{+})<f(\texttt{Round}(\widehat{R}_{k}))-{\epsilon\big{(}\texttt{Round}(\widehat{R}_{k})\big{)}},$
where $\epsilon\big{(}\texttt{Round}(\widehat{R}_{k})\big{)}>0$ is only
related to $\texttt{Round}(\widehat{R}_{k})$. If
$f(R^{+})<f(\widehat{R}_{k})$, set $R\leftarrow R^{+};$ Otherwise, set
$R\leftarrow\widehat{R}_{k}.$121212In this case, $R$ is also a smooth point.
This is because otherwise
$f(R^{+})<f(\texttt{Round}(\widehat{R}_{k}))=f(\widehat{R}_{k})$ and this case
won’t happen. Goto Step 3.

Step 3.

Set $\delta_{k+1}\leftarrow\delta_{k}/2$, $R_{k+1}\leftarrow R$, $k\leftarrow
k+1,$ goto Step 1.

In Step 1, many Riemannian optimization methods with linesearch can satisfy
$f(R)\leq f_{k}$ for any iteration $R.$ In practice, we will use a Riemannian
gradient descent method with Barzilai-Borwein (BB) step and nonmonotone
linesearch (see [2, 15, 19, 22, 30]) in Step 1 of Algorithm 1. For
$R\in\mathcal{B}_{n,r}$, the Riemannian gradient of a smooth function $f$ is
given by ${\rm grad}f(R):={\rm Proj}_{R}\left(\nabla f(R)\right).$

Riemannian gradient BB method:

Given $R_{0}\in\mathcal{B}_{n,r},\alpha_{0}$, integer $M\geq 0$,
$\gamma\in(0,1)$, $\beta>0,$ $\sigma>0$, $\epsilon>0$. Set $k=0.$

[Step 1] Compute the Riemannian gradient $g_{k}={\rm grad}f(R_{k})$. If
$\|g_{k}\|=0$, stop.

[Step 2] If $\alpha_{k}\leq\epsilon$ or $\alpha_{k}\geq 1/\epsilon,$ then set
$\alpha_{k}=\beta$

[Step 3] Set $\tau=1/\alpha_{k}.$

[Step 4] (nonmonotone line search)
If $f\big{(}{\rm Proj}_{\mathcal{B}_{n,r}}(R_{k}-\tau
g_{k})\big{)}\leq\max_{0\leq
j\leq\min\\{k,M\\}}(f_{k-j})-\gamma\tau\big{\langle}g_{k},\,g_{k}\big{\rangle}$,
set $\tau_{k}=\tau$, $R_{k+1}={\rm
Proj}_{\mathcal{B}_{n,r}}(R_{k}-\tau_{k}g_{k}),$ and goto Step 6; Otherwise,
set $\tau\leftarrow\sigma\tau$ and repeat Step 4.

[Step 5] Let $y_{k}=g_{k}-{\rm Proj}_{R_{k}}(g_{k-1}).$ Set
$\alpha_{k+1}=|\big{\langle}g_{k},\,y_{k}\big{\rangle}|/(\tau_{k}\big{\langle}g_{k},\,g_{k}\big{\rangle})$
or
$\big{\langle}y_{k},\,y_{k}\big{\rangle}/(\tau_{k}\big{\langle}g_{k},\,g_{k}\big{\rangle}),$
and goto Step 1.

In our code, we set $\gamma=10^{-4},\epsilon=10^{-10},\sigma=0.5$. For the
parameter $\delta$, we set

$\beta=\begin{cases}1&\|g_{k}\|_{F}>1,\\\ \|g_{k}\|_{F}^{-1}\
&10^{-5}\leq\|g_{k}\|_{F}\leq 1\\\ 10^{5}\ &\
\|g_{k}\|_{F}<10^{-5}\end{cases}.$ (59)

Given a tolerance ${\rm tol}>0$, we stop our algorithm when ${\|{\rm
grad}f(R_{k})\|_{F}}/({1+\|R_{k}\|_{F}})<{\rm tol}$.

Note that although we use nonmonotone linesearch in the above algorithm, the
condition $f(R)\leq f_{k}$ is still satisfied in Algorithm 1. The next theorem
shows that, our algorithm terminates in finite iterations131313Here finite
iteration means that the outer loop of algorithm 1 is finite. We consider the
whole Riemannian optimization method in step 1 as part of one outer
iteration..

######  Theorem 4.1.

For any initial point $R_{0}\in\mathcal{B}_{n,r}$ and any
$\delta_{0}\in(0,1),$ Algorithm 1 will terminate in finite iterations.

###### Proof.

Assume by contradiction that Algorithm 1 doesn’t terminate. Then we have that
$\delta_{k}\rightarrow 0.$ From Lemma A.2, for $k$ sufficiently large,
$\texttt{Round}\big{(}\widehat{R}_{k}\big{)}\in\mathcal{B}_{n,r}.$ Thus Case 1
and Case 2 in Step 2 will not be triggered when $k$ is large enough. Since
there are finite many singular points inside $E^{n}e_{1}^{\top}.$ One of them
will be visited infinitely many times, that is, there exists $a\in E^{n}$ and
$\\{i_{k}\\}_{k\in\mathbb{N}}\subset\mathbb{N}$ such that
$\texttt{Round}\big{(}\widehat{R}_{i_{k}}\big{)}=ae_{1}^{\top}$ for any
$k\in\mathbb{N}.$ Also, $ae_{1}^{\top}$ is non-optimal. From (69) in Lemma
A.2, we have that $f\big{(}\widehat{R}_{i_{k}}\big{)}\rightarrow
f\left(ae_{1}^{\top}\right).$ Note that from the condition $f(R)\leq f_{k}$ in
Step 1 and Case 3 in Step 2, we have that $f_{k+1}\leq f_{k}.$ Therefore, for
$k$ sufficiently large, $f_{k}\leq f(ae_{1}^{\top})-\epsilon(ae_{1}^{\top}).$
This contradicts that $f\big{(}\widehat{R}_{i_{k}}\big{)}\rightarrow
f\left(ae_{1}^{\top}\right).$

###  4.2 Graph multi-equipartition

In this section, we consider the case where $k\geq 3$. In this case, the
spectral upper bound in (6) already implies that any feasible point $R$ is a
smooth point of $\mathcal{B}_{n,r}$, then we simply treat $\mathcal{B}_{n,r}$
as a manifold. We use an augmented Lagrangian method on $\mathcal{B}_{n,r}$ to
solve the problem SDPLR1 in (6). Note that recently the convergence of ALM on
manifold has been studied in [24, 40, 45]. We reformulate problem (6) as
follows:

$\min\left\\{\frac{1}{2}\big{\langle}L,\,RR^{\top}\big{\rangle}+\delta_{\mathbb{S}_{+}^{r}}\left(Y\right):\
R^{\top}R+Y=\alpha I,\ R\in\mathcal{B}_{n,r}\right\\}.$ (60)

The Lagrangian function is:

$\mathcal{L}\left(R,Y,Z\right):=\frac{1}{2}\big{\langle}L,\,RR^{\top}\big{\rangle}-\big{\langle}Z,\,R^{\top}R+Y-\alpha
I\big{\rangle}+\delta_{\mathbb{S}^{r}_{+}}\left(Y\right)+\delta_{\mathcal{B}_{n,r}}\left(R\right).$
(61)

Let $\Pi_{\mathbb{S}^{r}_{-}}(\cdot)$ be the projection operator onto
$\mathbb{S}^{r}_{-}$, the cone of $r\times r$ symmetric negative semidefinite
matrices. The reduced augmented Lagrangian function is:

$\mathcal{L}_{\beta}\left(R,Z\right)=\frac{1}{2}\big{\langle}L,\,RR^{\top}\big{\rangle}+\frac{1}{2\beta}\left(\|\Pi_{\mathbb{S}_{-}^{r}}\left(Z-\beta\left(R^{\top}R-\alpha
I\right)\right)\|^{2}-\|Z\|^{2}\right).$ (62)
$\nabla_{R}\mathcal{L}_{\beta}\left(R,Z\right)=LR-2R\,\Pi_{\mathbb{S}_{-}^{r}}\left(Z-\beta\left(R^{\top}R-\alpha
I\right)\right).$

The template of the augmented Lagrangian method is as follows:

Algorithm 2:

Choose $R_{0}\in\mathcal{B}_{n,r},$ $Z_{0}\in\mathbb{S}^{r},$ $\beta_{0}>0.$
Set $k=0.$

Step 1.

$R_{k}\in\arg\min\left\\{\mathcal{L}_{\beta_{k}}\left(R,Z_{k-1}\right):\
R\in\mathcal{B}_{n,r}\right\\}.$

Step 2.

$Z_{k}=\Pi_{\mathbb{S}_{-}^{r}}\left(Z_{k-1}-\beta_{k}\left(R_{k}^{\top}R_{k}-\alpha
I\right)\right).$

Step 3.

Choose $\beta_{k+1}$, set $k\leftarrow k+1.$ Goto Step 1.

######  Proposition 4.2.

Any saddle point of the Lagrangian function (61) satisfies the first order KKT
condition of (6).

###### Proof.

Suppose $\left(R,Y,Z\right)$ is a saddle point of (61). We have

$R\in\arg\min_{R^{\prime}}\left\\{\mathcal{L}\left(R^{\prime},Y,Z\right)\right\\},\
Y\in\arg\min_{Y^{\prime}}\left\\{\mathcal{L}\left(R,Y^{\prime},Z\right)\right\\},\
Z\in\arg\max_{Z^{\prime}}\left\\{\mathcal{L}\left(R,Y,Z^{\prime}\right)\right\\},$

which implies that

$\displaystyle{\rm Proj}_{R}\left(LR-2RZ\right)=0$ (63) $\displaystyle
Z\in\mathcal{N}_{\mathbb{S}^{r}_{+}}\left(Y\right)$ (64) $\displaystyle
R^{\top}R+Y-\alpha I=0.$ (65)

The condition (64) implies that $Y\in\mathbb{S}_{+}^{r},$
$Z\in\mathbb{S}_{-}^{r}$ and $\left\langle Z,Y\right\rangle=0.$ From (65),
$RR^{\top}-\alpha I\preceq 0$ and so $R\in\Omega_{\alpha}^{2}.$ Since
$\left\langle Z,Y\right\rangle=0$, we have that $\left\langle
Z,R^{\top}R-\alpha I\right\rangle=0.$ Because $Z,\ R^{\top}R-\alpha I\preceq
0$, the last complementarity condition implies that $Z(R^{\top}R-\alpha I)=0$
and hence $Z$ and $R^{\top}R-\alpha I$ are simultaneously diagonalizable. By
considering the SVD of $R$ as in (16), we have $Z=V_{1}Z_{1}V_{1}^{\top}$ for
some $Z_{1}\preceq 0.$ So
$RZ=U_{1}\sqrt{\alpha}Z_{1}V_{1}^{\top}\in-\mathcal{N}_{\Omega_{\alpha}^{2}}\left(R\right).$
The condition (63) implies $R\in\mathcal{B}_{n,r}$ and

$J\left(LR-2RZ-{\rm diag}\left(\lambda_{LR-2RZ}^{R}\right)R\right)=0,$

where $\lambda_{LR-2RZ}^{R}$ is defined as in (36). Thus,

$\left(L-{\rm
diag}\left(\lambda_{LR-2RZ}^{R}\right)\right)R-e\mu^{\top}=2RZ\in-\mathcal{N}_{\Omega_{\alpha}^{2}}\left(R\right),$

which implies (20), where $\mu=R^{\top}\left(L-{\rm
diag}\left(\lambda_{LR-2RZ}^{R}\right)\right)e/n$. Therefore,
$R,\lambda_{LR-2RZ},\mu$ satisfies the first order KKT condition of (6).

We terminate the augmented Lagrangian method based on the relative dual
feasibility,

$\max\Bigg{\\{}\frac{\|{\rm
Proj}_{R}\left(LR-2RZ\right)\|_{F}}{1+\|R\|_{F}+\|Z\|_{F}},\ \frac{\|Y-{\rm
Proj}_{\mathbb{S}_{+}^{r}}\left(Y+Z\right)\|_{F}}{1+\|Y\|_{F}+\|Z\|_{F}}\Bigg{\\}},$

and the relative primal feasibility, $\|RR^{\top}+Y-\alpha
I\|_{F}/(1+\|Y\|_{F}+\|R\|_{F}).$ We stop the algorithm if both the primal and
dual feasibility are less than some given tolerance ${\rm tol}$.

##  5 Numerical experiment

In this section, we conduct numerical experiments to verify the efficiency of
our methods. Note that our goal is not to show that the Riemannian
optimization methods in section 4 are efficient since they are already well
studied the literature. Our purpose is to show that by using the geometric
structure of the algebraic variety $\mathcal{B}_{n,r}$, we are able to solve
the graph equipartition SDP problem (6) much more efficiently than other
methods that do not exploit the geometric structure of $\mathcal{B}_{n,r}$ or
only partially exploiting geometric structure like oblique manifold. All the
experiments were run in Matlab R2020b on a MacBook Pro with 1.4 GHz Quad-Core
Intel Core i5 processor and 8GB RAM.

###  5.1 A rank adaptive strategy

In our implementation, we choose the parameter $r$ adaptively to further
increase computational efficiency. We first set
$r=k-1+\lceil\sqrt{2(n+1)}\rceil,$ which is the theoretical upper bound.
Suppose $R\in\mathbb{R}^{n\times r}$ is an iteration point with singular
values $\sigma_{1}\geq\sigma_{2}\ldots\geq\sigma_{r}>0.$ We consider the ratio
of them $\big{\\{}\frac{\sigma_{i}}{\sigma_{i+1}}:\ i\in[r-1]\big{\\}}.$ If
there exists $i\in[r-1]$ such that $\frac{\sigma_{i}}{\sigma_{i+1}}>10$, then
we choose $j=\arg\max\big{\\{}\frac{\sigma_{i}}{\sigma_{i+1}}:\
i\in[r-1]\big{\\}}$ and drop the singular values
$\sigma_{j+1},\sigma_{j+2},\ldots,\sigma_{r}$ and their corresponding singular
vectors to save storage. We perform such rank reduction check in every 10
steps. Note that such a rank adaptive strategy is adopted from [15].

###  5.2 Experiments on minimum bisection SDPs

First, we consider the SDPs coming from minimum bisection problems. In [24],
Boumal et al. used an ALM on Riemannian manifold141414Source codes from
https://github.com/losangle/Optimization-on-manifolds-with-extra-constraints.
called RALM to solve the minimum bisection SDP problems. Their manifold is
${\rm OB}_{n,r}$ and they penalized the extra constraint $R^{\top}e=0.$ We
also compare with the well-known BM-method151515Source codes from
https://sburer.github.io/projects.html.(see [9, 10]) called SDPLR, which
directly apply an ALM on the factorized model (2) without using the underlying
manifold structure. We name our method AVBB, since we use gradient descent
with BB step on an algebraic variety. For RALM, the authors stop the algorithm
when the distance of two consecutive iterations is less than $10^{-10}.$ For
SDPLR, the authors stop their algorithm when the primal residue is less than
some tolerance, which we choose to be $10^{-6}$. For our method, we choose
$\delta_{0}=0.02$ in Algorithm 1. We stop the algorithm in step 1 if the
relative Riemannian gradient norm is less than ${\rm tol}=10^{-6}.$ Once we
get the solution $R$, define

$S_{\lambda}:=L-\mathcal{A}^{*}\left(\lambda\right).$ (66)

We check the KKT residues of computed solution $X=RR^{\top}$ for the SDP
problem (5) as follows:

$\displaystyle{\rm primal\ residue:}\ {\rm
Rp}=\frac{\|\mathcal{A}\left(RR^{\top}\right)-b\|_{2}}{1+\|b\|_{2}},\quad{\rm
dual\ residue:}\ {\rm
Rd}=\frac{\max\Big{\\{}0,-\lambda_{\min}\big{(}JS_{\lambda_{LR}^{R}}J\big{)}\Big{\\}}}{1+\|L\|_{F}},$
$\displaystyle{\rm\ complementarity:}\ {\rm
Rc}=\frac{\big{|}\big{\langle}RR^{\top},\,S_{\lambda_{LR}^{R}}\big{\rangle}\big{|}}{1+\|L\|_{F}},$
(67)

where $\lambda_{LR}^{R}$ is defined as in (36).

Because RALM and SDPLR do not use rank adaptive strategy. We first use our
method AVBB to solve the SDP problem and check the KKT residues of the output
$\hat{R}\in\mathbb{R}^{n\times\tilde{r}}.$ If we are able to solve SDP problem
successfully, we choose $r=\tilde{r}+5$ for SDPLR and RALM. Note that by doing
this, we already give the other solvers some advantage since $r=\tilde{r}+5$
is usually much smaller than the theoretical upper bound. Also, we choose $r$
to be slightly large than $\tilde{r}$ so that they are less likely to reach a
spurious local minimum. For the initialization, we generate a random matrix,
randn(n,r), and compute its retraction onto $\mathcal{B}_{n,r}$ to get an
initial point $R_{0}$ for our algorithm. For RALM and SDPLR, we use their own
initialization. When computing the dual residue, we use PROPACKmod in [39]
that is modified from PROCPACK in [23] to compute the smallest eigenvalue,
since it can make use of sparsity and low-rank property.

The graph data set we use in experiments comes from some small dense graphs
from [35] and large social networks from [33].

Table 1: Comparison of AVBB, RALM and SDPLR for minimum bisection SDP. problem | algorithm | Rp | Rd | Rc | obj | time
---|---|---|---|---|---|---
brock200-1 | AVBB | 1.97e-16 | 1.28e-10 | 8.62e-16 | 3.9857926e+03 | 2.95e-01
n=200 | RALM | 1.28e-06 | 5.46e-10 | 1.51e-06 | 3.9857915e+03 | 1.64e+01
m=5067 | SDPLR | 1.32e-07 | 4.75e-10 | 4.93e-07 | 3.9857923e+03 | 3.58e+00
brock200-4 | AVBB | 1.72e-16 | 1.31e-10 | 1.02e-15 | 5.6562152e+03 | 1.56e-01
n=200 | RALM | 1.21e-07 | 2.17e-09 | 8.35e-07 | 5.6562144e+03 | 1.90e+01
m=6812 | SDPLR | 1.32e-07 | 6.44e-11 | 8.56e-07 | 5.6562143e+03 | 7.52e+00
brock400-1 | AVBB | 1.95e-16 | 4.84e-10 | 1.72e-15 | 1.6933377e+04 | 3.43e-01
n=400 | RALM | 1.50e-07 | 1.10e-08 | 3.72e-07 | 1.6933376e+04 | 2.55e+01
m=20078 | SDPLR | 9.51e-08 | 3.40e-10 | 2.97e-07 | 1.6933377e+04 | 3.04e+01
c-fat200-1 | AVBB | 1.57e-16 | 1.09e-10 | 1.18e-15 | 1.7825896e+04 | 1.07e-01
n=200 | RALM | 1.09e-06 | 7.39e-12 | 3.24e-07 | 1.7825895e+04 | 7.62e+00
m=18367 | SDPLR | 1.31e-07 | 9.08e-11 | 1.05e-07 | 1.7825896e+04 | 3.06e+00
hamming-6-4 | AVBB | 1.72e-16 | 1.56e-09 | 5.78e-16 | 1.0240000e+03 | 6.98e-02
n=64 | RALM | 1.81e-16 | 2.97e-10 | 7.15e-16 | 1.0240000e+03 | 3.41e-01
m=1313 | SDPLR | 1.08e-07 | 1.47e-06 | 2.40e-07 | 1.0240000e+03 | 3.20e-01
hamming-7-5-6 | AVBB | 2.16e-16 | 4.23e-09 | 1.12e-15 | 1.5360000e+03 | 5.51e-02
n=128 | RALM | 2.14e-16 | 1.26e-10 | 9.78e-16 | 1.5360000e+03 | 1.02e+00
m=1793 | SDPLR | 1.50e-07 | 2.24e-07 | 1.09e-07 | 1.5360000e+03 | 3.44e-01
hamming-8-3-4 | AVBB | 1.83e-16 | 2.16e-10 | 1.36e-15 | 1.4336000e+04 | 5.06e-02
n=256 | RALM | 2.41e-13 | 1.48e-10 | 2.92e-15 | 1.4336000e+04 | 2.08e+00
m=16129 | SDPLR | 3.18e-08 | 9.79e-08 | 1.80e-07 | 1.4336000e+04 | 4.20e-01
hamming-8-4 | AVBB | 2.01e-16 | 9.02e-11 | 8.05e-16 | 7.4240000e+03 | 6.57e-02
n=256 | RALM | 3.35e-16 | 1.67e-11 | 1.15e-15 | 7.4240000e+03 | 1.40e+00
m=11777 | SDPLR | 5.76e-08 | 1.95e-07 | 2.22e-08 | 7.4240000e+03 | 4.50e-01
hamming-9-5-6 | AVBB | 2.22e-16 | 1.81e-10 | 2.42e-15 | 5.0176000e+04 | 1.37e-01
n=512 | RALM | 2.69e-14 | 2.87e-11 | 3.12e-15 | 5.0176000e+04 | 5.41e+00
m=53761 | SDPLR | 2.63e-08 | 8.32e-08 | 4.17e-08 | 5.0176000e+04 | 7.59e-01
hamming-9-8 | AVBB | 0.00e+00 | -0.00e+00 | 0.00e+00 | 0.0000000e+00 | 3.56e-01
n=512 | RALM | 2.22e-16 | 4.01e-19 | 1.40e-15 | 1.1220542e-14 | 3.24e+00
m=2305 | SDPLR | 3.32e-08 | 0.00e+00 | 1.11e-10 | 1.1938166e-08 | 6.37e-01
hamming-10-2 | AVBB | 0.00e+00 | -0.00e+00 | 0.00e+00 | 0.0000000e+00 | 6.32e-01
n=1024 | RALM | 2.22e-16 | 0.00e+00 | 3.04e-15 | 1.8270875e-14 | 5.61e+00
m=23041 | SDPLR | 1.01e-08 | 5.07e-14 | 4.43e-11 | 3.2290904e-08 | 1.26e+00
hamming-11-2 | AVBB | 0.00e+00 | -0.00e+00 | 0.00e+00 | 0.0000000e+00 | 6.95e-01
n=2048 | RALM | 2.27e-16 | 7.55e-18 | 1.56e-14 | 2.7414197e-12 | 2.28e+01
m=56321 | SDPLR | 2.96e-08 | 1.56e-18 | 3.83e-11 | 4.8081223e-08 | 5.05e+00
hamming6-2 | AVBB | 1.28e-16 | 3.03e-09 | 5.47e-17 | 6.4000000e+01 | 4.40e-02
n=64 | RALM | 4.47e-16 | 1.31e-08 | 2.33e-16 | 6.4000000e+01 | 1.67e+00
m=193 | SDPLR | 1.85e-07 | 8.43e-08 | 2.39e-07 | 6.4000019e+01 | 4.33e-01
hamming8-2-G | AVBB | 1.96e-16 | 4.60e-09 | 3.91e-16 | 2.5600000e+02 | 6.54e-02
n=256 | RALM | 1.84e-16 | 2.00e-10 | 4.27e-16 | 2.5600000e+02 | 2.28e+00
m=1024 | SDPLR | 9.34e-08 | 1.02e-06 | 1.07e-08 | 2.5600000e+02 | 5.53e-01
MANN-a27 | AVBB | 2.36e-16 | 3.63e-09 | 2.49e-16 | 1.2219459e+02 | 2.95e-01
n=378 | RALM | 2.22e-16 | 1.73e-10 | 2.24e-16 | 1.2219459e+02 | 2.36e+01
m=703 | SDPLR | 8.08e-08 | 7.13e-07 | 1.13e-07 | 1.2219460e+02 | 7.51e-01
johnson8-4-4 | AVBB | 1.86e-16 | 1.18e-09 | 3.92e-16 | 2.8000000e+02 | 2.86e-02
n=70 | RALM | 1.73e-16 | 8.38e-10 | 4.85e-16 | 2.8000000e+02 | 6.10e-01
m=561 | SDPLR | 1.12e-07 | 1.55e-06 | 5.82e-08 | 2.8000000e+02 | 4.20e-01
johnson16-2-4 | AVBB | 1.98e-16 | 1.34e-11 | 4.06e-16 | 9.6000000e+02 | 4.86e-02
n=120 | RALM | 1.85e-16 | 1.01e-10 | 8.23e-16 | 9.6000000e+02 | 1.09e+00
m=1681 | SDPLR | 6.41e-08 | 7.84e-07 | 1.62e-07 | 9.5999998e+02 | 3.38e-01
keller4 | AVBB | 1.88e-16 | 1.42e-09 | 4.34e-16 | 3.2898566e+03 | 1.12e-01
n=171 | RALM | 7.57e-06 | 1.56e-11 | 1.08e-05 | 3.2898481e+03 | 4.66e+00
m=5101 | SDPLR | 1.42e-07 | -0.00e+00 | 1.47e-06 | 3.2898555e+03 | 7.43e+00
keller5-G | AVBB | 1.97e-16 | 4.69e-10 | 2.17e-15 | 4.7829110e+04 | 5.65e-01
n=776 | RALM | 1.76e-07 | 4.30e-10 | 3.51e-07 | 4.7829108e+04 | 3.16e+01
m=74710 | SDPLR | 6.92e-08 | 3.79e-11 | 2.20e-07 | 4.7829109e+04 | 4.45e+01
p-hat300-1 | AVBB | 1.75e-16 | 4.13e-11 | 1.49e-15 | 3.2077515e+04 | 1.71e-01
n=300 | RALM | 7.60e-07 | 1.00e-09 | 1.10e-06 | 3.2077511e+04 | 1.66e+01
m=33918 | SDPLR | 1.09e-07 | 6.84e-11 | 4.14e-07 | 3.2077514e+04 | 2.56e+01
san200-07-1 | AVBB | 1.84e-16 | 1.04e-09 | 1.19e-15 | 4.8782920e+03 | 1.84e-01
n=200 | RALM | 1.99e-07 | 3.22e-09 | 8.08e-07 | 4.8782913e+03 | 1.65e+01
m=5971 | SDPLR | 1.31e-07 | 2.78e-10 | 6.54e-07 | 4.8782914e+03 | 6.26e+00
musae-PTBR | AVBB | 1.67e-16 | 1.21e-10 | 1.09e-14 | 1.0504381e+04 | 4.64e+00
n=1912 | RALM | 5.45e-08 | 4.26e-06 | 5.48e-06 | 1.0504640e+04 | 6.85e+01
m=31299 | SDPLR | 5.98e-07 | 1.09e-10 | 1.84e-05 | 1.0504329e+04 | 4.70e+02
musae-chameleon | AVBB | 1.42e-16 | 1.28e-10 | 3.18e-15 | 4.7543522e+02 | 2.40e+01
n=2277 | RALM | 3.19e-09 | 2.04e-05 | 1.03e-07 | 4.7833683e+02 | 8.82e+01
m=36101 | SDPLR | 4.79e-08 | 4.77e-19 | 4.92e-07 | 4.7543371e+02 | 4.76e+02
musae-RU | AVBB | 1.59e-16 | 5.82e-11 | 9.31e-16 | 1.1774237e+04 | 1.49e+01
n=4385 | RALM | - | - | - | - | -
m=37304 | SDPLR | 1.97e-05 | 9.37e-09 | 5.47e-05 | 1.1774072e+04 | 4.77e+02
musae-ES | AVBB | 1.67e-16 | 1.55e-11 | 6.02e-16 | 1.5038307e+04 | 6.45e+00
n=4648 | RALM | - | - | - | - | -
m=59482 | SDPLR | 7.30e-05 | 7.19e-07 | 9.63e-05 | 1.5038055e+04 | 5.02e+02
musae-squirrel | AVBB | 1.66e-16 | 3.65e-10 | 8.89e-15 | 1.0907331e+04 | 3.77e+01
n=5201 | RALM | - | - | - | - | -
m=217073 | SDPLR | 3.32e-04 | 1.78e-07 | 9.43e-04 | 1.0899471e+04 | 5.06e+02
musae-FR | AVBB | 1.55e-16 | 2.56e-05 | 1.62e-15 | 3.6960466e+04 | 2.04e+01
n=6549 | RALM | - | - | - | - | -
m=112666 | SDPLR | 4.01e-02 | 4.41e-05 | 2.24e-01 | 3.4091560e+04 | 5.01e+02
musae-ENGB | AVBB | 1.71e-16 | 7.65e-10 | 1.06e-15 | 8.9985089e+03 | 1.04e+01
n=7126 | RALM | - | - | - | - | -
m=35324 | SDPLR | 9.42e-03 | 4.88e-06 | 8.45e-02 | 8.8164546e+03 | 4.96e+02
musae-DE | AVBB | 1.53e-16 | 4.69e-11 | 4.17e-16 | 4.9503018e+04 | 2.88e+01
n=9498 | RALM | - | - | - | - | -
m=153138 | SDPLR | - | - | - | - | -
musae-crocodile | AVBB | 1.83e-16 | 9.02e-12 | 1.70e-15 | 1.3265035e+04 | 3.32e+01
n=11631 | RALM | - | - | - | - | -
m=180020 | SDPLR | - | - | - | - | -
musae-facebook | AVBB | 1.96e-16 | 1.20e-09 | 3.17e-15 | 6.0495334e+03 | 1.12e+02
n=22470 | RALM | - | - | - | - | -
m=171002 | SDPLR | - | - | - | - | -

From Table 1, we can see that if the graph size is moderate, all the three
methods can find optimal solution for the SDP (5) successfully since its
corresponding KKT residues (Rp,Rd,Rc) are all smaller than the required
tolerance. However, when the graph size is too large, such as musae-DE, musae-
crocodile and musae-facebook, only AVBB can solve these problems. In
particular, we are able to solve the largest instance musae-facebook with
$n=22470$ and $m=171002$ in about two minutes.

Both RALM and SDPLR will terminate if there is little progress for many
iterations. In the tables, “-” means that they terminate prematurely and
return a solution that is far from the optimal solution. Among all the
problems, AVBB is clearly faster and more accurate than RALM and SDPLR. Note
that for hamming-9-8, hamming-10-2 and hamming-11-2, the optimal solution is a
singular point and AVBB finds them successfully.

###  5.3 Experiments on graph equipartition SDPs

In this section, we test on graph equipartition SDP problems with more than 2
partitions. Since RALM and SDPLR cannot handle the extra SDP upper bound
constraint $X\preceq\alpha I$, we do not test them in the experiments.
Instead, we compare our method with the interior point method (which we denote
as IPM) in section 6 of [38] that can handle the SDP upper bound constraint
directly. Note that applying the IPM in [38] to handle the semidefinite upper
bound constraint directly is much more efficient than applying an IPM to the
reformulated problem of converting the bound constraint to an affine
constraint with an additional slack SDP variable.

We call our method AVALM where “AV” stands for “algebraic variety”. We use a
Riemannian gradient method with BB step and nonmonotone linesearch to solve
the ALM subproblem to the accuracy of $10^{-6}$ in the relative norm of the
Riemannian gradient. We also set the maximum number of iterations for solving
the subproblem to be 200. We set the initial penalty parameter $\beta=0.1$ and
set

$\beta_{k+1}=\begin{cases}\max\\{\beta_{k}/1.2,0.1\\}&{\rm
pfeas<dfeas/1000}\\\ \min\\{1.2\beta_{k},10\\}&{\rm
pfeas\geq\max\\{dfeas/1000,10\cdot tol\\}}.\end{cases}$

We stop our ALM algorithm if both the primal feasibility (pfeas) and dual
feasibility (dfeas) are less than $10^{-6}.$ We also stop our algorithm if
there is no progress for many iterations. After we get a solution $(R,Z)$, let
${\bf S}=S_{\lambda_{LR-2RZ}}$ as in (66), then we check the following
relative KKT residues for the original convex SDP problem (5):

$\displaystyle{\rm primal\ residue:}\ {\rm
Rp}=\max\left\\{\frac{\|\mathcal{A}\left(RR^{\top}\right)-b\|_{2}}{1+\|b\|_{2}},\|R\|_{2}-\sqrt{\alpha}\right\\},$
$\displaystyle{\rm dual\ residue:}\ {\rm
Rd}=\frac{\max\left\\{\lambda_{\max}\left(U_{1}^{\top}{\bf S}U_{1}\right),\
\max\left\\{-\lambda_{\min}\left(J\left({\bf S}-U_{1}U_{1}^{\top}{\bf
S}U_{1}U_{1}^{\top}\right)J\right),0\right\\}\right\\}}{1+\|L\|_{F}},$
$\displaystyle{\rm complementarity:}\ {\rm Rc}=\frac{\Big{|}\left\langle{\bf
S}-U_{1}U_{1}^{\top}{\bf
S}U_{1}U_{1}^{\top},RR^{\top}\right\rangle\Big{|}}{1+\|L\|_{F}},$ (68)

where $U_{1}$ is defined in (16). If $\max\\{{\rm Rp},{\rm Rd},{\rm Rc}\\}$ is
small, then we have obtained an accurate approximate optimal solution to the
SDP problem (5).

Since the interior point method IPM cannot handle graph equipartition SDP
problems arising from large graphs, we only compare AVALM with IPM on small
datasets in [35], where we set the number of partitions $k=5.$ From Table
LABEL:TEQP1, we see that both AVALM and IPM can solve all of these problems
accurately. However, AVALM is much more efficient than IPM, and for some
instances AVALM can be 20-100 times faster.

Table 2: Comparison of AVALM with IPM for graph equipartition SDP with $k=5$. problem | algorithm | Rp | Rd | Rc | obj | time
---|---|---|---|---|---|---
brock200-1 | AVALM | 6.22e-06 | 9.53e-09 | 2.12e-08 | 3.9860693e+03 | 4.54e-01
n=200,m=5067 | IPM | 1.21e-07 | 2.37e-07 | 2.14e-06 | 3.9860704e+03 | 3.90e+00
brock200-4 | AVALM | 4.15e-06 | 4.47e-08 | 9.17e-11 | 5.6564739e+03 | 4.40e-01
n=200,m=6812 | IPM | 1.34e-07 | 5.56e-09 | 1.37e-06 | 5.6564750e+03 | 3.88e+00
brock400-1 | AVALM | 2.15e-16 | 1.12e-09 | 2.15e-15 | 1.6933377e+04 | 3.55e-01
n=400,m=20078 | IPM | 1.93e-07 | 8.08e-10 | 6.99e-06 | 1.6933383e+04 | 7.52e+00
c-fat200-1 | AVALM | 1.47e-06 | 8.34e-16 | 5.45e-18 | 1.7842168e+04 | 5.48e-01
n=200,m=18367 | IPM | 7.51e-08 | 1.22e-10 | 1.92e-06 | 1.7842169e+04 | 1.96e+00
hamming-6-4 | AVALM | 1.63e-16 | 1.01e-09 | 3.73e-16 | 1.0240000e+03 | 1.94e-02
n=64,m=1313 | IPM | 2.41e-09 | 1.63e-08 | 4.39e-08 | 1.0240002e+03 | 2.77e-01
hamming-7-5-6 | AVALM | 2.23e-16 | 2.99e-09 | 9.81e-16 | 1.5360000e+03 | 3.11e-02
n=128,m=1793 | IPM | 6.72e-09 | 1.12e-10 | 3.73e-07 | 1.5360001e+03 | 6.52e-01
hamming-8-3-4 | AVALM | 2.01e-16 | 6.37e-10 | 1.55e-15 | 1.4336000e+04 | 5.10e-02
n=256,m=16129 | IPM | 2.99e-08 | 4.08e-11 | 8.95e-08 | 1.4336000e+04 | 1.31e+00
hamming-8-4 | AVALM | 2.05e-16 | 6.60e-11 | 8.56e-16 | 7.4240000e+03 | 5.39e-02
n=256 | IPM | 1.74e-07 | 1.66e-10 | 1.61e-07 | 7.4240001e+03 | 1.32e+00
hamming-9-5-6 | AVALM | 2.55e-16 | 1.70e-10 | 2.23e-15 | 5.0176000e+04 | 1.12e-01
n=512,m=53761 | IPM | 8.70e-09 | -0.00e+00 | 2.15e-06 | 5.0176005e+04 | 3.74e+00
hamming-9-8 | AVALM | 2.38e-08 | 5.54e-09 | 1.41e-09 | 7.6800000e+02 | 3.08e-01
n=512,m=2305 | IPM | 3.50e-07 | 1.00e-13 | 5.65e-08 | 7.6800002e+02 | 5.78e+00
hamming-10-2 | AVALM | 7.25e-09 | 3.93e-10 | 8.14e-11 | 6.9120000e+03 | 4.19e-01
n=1024,m=23041 | IPM | 5.19e-08 | 0.00e+00 | 3.72e-07 | 6.9120002e+03 | 2.69e+01
hamming-11-2 | AVALM | 2.12e-16 | 3.04e-09 | 1.37e-08 | 1.5360000e+04 | 8.34e-01
n=2048,m=56321 | IPM | 4.05e-07 | 1.00e-12 | 2.87e-07 | 1.5360001e+04 | 2.10e+02
hamming6-2 | AVALM | 4.75e-06 | 4.50e-08 | 6.57e-08 | 6.4711622e+01 | 2.12e-01
n=64,m=193 | IPM | 5.12e-08 | 9.33e-09 | 9.55e-08 | 6.4711622e+01 | 6.46e-01
hamming8-2-G | AVALM | 1.89e-16 | 1.16e-08 | 4.04e-16 | 2.5600000e+02 | 5.53e-02
n=256,m=1024 | IPM | 1.63e-08 | 2.02e-10 | 2.76e-07 | 2.5600002e+02 | 1.27e+00
MANN-a27 | AVALM | 2.30e-16 | 2.58e-08 | 3.89e-16 | 1.2219459e+02 | 3.97e-01
n=378,m=703 | IPM | 4.14e-08 | -0.00e+00 | 9.37e-07 | 1.2219464e+02 | 3.05e+00
johnson8-4-4 | AVALM | 1.85e-16 | 1.16e-09 | 9.12e-09 | 2.8000000e+02 | 2.13e-02
n=70,m=561 | IPM | 4.01e-09 | 9.61e-10 | 1.12e-08 | 2.8000001e+02 | 2.05e-01
johnson16-2-4 | AVALM | 2.21e-16 | 5.68e-10 | 9.10e-16 | 9.6000000e+02 | 2.68e-02
n=120,m=1681 | IPM | 1.73e-08 | 2.92e-18 | 1.14e-06 | 9.6000017e+02 | 4.28e-01
keller4 | AVALM | 5.62e-06 | 2.77e-13 | 1.39e-08 | 3.2960645e+03 | 1.82e-01
n=171,m=5101 | IPM | 1.03e-07 | 1.39e-09 | 1.49e-06 | 3.2960652e+03 | 1.91e+00
keller5-G | AVALM | 1.94e-16 | 5.93e-10 | 2.26e-15 | 4.7829110e+04 | 9.28e-01
n=776,m=74710 | IPM | 3.37e-03 | 1.04e-13 | 5.75e-08 | 4.7855470e+04 | 3.56e+01
p-hat300-1 | AVALM | 4.07e-06 | 4.19e-10 | 5.46e-10 | 3.2077517e+04 | 2.83e-01
n=300,m=33918 | IPM | 3.43e-07 | 2.84e-09 | 3.57e-06 | 3.2077523e+04 | 5.83e+00
san200-07-1 | AVALM | 2.06e-06 | 8.32e-07 | 3.42e-09 | 4.8786431e+03 | 4.26e-01
n=200,m=5971 | IPM | 8.72e-09 | 8.93e-09 | 3.42e-06 | 4.8786445e+03 | 4.61e+00

We also test AVALM alone on graph equipartition SDPs arising large social
network graphs in [33], where we choose $k=3,5,10$ respectively. Table
LABEL:TEQP3 clearly shows that our algorithm can solve all of these large
scale SDP problems extremely efficiently and accurately with certified global
optimality for the SDP problem (5) based on ${\rm Rp}$, ${\rm Rd}$ and ${\rm
Rc}$. For the largest instance musae-facebook with $n=22470$ and $m=171002$,
we are able to solve the SDP problem in under 5 minutes when $k=3$ or $5$, and
in about 11 minutes when $k=10$.

Table 3: Test results of AVALM on large graph equipartition SDP. problem | k | Rp | Rd | Rc | obj | time
---|---|---|---|---|---|---
musae-PTBR | 3 | 1.26e-05 | 1.08e-09 | 3.88e-10 | 1.1097575e+04 | 2.84e+00
n=1912 | 5 | 3.53e-08 | 2.30e-10 | 3.67e-08 | 1.3089330e+04 | 6.41e+00
m=31299 | 10 | 2.09e-07 | 1.23e-08 | 1.20e-06 | 1.5773687e+04 | 1.25e+01
musae-chameleon | 3 | 9.60e-08 | 9.21e-10 | 2.61e-08 | 1.2425855e+03 | 3.94e+00
n=2277 | 5 | 2.07e-06 | 1.21e-09 | 5.51e-09 | 2.2463854e+03 | 3.77e+00
m=36101 | 10 | 3.17e-07 | 2.08e-07 | 7.58e-07 | 4.2365961e+03 | 2.85e+01
musae-RU | 3 | 1.69e-16 | 5.20e-11 | 2.36e-09 | 1.2051705e+04 | 7.31e+00
n=4385 | 5 | 1.18e-07 | 2.63e-08 | 8.56e-09 | 1.3456477e+04 | 1.07e+01
m=37304 | 10 | 1.75e-07 | 9.25e-10 | 1.41e-09 | 1.5615746e+04 | 1.43e+01
musae-ES | 3 | 6.99e-06 | 4.86e-10 | 2.68e-09 | 1.5267811e+04 | 4.89e+00
n=4648 | 5 | 8.00e-08 | 1.04e-09 | 3.05e-08 | 1.8335808e+04 | 5.57e+00
m=59482 | 10 | 1.25e-07 | 1.31e-08 | 6.58e-08 | 2.3154327e+04 | 3.18e+01
musae-squirrel | 3 | 1.61e-07 | 4.92e-08 | 1.63e-09 | 1.6163193e+04 | 2.71e+01
n=5201 | 5 | 3.27e-06 | 2.50e-09 | 5.44e-07 | 2.8359632e+04 | 5.00e+01
m=217073 | 10 | 4.14e-06 | 8.00e-08 | 2.39e-07 | 5.8591875e+04 | 2.54e+02
musae-FR | 3 | 4.52e-08 | 9.30e-09 | 8.35e-08 | 3.7199157e+04 | 1.38e+01
n=6549 | 5 | 8.64e-07 | 9.50e-11 | 4.35e-10 | 4.1936390e+04 | 9.01e+00
m=112666 | 10 | 7.28e-07 | 7.67e-10 | 2.52e-09 | 5.0585680e+04 | 3.55e+01
musae-ENGB | 3 | 1.90e-16 | 2.35e-09 | 7.57e-16 | 8.9985089e+03 | 8.38e+00
n=7126 | 5 | 3.56e-06 | 4.56e-08 | 3.70e-07 | 9.3864146e+03 | 1.67e+01
m=35324 | 10 | 1.04e-04 | 3.05e-08 | 1.70e-06 | 1.0792772e+04 | 2.85e+01
musae-DE | 3 | 1.95e-16 | 5.19e-10 | 1.99e-06 | 5.1971122e+04 | 2.26e+01
n=9498 | 5 | 6.75e-07 | 1.06e-10 | 1.88e-09 | 6.0468739e+04 | 1.73e+01
m=153138 | 10 | 1.58e-07 | 5.76e-10 | 3.24e-07 | 7.1324948e+04 | 1.06e+02
musae-crocodile | 3 | 5.28e-07 | 4.78e-11 | 5.25e-10 | 1.3265057e+04 | 3.11e+01
n=11631 | 5 | 2.73e-07 | 3.69e-10 | 1.04e-07 | 1.6310440e+04 | 5.30e+01
m=180020 | 10 | 3.56e-05 | 9.59e-10 | 2.30e-06 | 2.8584887e+04 | 2.32e+02
musae-facebook | 3 | 1.82e-16 | 8.27e-09 | 1.59e-15 | 6.0495349e+03 | 1.19e+02
n=22470 | 5 | 2.50e-06 | 5.50e-09 | 1.89e-08 | 6.7587115e+03 | 2.42e+02
m=171002 | 10 | 1.98e-16 | 1.62e-09 | 2.04e-09 | 9.6414771e+03 | 6.53e+02

##  6 Conclusion

In this paper, we study the SDP relaxation of a graph equipartition problem.
We study the optimality condition of the low-rank factorization model of the
SDP problem, which contains an additional SDP upper bound as compared to the
standard linear SDP problem. We prove the equivalence between the SDP problem
and its factorized version. In addition, we prove that under the constraint
nondegeneracy condition, any rank-deficient second order stationary point of
the factorized problem is a global optimal solution. To solve the SDP problem,
we study the properties of a special algebraic variety $\mathcal{B}_{n,r}$,
which comes from adding an extra coupling affine constraint to the oblique
manifold ${\rm OB}_{n,r}.$ We find a closed form solution of the projection
mapping onto $\mathcal{B}_{n,r}$, and more importantly, we prove that the
retraction is equivalent to a geometric median problem under some condition.
We also study the local geometric properties of singular points on
$\mathcal{B}_{n,r}.$ With these properties, we are able to use a gradient
descent method with BB step and nonmonotone line search and an augmented
Lagrangian method on the algebraic variety $\mathcal{B}_{n,r}$ to solve the
factorized SDP problem efficiently. The results from our numerical experiments
verify the high efficiency of our method. Our work shows that by making use of
the geometric property of the feasible region, one can design algorithms that
are significantly faster than other algorithms that do not rely on the
geometric property like SDPLR or only partially rely on its geometric property
like RALM.

## Acknowledgments

We thank the editors and reviewers for their valuable suggestions and
comments, which have helped to improve the quality of this paper.

## References

  * [1] A. I. Barvinok, Problems of distance geometry and convex properties of quadratic maps. Discrete Computational Geometry, 13 (1995): 189–202.
  * [2] J. Barzilai and J.M. Borwein, Two-point step size gradient methods. IMA J. Numerical Analysis, 8 (1988): 141–148.
  * [3] S. Benson, Y. Ye, and X. Zhang, Solving large-scale sparse semidefinite programs for combinatorial optimization. SIAM J. Optimization, 10 (2000): 443–461.
  * [4] S. Bhojanapalli, N. Boumal, P. Jain and P. Netrapalli, Smoothed analysis for low-rank solutions to semidefinite programs in quadratic penalty form. Conference On Learning Theory. PMLR, 2018.
  * [5] J.F. Bonnans and A. Shapiro, Perturbation analysis of optimization problems. Springer Science & Business Media, 2013.
  * [6] N. Boumal, V. Voroninski and A. S. Bandeira, The non-convex Burer-Monteiro approach works on smooth semidefinite programs. Advances in Neural Information Processing Systems, 29 (2016): 2757–2765.
  * [7] N. Boumal, V. Voroninski and A. S. Bandeira, Deterministic guarantees for Burer-Monteiro factorizations of smooth semidefinite programs. Communications on Pure and Applied Mathematics, 73 (2020): 581–608.
  * [8] N. Boumal, B. Mishra, P.-A. Absil and R. Sepulchre, Manopt, a matlab toolbox for optimization on manifolds. J. Machine Learning Research, 15 (2014): 1455–1459.
  * [9] S. Burer and R.D.C. Monteiro, A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization. Mathematical Programming, 95 (2003): 329–357.
  * [10] S. Burer and R.D.C. Monteiro, Local minima and convergence in low-rank semidefinite programming. Mathematical Programming, 103 (2005): 427–444.
  * [11] D. Cifuentes, On the Burer–Monteiro method for general semidefinite programs, Optimization Letters, 15 (2021): 2299–2309.
  * [12] D.A. Cox, J. Little and D. O’Shea, Ideals, varieties, and algorithms: an introduction to computational algebraic geometry and commutative algebra. Springer Science & Business Media, 2013.
  * [13] C. Ding, D. Sun and K-C. Toh, An introduction to a class of matrix cone programming. Mathematical Programming, 144 (2014): 141–179.
  * [14] S. Friedland and R. Loewy, Subspaces of symmetric matrices containing matrices with a multiple first eigenvalue. Pacific Journal of Mathematics, 62 (1976): 389-399.
  * [15] B. Gao and P.-A. Absil, A Riemannian rank adaptive method for low-rank matrix completion. arXiv preprint arXiv:2103.14768 (2021).
  * [16] M.-X. Goemans, Semidefinite programming in combinatorial optimization. Mathematical Programming, 79 (1997): 143–161.
  * [17] U. Görtz and T. Wedhorn, Algebraic Geometry I: Schemes. Vieweg+Teubner Verlag, (2010).
  * [18] S. Helmberg and F. Rendl, A spectral bundle method for semidefinite programming. SIAM J. Optimization, 10 (2000): 673–696.
  * [19] J. Hu, X. Liu, Z-W. Wen and Y-X. Yuan, A brief introduction to manifold optimization. J. Operations Research Society of China, 8 (2020): 199–248.
  * [20] D. Karger, R. Motwani and M. Sudan, Approximate graph coloring by semidefinite programming. J. ACM, 45 (1998): 246–265.
  * [21] S.E. Karisch and F. Rendl, Semidefinite programming and graph equipartition. In Topics in Semidefinite and Interior-Point Methods, American Mathematical Society, (1998): 77–95.
  * [22] B. Lannazzo and M. Porcelli, The Riemannian Barzilai-Borwein method with nonmonotone line search and the matrix geometric mean computation. IMA J. Numerical Analysis, 38 (2018): 495–517.
  * [23] R.M. Larsen, PROPACK – Software for large and sparse SVD calculations. Available at http://sun. stanford. edu/rmunk/PROPACK (2004): 2008–2009.
  * [24] C. Liu and N. Boumal, Simple algorithms for optimization on Riemannian manifolds with constraints. Applied Mathematics & Optimization, 82 (2020): 949–981.
  * [25] https://www.mosek.com/
  * [26] G. Pataki, On the rank of extreme matrices in semidefinite programs and the multiplicity of optimal eigenvalues. Mathematics of Operations Research, 23 (1998): 339–358.
  * [27] J. Peng and Y. Wei, Approximating k-means-type clustering via semidefinite programming. SIAM J. Optimization, 18 (2007): 186–205.
  * [28] K. B. Petersen and M. S. Pedersen, The matrix cookbook. Technical University of Denmark (2008): 7(15), 510.
  * [29] T. Pumir, S. Jelassi and N. Boumal, Smoothed analysis of the low-rank approach for smooth semidefinite programs. arXiv preprint arXiv:1806.03763 (2018).
  * [30] M. Raydan, The Barzilai and Borwein gradient method for the large scale unconstrained minimization problem. SIAM J. Optimization, 7 (1997): 26–33.
  * [31] F. Rendl, Semidefinite programming and combinatorial optimization. Applied Numerical Mathematics, 29 (1999): 255–281.
  * [32] S.M. Robinson, First order conditions for general nonlinear optimization. SIAM J. Applied Mathematics, 30 (1976): 597–607.
  * [33] B. Rozemberczki, C. Allen and R. Sarkar, Multi-scale Attributed Node Embedding. J. Complex Networks 9 (2021): cnab014.
  * [34] A. Shapiro, Sensitivity analysis of generalized equations. J. of Mathematical Sciences, 115 (2003).
  * [35] M. Trick, V. Chvatal, W. Cook, D. Johnson, C. McGeoch, and R. Tarjan, The Second DIMACS Implementation Challenge: NP Hard Problems: Maximum Clique, Graph Coloring, and Satisfiability. Rutgers University, New Brunswick, NJ, 1992; available at http://dimacs.rutgers.edu/Challenges/.
  * [36] K.-C. Toh, M. Todd and R. Tütüncü, SDPT3 – a MATLAB software package for semidefinite programming, version 1.3. Optimization Methods and Software, 11 (1999): 545–581.
  * [37] K.-C. Toh, M. Todd and R. Tütüncü, Solving semidefinite-quadratic-linear programs using SDPT3. Mathematical Programming, 95 (2003): 189–217.
  * [38] K.-C. Toh, R. Tütüncü and M. Todd, Inexact primal-dual path-following algorithms for a special class of convex quadratic SDP and related problems. Pacific J. Optimization, 3 (2007):135–164.
  * [39] K.C. Toh and S.W. Yun, An accelerated proximal gradient algorithm for nuclear norm regularized least squares problems. Pacific J. Optimization, 6 (2010): 615–640.
  * [40] Y. Wang, K. Deng, H. Liu and Z. Wen, A decomposition augmented Lagragian method for low-rank semidefinite programming. arXiv preprint arXiv:2109.11707 (2021).
  * [41] E. Weiszfeld, Sur le point lequel la Somme des distances de n points donnés est minimum. Tohoku Mathematical Journal, First Series 43 (1937): 355–386.
  * [42] L. Yang, D. Sun and K-C. Toh, SDPNAL+: a majorized semismooth Newton-CG augmented Lagrangian method for semidefinite programming with nonnegative constraints. Mathematical Programming Computation, 7 (2015): 331–366.
  * [43] L. Zhang, N. Zhang and X. Xiao, On the second-order directional derivatives of singular values of matrices and symmetric matrix-valued functions. Set-Valued and Variational Analysis, 21 (2013): 557–586.
  * [44] X-Y, Zhao, D. Sun and K-C. Toh, A Newton-CG augmented Lagrangian method for solving semidefinite programming. SIAM J. Optimization, 20 (2010): 1737–1765.
  * [45] Y. Zhou, C. Bao, C. Ding and J. Zhu, A semismooth Newton based augmented Lagragian method for nonsmooth optimization on matrix manifold. arXiv preprint arXiv:2103.02855 (2021).

##  Appendix A Useful auxiliary results

###  A.1 A corollary of the implicit function theorem

######  Lemma A.1.

Let $F:\mathbb{R}^{n}\times\mathbb{R}^{m}\rightarrow\mathbb{R}^{k}$ be a
continuously differentiable mapping in a neighbourhood of
$\left(x_{0},y_{0}\right)\in\mathbb{R}^{n}\times\mathbb{R}^{m}.$ Suppose
$F(x_{0},y_{0})=0$ and ${\rm
D}_{y}F(x_{0},y_{0}):\mathbb{R}^{m}\rightarrow\mathbb{R}^{k}$ is a surjective
linear mapping. Then there exists $\delta>0$ and a continuously differentiable
mapping $y:B_{\delta}(x_{0})\rightarrow\mathbb{R}^{k}$ such that
$y(x_{0})=y_{0}$ and $F(x,y(x))=0.$ Here,
$B_{\delta}(x_{0}):=\left\\{x\in\mathbb{R}^{n}:\ {\rm
dist}(x,x_{0})<\delta\right\\}.$

###### Proof.

Since ${\rm D}_{y}F(x_{0},y_{0})$ is surjective, we have that $m\geq k.$ If
$m=k,$ then Lemma A.1 directly follows from the implicit function theorem and
the mapping $y$ is unique. If $m>k$, we may find a linear mapping
$\mathcal{A}:\mathbb{R}^{m}\rightarrow\mathbb{R}^{m-k}$ such that the linear
operator $\triangle y\rightarrow\left({\rm D}_{y}F(x_{0},y_{0})[\triangle
y],\mathcal{A}(\triangle y)\right)$ is a bijective mapping from
$\mathbb{R}^{m}$ to $\mathbb{R}^{m}.$ Thus, we can define a new mapping
$\widehat{F}:\mathbb{R}^{n}\times\mathbb{R}^{m}\rightarrow\mathbb{R}^{m}$ such
that $\widehat{F}(x,y):=\left(F(x,y),\mathcal{A}(y)-\mathcal{A}(y_{0})\right)$
and apply the implicit function theorem to $\widehat{F}.$

###  A.2 Properties of the rounding procedure

######  Lemma A.2.

Suppose $n\geq 2$ and $\delta<\frac{1}{2}.$ Then for any
$R\in\mathcal{B}_{n,r}^{\delta+},$

$\left\|\texttt{Round}\left(R\right)\texttt{Round}\left(R\right)^{\top}-RR^{\top}\right\|\leq
2\sqrt{\delta}n.$ (69)

Moreover, if $\delta<\frac{1}{n^{2}},$ then
$\texttt{Round}(R)\in\mathcal{B}_{n,r}.$

###### Proof.

Let $s\in\mathbb{R}$ be the unique largest singular value of $R.$ Its
corresponding singular vectors are $u\in\mathbb{R}^{n}$ and
$v\in\mathbb{R}^{r}.$ Then we have that

$\|usv^{\top}-R\|^{2}=\|R\|^{2}-s^{2}=n-\|R\|_{2}^{2}\leq\delta n.$ (70)

Thus,

$\displaystyle\|\texttt{sgn}(u)-us\|=\|\texttt{sgn}(u)v^{\top}-usv^{\top}\|\leq{\rm
dist}\left(usv^{\top},{\rm OB}_{n,r}\right)$ $\displaystyle\leq{\rm
dist}\left(usv^{\top},\mathcal{B}_{n,r}\right)\leq\|usv^{\top}-R\|\leq\sqrt{\delta}\sqrt{n},$
(71)

where the first inequality comes from the retraction of oblique manifold, the
second inequality comes from $\mathcal{B}_{n,r}\subset{\rm OB}_{n,r},$ the
third inequality comes from $R\in\mathcal{B}_{n,r}.$ We have that,

$\displaystyle\left\|\texttt{Round}\left(R\right)\texttt{Round}\left(R\right)^{\top}-RR^{\top}\right\|=\left\|\texttt{sgn}(u)\texttt{sgn}(u)^{\top}-us^{2}u^{\top}\right\|$
$\displaystyle\leq\left\|\texttt{sgn}(u)\texttt{sgn}(u)^{\top}-\texttt{sgn}(u)su^{\top}\right\|+\left\|\texttt{sgn}(u)su^{\top}-us^{2}u^{\top}\right\|$
$\displaystyle=\left\|\texttt{sgn}(u)\right\|\left\|\left(\texttt{sgn}(u)^{\top}-su^{\top}\right)\right\|+\left\|\left(\texttt{sgn}(u)-us\right)\right\|\left\|su^{\top}\right\|$
$\displaystyle\leq 2\sqrt{n}.\sqrt{\delta}\sqrt{n}=2\sqrt{\delta}n.$ (72)

Now suppose that $\delta<\frac{1}{n^{2}}.$ From (70),
$\|usv^{\top}-R\|^{2}<1.$ Since every row of $R$ has unit length, we have that
every entry of $u$ is nonzero. Thus, $\texttt{sgn}(u)\in\\{-1,1\\}^{n}$ and
$\texttt{Round}(R)=\texttt{sgn}(u)e_{1}^{\top}\in{\rm OB}_{n,r}.$ From (A.2),
we have that

$\left|e^{\top}\texttt{sgn}(u)\right|=\left|e^{\top}\left(\texttt{sgn}(u)-us\right)\right|\leq\sqrt{n}\|\texttt{sgn}(u)-us\|\leq
n\sqrt{\delta}<1,$ (73)

where the first equality comes from that because $R^{\top}e=0$ and $u$ is an
singular vector of $R,$ $e^{\top}u=0.$ Because
$\left|e^{\top}\texttt{sgn}(u)\right|$ is an integer, we have that
$\left|e^{\top}\texttt{sgn}(u)\right|=0$ and so $\texttt{sgn}(u)\in E^{n}.$
Therefore,
$\texttt{Round}(R)=\texttt{sgn}(u)e_{1}^{\top}\in\mathcal{B}_{n,r}.$

##  Appendix B Escaping from a non-optimal singular point of
$\mathcal{B}_{n,r}$

In this section, we consider problem (5) with $\mathcal{A}={\rm diag}(.)$ and
its low rank version (7) without the spectral upper bound, i.e., the problems
(55) and (56) in subsection 4.1.

Consider some singular point $R=ae_{1}^{\top},$ where $a\in E^{n}.$ For any
$[0,H]\in\mathcal{T}_{\mathcal{B}_{n,r}}(R)$ and any $[-a\circ{\rm
diag}\left(HH^{\top}\right),W]\in\mathcal{T}^{2}_{\mathcal{B}_{n,r}}\left(R,[0,H]\right).$
Define

$\widehat{R}(t,H,W)=\left[a-\frac{t^{2}}{2}a\circ{\rm
diag}\big{(}HH^{\top}\big{)},\ tH+\frac{t^{2}}{2}W\right].$ (74)

From the property of the tangent cone and second order tangent set, we have
that

${\rm
Proj}_{\mathcal{B}_{n,r}}\big{(}\widehat{R}(t,H,W)\big{)}=\widehat{R}(t,H,W)+o(t^{2}).$
(75)

Substituting $\widehat{R}(t,H,W)$ into $f,$ we have that

$f\big{(}\widehat{R}(t,H,W)\big{)}=\frac{a^{\top}Ca}{2}-\frac{t^{2}a^{\top}C\left(a\circ{\rm
diag}\left(HH^{\top}\right)\right)}{2}+\frac{t^{2}\left\langle
CH,H\right\rangle}{2}+O(t^{3}).$ (76)

(75) implies that $W$ has little influence on the function value and we can
simply choose $W=0.$ Combining (75) and (76), we have that

$f\left({\rm
Proj}_{\mathcal{B}_{n,r}}\big{(}\widehat{R}(t,H,0)\big{)}\right)=f(R)+t^{2}F(H)+o(t^{2}),$
(77)

where

$F(H):=-\frac{a^{\top}C\left(a\circ{\rm
diag}\left(HH^{\top}\right)\right)}{2}+\frac{\left\langle
CH,H\right\rangle}{2}=\frac{1}{2}\left\langle\left(C-{\rm
Diag}\big{(}Caa^{\top}\big{)}\right)H,H\right\rangle.$ (78)

From (77), we know that in order to reduce the function value along a certain
direction as fast as possible, we should choose
$H\in\mathcal{T}_{\mathcal{B}_{n,r}}(R)$ such that the coefficient of $t^{2}$
in (77) is as small as possible. This introduces the following problem:

$\min\Bigg{\\{}F(H):\ e^{\top}H=0,\ a^{\top}{\rm
diag}\big{(}HH^{\top}\big{)}=0,\ \|H\|^{2}=n,\
H\in\mathbb{R}^{n\times(r-1)}\Bigg{\\}}.$ (79)

The constraint ensures that $[0,H]\in\mathcal{T}_{\mathcal{B}_{n,r}}(R).$
Also, we additionally fix the norm of $H$ to be $n$ because we only care about
the direction rather than the length of the vector. The following lemma says
that any feasible solution with a negative value provides us with a descent
direction that can escape from a singular point.

######  Lemma B.1.

Suppose $r>1,$ $R=ae_{1}^{\top}\in\mathbb{R}^{n\times r}$ is a singular point
of $\mathcal{B}_{n,r}$ for some $a\in E^{n}$ and $f$ is defined as in (56). If
problem (79) has some feasible solution $H\in\mathbb{R}^{n\times(r-1)}$ such
that its objective function value is negative. Then there exists some
$\alpha<0$ such that

$f\left({\rm
Proj}_{\mathcal{B}_{n,r}}\big{(}\widehat{R}(t,H,0)\big{)}\right)=f(R)+\alpha
t^{2}+o(t^{2}).$ (80)

Moreover, ${\rm Proj}_{\mathcal{B}_{n,r}}\big{(}\widehat{R}(t,H,0)\big{)}$ is
a smooth point of $\mathcal{B}_{n,r}$ for any nonzero $t$ that is sufficiently
small.

###### Proof.

For any feasible solution $H$ of problem (79), the constraints of (79) and
(34) implies that $[0,H]\in\mathcal{T}_{\mathcal{B}_{n,r}}(R).$ Therefore,
from (77) and (78), we have that

$f\left({\rm
Proj}_{\mathcal{B}_{n,r}}\big{(}\widehat{R}(t,H,0)\big{)}\right)=f(R)+F(H)t^{2}+o(t^{2}).$
(81)

Thus, if $F(H)<0$, then (80) holds by defining $\alpha:=F(H).$

Now we move on to prove the smoothness of ${\rm
Proj}_{\mathcal{B}_{n,r}}\big{(}\widehat{R}(t,H,0)\big{)}$ on the condition
that $F(H)<0.$ First, since $F(H)<0,$ from (80), we have that there exists
$\epsilon>0$ such that $f\big{(}{\rm
Proj}_{\mathcal{B}_{n,r}}\big{(}\widehat{R}(t,H,0)\big{)}\big{)}<f(R)$ for any
$0<|t|<\epsilon.$ Note that for any singular point $R_{1}$ of
$\mathcal{B}_{n,r},$ if $R_{1}=ab^{\top}$ for some $b\in{\rm S}^{r-1},$ then
we have that $f(R_{1})=f(R).$ This means that ${\rm
Proj}_{\mathcal{B}_{n,r}}\big{(}\widehat{R}(t,H,0)\big{)}$ can not be $R_{1}$
for any $0<|t|<\epsilon.$ On the other hand, if $R_{1}\notin a({\rm
S^{r-1}})^{\top}$, then we have that $R_{1}=\hat{a}b^{\top}$ for some
$\hat{a}\in E^{n}$ and $b\in{\rm S}^{r-1}$ such that $|a^{\top}\hat{a}|\leq
n-1.$ Hence

$\displaystyle\|R-R_{1}\|_{F}^{2}=\|R\|^{2}+\|R_{1}\|^{2}-2\left\langle
R,R_{1}\right\rangle=2n-2\left\langle
ae_{1}^{\top},\hat{a}b^{\top}\right\rangle=2n-2a^{\top}\hat{a}e_{1}^{\top}b$
$\displaystyle\geq 2n-2|a^{\top}\hat{a}|\geq 2n-2(n-1)=2.$ (82)

From (74) and (75), we have that ${\rm
Proj}_{\mathcal{B}_{n,r}}\big{(}\widehat{R}(t,H,0)\big{)}=R+O(t).$ This
together with (82) implies that there exists $\epsilon_{1}>0$ such that ${\rm
Proj}_{\mathcal{B}_{n,r}}\big{(}\widehat{R}(t,H,0)\big{)}$ can not be $R_{1}$
for any $|t|<\epsilon_{1}.$ Thus, ${\rm
Proj}_{\mathcal{B}_{n,r}}\big{(}\widehat{R}(t,H,0)\big{)}$ is smooth for any
$0<|t|<\min\left\\{\epsilon_{1},\epsilon_{2}\right\\}.$

Lemma B.1 tells us that we can find a descent direction at a singular point by
solving (79). With Lemma B.1, we make the following definition.

######  Definition B.2.

Suppose $r>1,$ $R=ae_{1}^{\top}\in\mathbb{R}^{n\times r}$ is a singular point
of $\mathcal{B}_{n,r}$ for some $a\in E^{n}$ and $f$ is defined as in (56).
$H\in\mathbb{R}^{n\times(r-1)}$ is called an escaping direction of
$(\ref{BCsmo})$ at $R=ae_{1}^{\top}$ if it is a feasible solution of (79) with
a negative objective value.

Note that (79) is a non-convex problem which may have spurious local minima.
Moreover, if there does not exist a feasible solution with a negative
objective function value for (79), then what can we tell about $R$? In order
to handle these two problems, we consider the following SDP relaxation of
(79):

$\min\Bigg{\\{}\left\langle\left(C-{\rm
Diag}\big{(}Caa^{\top}\big{)}\right),X\right\rangle:\begin{array}[]{l}Xe=0,\
\left\langle X,ee^{\top}\right\rangle=0,\ \left\langle{\rm
diag}(a),X\right\rangle=0,\\\ \left\langle
I,X\right\rangle=n,X\in\mathbb{S}^{n}_{+}\end{array}\Bigg{\\}}.$ (83)

Note that we add the redundant constraint $Xe=0$ similarly as in problem (5).
Let $J=PP^{\top},$ where $P\in{\rm St}(n,n-1).$ For any feasible solution $X$
of (83), we have that $X=JXJ=PP^{\top}XPP^{\top}.$ By introducing
$Y:=P^{\top}XP,$ problem (83) is equivalent to the following SDP problem with
a smaller size.

$\min\left\\{\left\langle C_{P},Y\right\rangle:\ \left\langle P^{\top}{\rm
diag}(a)P,Y\right\rangle=0,\ \left\langle I,Y\right\rangle=n,\
Y\in\mathbb{S}^{n-1}_{+}\right\\},$ (84)

where $C_{P}:=P^{\top}\left(C-{\rm
Diag}\left(Caa^{\top}\right)\right)P\in\mathbb{S}^{n-1}.$ The Lagrangian dual
problem of (84) is as follows:

$\max\left\\{ny_{1}:\ C_{P}-y_{1}I-y_{2}P^{\top}{\rm
diag}(a)P\in\mathbb{S}^{n-1}_{+},\ y_{1},y_{2}\in\mathbb{R}\right\\}.$ (85)

Note that we have

$\left\langle P^{\top}{\rm
diag}\left(a\right)P,I\right\rangle=\left\langle{\rm
diag}(a),J\right\rangle=\left\langle{\rm diag}\left(a\right),{\rm
diag}\left(J\right)\right\rangle=\frac{(n-1)}{n}\left\langle
a,e\right\rangle=0.$ (86)

(86) implies that $Y=\frac{nI}{n-1}$ is a strict feasible solution of (84).
Also, we may choose $y_{1}$ to be sufficiently negative to make
$\left(y_{1},y_{2}\right)$ strictly feasible for (85). The Slater’s conditions
for primal and dual problem imply that the duality gap is zero and the optimal
solution sets are non-empty for (84) and (85). Thus, problem (84) has a KKT
solution. It is easy to check that $P^{\top}aa^{\top}P$ is a feasible solution
of $(\ref{singSDP1})$ with optimal value $0.$ Thus, the optimal value of (84)
is non-positive. The following theorem tells that we can check the global
optimality of a singular point and escape from a non-optimal singular point of
(56) by solving the SDP problem (84).

######  Theorem B.3.

Suppose $r>1,$ $R=ae_{1}^{\top}\in\mathbb{R}^{n\times r}$ is a singular point
of $\mathcal{B}_{n,r}$ for some $a\in E^{n}.$ Let
$\left(Y,\lambda_{1},\lambda_{2}\right)$ be a KKT solution of (84).

  * (i)

If $\left\langle C_{P},Y\right\rangle=0,$ then $aa^{\top}$ is a global optimal
solution of (55).

  * (ii)

If $\left\langle C_{P},Y\right\rangle<0,$ and $Y=HH^{\top}$ for some
$H\in\mathbb{R}^{(n-1)\times(r-1)},$ then $PH$ is an escaping direction of
(56) at $R=ae_{1}^{\top}.$

###### Proof.

We first prove (i). $\left\langle C_{P},Y\right\rangle=0$ implies that the
optimal value of (84) is zero. Because $P^{\top}aa^{\top}P$ is also a feasible
solution of $(\ref{singSDP1})$ with optimal value $0,$
$\left(P^{\top}aa^{\top}P,\lambda_{1},\lambda_{2}\right)$ is a KKT solution of
(84). We have the following KKT conditions hold for
$\left(P^{\top}aa^{\top}P,\lambda_{1},\lambda_{2}\right):$

$\displaystyle\left\langle P^{\top}{\rm
diag}(a)P,P^{\top}aa^{\top}P\right\rangle=0,\ \left\langle
I,P^{\top}aa^{\top}P\right\rangle=n,\
P^{\top}aa^{\top}P\in\mathbb{S}^{n-1}_{+}.$ (87) $\displaystyle
C_{P}-y_{1}I-y_{2}P^{\top}{\rm diag}(a)P\in\mathbb{S}^{n-1}_{+}.$ (88)
$\displaystyle\left\langle C_{P}-y_{1}I-y_{2}P^{\top}{\rm
diag}(a)P,P^{\top}aa^{\top}P\right\rangle=0.$ (89)

Note that (87) holds for any $a\in E^{n}.$ Because $C_{P}=P^{\top}\left(C-{\rm
Diag}\left(Caa^{\top}\right)\right)P,$ (88) implies that:

$P\left(P^{\top}\left(C-{\rm
Diag}\big{(}Caa^{\top}\big{)}\right)P-y_{1}I-y_{2}P^{\top}{\rm
diag}(a)P\right)P^{\top}\in\mathbb{S}^{n}_{+}.$ (90)

After using $PP^{\top}=J$ and $J^{2}=J$ in (90), we get

$J\left(C-{\rm Diag}\big{(}Caa^{\top}\big{)}-y_{1}I-y_{2}{\rm
diag}(a)\right)J\in\mathbb{S}^{n}_{+},$ (91)

and (89) implies that

$\left\langle J\left(C-{\rm Diag}\big{(}Caa^{\top}\big{)}-y_{1}I-y_{2}{\rm
diag}(a)\right)J,aa^{\top}\right\rangle=0.$ (92)

Define $\mu:={\rm diag}\left(Caa^{\top}\right)-y_{1}e-y_{2}a.$ Combining (91),
(92) with feasibility of $aa^{\top}$, we have that:

$\displaystyle{\rm diag}\big{(}aa^{\top}\big{)}=e,\ aa^{\top}e=0,\
\left\langle aa^{\top},ee^{\top}\right\rangle=0,\
aa^{\top}\in\mathbb{S}^{n}_{+}.$ (93) $\displaystyle J\left(C-{\rm
diag}(\mu)\right)J\in\mathbb{S}^{n}_{+}.$ (94) $\displaystyle\left\langle
J\left(C-{\rm diag}(\mu)\right)J,aa^{\top}\right\rangle=0.$ (95)

The above equations are exactly the KKT condition of (55). Therefore,
$aa^{\top}$ is an optimal solution of (55).

Now we move on to prove (ii). Because $Y$ is feasible for (84), $PYP^{\top}$
is feasible for (83). Also, because $\left\langle C_{P},Y\right\rangle<0,$ we
have that $\left\langle C-{\rm
Diag}\big{(}Caa^{\top}\big{)},\,PYP^{\top}\right\rangle<0.$ Therefore, $PH$ is
a feasible solution of (79) with a negative function value. From Definition
B.2, $PH$ is an escaping direction of (56) at $R=ae_{1}^{\top}.$

Theorem B.3 says that when we come to a singular point $ae_{1}^{\top}$ for
some $a\in E^{n}.$ We can first solve another SDP problem (84), which always
has a KKT solution. If the optimal value for (84) is zero, then $aa^{\top}$ is
already an optimal solution of (55). Moreover, from the proof of Theorem B.3,
we can get the optimal dual variable of (55) directly from the optimal dual
variable of (84). If the optimal value for (84) is negative and the optimal
solution has rank $\leq r-1$, then we can construct an escaping direction from
the optimal solution of (84).

One may wonder whether it is too expensive to escape from a singular point by
solving an SDP with the matrix variable of size $(n-1)\times(n-1).$ Actually,
problem (84) is simple because there are only 2 constraints and there exists
an optimal solution of rank $\sqrt{2\cdot 2}=2.$ This implies that we can use
the Burer and Monteiro factorization to solve (84) efficiently. Moreover, we
can terminate the algorithm in advance when its primal feasibility is small
enough and its function value is negative because an escaping direction does
not require that $H$ to be a minimizer of (79).

Next we discuss how to compute the retraction ${\rm
Proj}_{\mathcal{B}_{n,r}}\big{(}\widehat{R}(t,H,0)\big{)}$ for a singular
point $R=ae_{1}^{\top}$ such that $a\in E^{n}.$ In section 3.3, we have only
discussed the retraction around a smooth point. When it comes to a singular
point, suppose $H\in\mathbb{R}^{n\times(r-1)}$ is an escaping direction of
(56). For the maximum eigenvalue function $\lambda_{\max}(\cdot),$ since $n$
is the only nonzero eigenvalue of $RR^{\top}=aa^{\top},$ $\lambda_{\max}$ is
differentiable at $RR^{\top}$ and the gradient is given by:

$\nabla\lambda_{\max}\big{(}RR^{\top}\big{)}=\frac{1}{n}aa^{\top}.$ (96)

The following proposition tells that the geometric median approach in
subsection 3.3 also works for a singular point.

######  Proposition B.4.

Suppose $r>1,$ $R=ae_{1}^{\top}\in\mathbb{R}^{n\times r}$ is a singular point
of $\mathcal{B}_{n,r}$ for some $a\in E^{n}$, and
$H\in\mathbb{R}^{n\times(r-1)}$ is an escaping direction of $(\ref{BCsmo})$ at
$R=ae_{1}^{\top}.$ There exists $\delta>0$ such that for any $0<|t|<\delta,$
there exists $0<\beta<1$ such that $\|\widehat{R}(t,H,0)\|<(1-\beta)\sqrt{n}$
and ${\rm
dist}\big{(}\widehat{R}(t,H,0),\mathcal{B}_{n,r}\big{)}<\frac{\beta\sqrt{n}}{\sqrt{n}+1}.$

###### Proof.

From (74),

$\displaystyle\widehat{R}(t,H,0)\widehat{R}(t,H,0)^{\top}=$
$\displaystyle\;\;aa^{\top}-\frac{t^{2}}{2}\left(a\,{\rm
diag}\big{(}HH^{\top}\big{)}^{\top}{\rm diag}(a)+{\rm diag}(a){\rm
diag}\big{(}HH^{\top}\big{)}a^{\top}\right)$
$\displaystyle+t^{2}HH^{\top}+O(t^{3}).$ (97)

Thus, from (96), we have that

$\lambda_{\max}\big{(}\widehat{R}(t,H,0)\widehat{R}(t,H,0)^{\top}\big{)}=n-t^{2}\big{(}n\|H\|^{2}-\|a^{\top}H\|^{2}\big{)}+O(t^{3}),$
(98)

where we have used the fact that $a\in E^{n}$ and so $a^{\top}a=n,$ ${\rm
diag}(a)a=e.$ Using the Cauchy-Schwarz inequality on every entry of
$a^{\top}H$, we have that

$\|a^{\top}H\|^{2}\leq\|a\|^{2}\|H\|^{2}=n\|H\|^{2}.$ (99)

The equality is attained in (99) if and only if $H=ab^{\top}$ for some vector
$b\in\mathbb{R}^{r-1}.$ However, this implies that the objective function
value of problem (79) is $0,$ which contradicts to that $H$ is an escaping
direction. Thus, we have that $\|a^{\top}H\|^{2}<n\|H\|^{2}.$ From the (98)
and relation between eigenvalue and singular value, we have that

$\|\widehat{R}(t,H,0)\|_{2}=\sqrt{n}\left(1-\frac{t^{2}}{2n}\left(n\|H\|^{2}-\|a^{\top}H\|^{2}\right)\right)+O(t^{3}).$
(100)

Also, from (75), we have that

${\rm dist}\big{(}\widehat{R}(t,H,0),\mathcal{B}_{n,r}\big{)}=o(t^{2}).$ (101)