# An Overview and Comparison of Spectral Bundle Methods for Primal and Dual
Semidefinite Programs††thanks: This work is supported by NSF ECCS-2154650.
Corresponding author: Yang Zheng (zhengy@eng.ucsd.edu).

Feng-Yi Liao  Department of Electrical and Computer Engineering, University of
California, San Diego<EMAIL_ADDRESS><EMAIL_ADDRESS> Lijun Ding
Wisconsin Institute for Discovery, University of Wisconsin–Madison, Madison
<EMAIL_ADDRESS> Yang Zheng  Department of Electrical and Computer
Engineering, University of California, San Diego<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The spectral bundle method developed by Helmberg and Rendl is well-established
for solving large-scale semidefinite programs (SDPs) in the dual form,
especially when the SDPs admit low-rank primal solutions. Under mild
regularity conditions, a recent result by Ding and Grimmer has established
fast linear convergence rates when the bundle method captures the rank of
primal solutions. In this paper, we present an overview and comparison of
spectral bundle methods for solving both primal and dual SDPs. In particular,
we introduce a new family of spectral bundle methods for solving SDPs in the
primal form. The algorithm developments are parallel to those by Helmberg and
Rendl, mirroring the elegant duality between primal and dual SDPs. The new
family of spectral bundle methods also achieves linear convergence rates for
primal feasibility, dual feasibility, and duality gap when the algorithm
captures the rank of the dual solutions. Therefore, the original spectral
bundle method by Helmberg and Rendl is well-suited for SDPs with low-rank
primal solutions, while on the other hand, our new spectral bundle method
works well for SDPs with low-rank dual solutions. These theoretical findings
are supported by a range of large-scale numerical experiments. Finally, we
demonstrate that our new spectral bundle method achieves state-of-the-art
efficiency and scalability for solving polynomial optimization compared to a
set of baseline solvers SDPT3, MOSEK, CDCS, and SDPNAL+.

##  1 Introduction

Semidefinite programs (SDPs) are an important class of convex optimization
problems that minimize a linear function in the space of positive semidefinite
(PSD) matrices subject to linear equality constraints [1]. Mathematically, the
standard primal and dual SDPs are in the form of

$\displaystyle\min_{X}$ $\displaystyle\langle C,X\rangle$ (P)
$\displaystyle\mathrm{subject\leavevmode\nobreak\ to}$ $\displaystyle\langle
A_{i},X\rangle=b_{i},\quad i=1,\ldots,m,$ $\displaystyle
X\in\mathbb{S}^{n}_{+},$

and

$\displaystyle\max_{y,Z}$ $\displaystyle b^{\mathsf{T}}y$ (D)
$\displaystyle\mathrm{subject\leavevmode\nobreak\ to}$ $\displaystyle
Z+\sum_{i=1}^{m}A_{i}y_{i}=C,$ $\displaystyle Z\in\mathbb{S}^{n}_{+},$

where $b\in\mathbb{R}^{m},C,A_{1},\ldots,A_{m}\in\mathbb{S}^{n}$ are the
problem data, $\mathbb{S}^{n}_{+}$ denotes the set of $n\times n$ PSD matrices
(we also write $X\succeq 0$ to denote $X\in\mathbb{S}^{n}_{+}$ when the
dimension is clear from the context or not important), and
$\langle\cdot,\cdot\rangle$ denotes the standard trace inner product on the
space of symmetric matrices.

SDPs offer a powerful mathematical framework that has gained significant
attention for decades [1, 2] and still receives strong research interests
today [3, 4, 5]. Indeed, SDP provides a versatile and robust modeling and
optimization approach for solving a wide range of problems in different
fields. Undoubtedly, SDPs have become powerful tools in control theory [6],
combinatorial optimization [7], polynomial optimization [2], machine learning
[8], and beyond [9]. In theory, one can solve any SDP instance up to arbitrary
precision in polynomial time using second-order interior point methods
(IPMs)[1, 10]. At each iteration of second-order IPMs, one usually needs to
solve a linear system with a coefficient matrix (i.e., the Schur complement
matrix) being generally dense and ill-conditioned. Consequently, IPMs often
suffer from both computational and memory issues when solving SDPs from large-
scale practical applications.

Improving the scalability of SDPs has gained significant attention in recent
years; see [5, 3, 11] for surveys. In particular, first-order methods (FOMs)
are at the forefront of developing scalable algorithms for solving large-scale
SDPs thanks to their low complexity per iteration. For instance, the
alternating direction method of multipliers (ADMM) is used to solve large-
scale SDPs in the dual form D [12]. The ADMM framework has been extended to
solve the homogenous self-dual embedding of SDPs P and D in [13]. In [14],
ADMM has been applied to solving SDPs with a quadratic cost function. It is
known that augmented Lagrangian methods (ALMs) are also suitable for solving
large-scale optimization problems. Some efficient ALM-based algorithms have
recently been developed to solve large-scale SDPs. For example, a Newton-CG
augmented Lagrangian method is proposed to solve SDPs with a large number of
affine constraints. An enhanced version is developed in [15] to further tackle
degenerate SDPs by employing a semi-smooth Newton-CG scheme coupled with a
warm start strategy [16]. The algorithms [15, 16] have been implemented in a
MATLAB package, SDPNAL+, which has shown promising numerical performance. To
tackle the storage issue, the sketching idea, approximating a large matrix $X$
without explicitly forming it, is exploited in the ALM framework together with
a conditional gradient method for SDPs [17]. In [18], an optimal storage
scheme is developed to solve SDPs by using a first-order method to solve the
dual SDP D and recovering the primal solution in P. Finally, a class of
efficient first-order spectral bundle methods has been developed to solve an
equivalent eigenvalue problem when primal SDPs enjoy a constant trace property
[19, 20].

Another important idea in designing efficient algorithms is to exploit the
underlying sparsity and structures in SDPs [3, 4, 21]. When SDPs have an
aggregate sparsity pattern, chordal decomposition [4, 3] has been exploited to
reduce the dimension of PSD constraints in the design of both IPMs [22] and
ADMM [23]. In [24], partially separable properties in conic programs
(including chordal decomposition) have been investigated to design efficient
first-order algorithms. On the other hand, when the SDPs have low-rank
solutions, low-rank factorization decomposing a big PSD matrix
$X\in\mathbb{S}^{n}_{+}$ into $VV^{\mathsf{T}}$, where
$V\in\mathbb{R}^{n\times r}$ with $r\ll n$, has been utilized reduce the
searching space in [25]. This low-rank factorization leads to a nonconvex
optimization problem, and there are significant efforts in advancing
theoretical understanding of the factorization approach [26, 27, 28]. Similar
to the low-rank factorization, one can approximate the PSD constraint
$X\in\mathbb{S}^{n}_{+}$ by $X=VSV^{\mathsf{T}}$ where
$S\in\mathbb{S}^{r}_{+}$ and the factor $V\in\mathbb{R}^{n\times r}$ is fixed.
This is one of the main ideas in the design of the spectral bundle method [19,
29] and the spectral Frank Wolf algorithm [30]. The core strategy in [19, 29,
30] is to iteratively search for the factor $V$ such that it spans the range
space of the optimal solution in P. Another similar approximation strategy is
the basis pursuit techniques in [31, 32].

In this paper, we focus on the development of the low-rank approximation
$X=VSV^{\mathsf{T}}$ in spectral bundle methods [19, 20, 33, 29]. The spectral
bundle method, originally developed by Helmberg and Rendl in [19], is well-
established to solve large-scale SDPs, thanks to its low per-iteration
complexity and fast practical convergence. Further developments of spectral
bundle methods appear in [20, 33]. Very recently, Ding and Grimmer established
sublinear convergence rates of the spectral bundle method in terms of primal
feasibility, dual feasibility, and duality gap, and further proved a linear
convergence rate when the algorithm captures a rank condition [29]. To the
best of our knowledge, all existing spectral bundle methods [29, 19, 34, 20,
33] focus on solving dual SDPs in the form of D. As shown in [29], these
spectral bundle methods are more desirable when the primal SDP P admits low-
rank solutions in which it is easier to enforce the rank condition to
guarantee linear convergence. On the other hand, when the dual SDP D admits
low-rank solutions, the existing spectral bundle methods may offer less
benefit in terms of convergence and efficiency. Indeed, SDPs arising from
moment/sum-of-squares (SOS) optimization problems and their applications [35,
36] are likely to admit low-rank solutions in the dual SDP (D) (this low-rank
property is consistent with the flat extension theory on the moment side [36,
Theorem 3.7] when it is formulated as a dual SDP).

In this work, we present an overview and comparison of spectral bundle methods
for solving both primal and dual SDPs. In particular, we introduce a new
family of spectral bundle methods for solving SDPs in the primal form (P). Our
algorithm developments are parallel to those by Helmberg and Rendl [19] which
focuses on solving dual SDPs (D), mirroring the elegant duality between primal
and dual SDPs. In particular, our contributions are as follows.

  * •

We propose a new family of spectral bundle methods, called
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP, for solving primal SDPs P, while all
existing methods [29, 19, 34, 20, 33] focus on dual SDPs (D). We first
translate the primal SDP (P) into an eigenvalue optimization problem using the
exact penalty method (Proposition 3.1). Then, each iteration of
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP solves a small subproblem formulated
from $r_{\mathrm{p}}$ past eigenvectors and $r_{\mathrm{c}}$ current
eigenvectors of the primal variable $X$ evaluated at the past and current
iterates respectively (Proposition 4.1).

  * •

We show that any configuration of $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP
admits $\mathcal{O}(1/\epsilon^{3})$ convergence rate in primal feasibility,
dual feasibility, and duality gap. Similar to [29],
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP has a faster convergence rate,
$\mathcal{O}(1/\epsilon)$, when the SDPs P and D satisfy strict
complementarity (Theorem 4.1). In Theorem 4.2, we further show linear
convergence of $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP if $(i)$ strict
complementarity holds, $(ii)$ the number of eigenvectors $r_{\mathrm{c}}$ is
larger than the rank of dual optimal solutions. Our proofs largely follow the
strategies in [29, Section 3], [37, Section 4] and [38, Section 7], and we
complete some detailed calculations and handle the constrained case for the
primal SDPs. As a byproduct, we revisit the results for generic bundle methods
in [39, Theorem 2.1 and 2.3] for constrained convex optimization (see Lemma
2.2).

  * •

We present a detailed comparison between the primal and dual formulations of
the spectral bundle methods by showing the symmetry of the parameters and
convergence behaviors from both sides. It becomes clear that the existing dual
formulation in [29, 19, 34, 20, 33] is advantageous when the primal SDP P
admits low-rank solutions. On the other hand, our primal formulation
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP is more suitable when the dual SDP D
admits low-rank solutions. These theoretical findings are supported by a range
of large-scale numerical experiments.

  * •

Finally, we present an open-source implementation of the spectral bundle
algorithms for both P and D, while the existing implementations of spectral
bundle algorithms for D are not open-source or not easily accessible. We
demonstrate that our new spectral bundle method
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP achieves state-of-the-art efficiency
and scalability for solving polynomial optimization compared to a set of
baseline solvers SDPT3 [40], MOSEK [41], CDCS [42], and SDPNAL+ [16].

The rest of the paper is structured as follows. Section 2 covers some
preliminaries on SDPs and nonsmooth optimization. Section 3 presents the exact
penalty formulations for primal and dual SDPs. This is followed by our new
family of spectral bundle methods $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP and
the convergence results in Section 4. Section 5 reviews the classical spectral
bundle methods, called $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD, and the
connections and differences between $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP and
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD are clarified. Our open-source
implementation and numerical experiments are presented in Section 6. Section 7
concludes the paper. Some detailed calculations and technical proofs are
postponed in the appendix.

Notation. We use $\langle\cdot,\cdot\rangle$ to denote the dot product and the
trace inner product on the space of $\mathbb{R}^{n}$ and $\mathbb{S}^{n}$
respectively. For a symmetric matrix $A\in\mathbb{S}^{n}$, we denote its
eigenvalues in the decreasing order as
$\lambda_{\max}(A)=\lambda_{1}(A)\geq\cdots\geq\lambda_{n}(A)$. Given a vector
on $\mathbb{R}^{n}$, we use $\|\cdot\|$ to denote its two norm. For a matrix
$M\in\mathbb{R}^{m\times n}$, its Frobenius norm, operator two norm, and
nuclear norm are denoted by $\|\cdot\|,\|\cdot\|_{\mathrm{op}},$ and
$\|\cdot\|_{*}$ respectively. In the primal and dual SDPs P and D, for
notational simplicity, we will also denote a linear map
$\mathcal{A}:\mathbb{S}^{n}\to\mathbb{R}^{m}$ as
$\mathcal{A}(X):=\begin{bmatrix}\langle A_{1},X\rangle,\ldots,\langle
A_{m},X\rangle\end{bmatrix}^{\mathsf{T}},$ and its adjoint map
$\mathcal{A}^{*}$ that is a linear mapping from $\mathbb{R}^{m}$ to
$\mathbb{S}^{n}$ as $\mathcal{A}^{*}(y):=\sum_{i=1}^{m}A_{i}y_{i}$. The
optimal cost value of SDPs P and D is denoted as $p^{\star}$ and $d^{\star}$
respectively. Finally, given a closed set $\mathcal{C}\subset\mathbb{R}^{n}$
and a point $Y\in\mathbb{R}^{n}$, the distance of $Y$ to $\mathcal{C}$ is
defined as $\mathrm{dist}(Y,\mathcal{C})=\inf_{X\in\mathcal{C}}\|X-Y\|$.

##  2 Preliminaries

In this section, we first introduce standard assumptions and an important
notion of strict complementarity for P and D. We then briefly overview the
exact penalization for constrained nonsmooth convex optimization and the
generic bundle method.

###  2.1 Strict complementarity of SDPs

Throughout this paper, we make the following standard assumptions for well-
behaved SDPs P and D.

######  Assumption 2.1.

The matrices $A_{i},i=1,\ldots,m$ in P and D are linearly independent.

######  Assumption 2.2.

The SDPs P and D satisfy Slater’s constraint qualification, i.e., they are
both strictly feasible.

Assumption 2.1 allows us to uniquely determine $y$ from a given dual feasible
$Z$, i.e., the feasible point $y$ is unique in $Z+\mathcal{A}^{*}(y)=C$ when
giving a feasible $Z$. Under Assumption 2.2, the strong duality holds for P
and D (i.e., $p^{\star}=d^{\star}$), and both P and D are solvable (i.e.,
there exist at least a primal minimizer $X^{\star}$ and a dual maximizer
$(y^{\star},Z^{\star})$ that achieve the optimal cost) [43, Section 5.2.3]. We
denote the set of primal optimal solutions to P as $\mathcal{P}^{\star}$ and
the set of dual optimal solutions to D as $\mathcal{D}^{\star}$, i.e.,

$\displaystyle\mathcal{P}^{\star}$
$\displaystyle=\mathopen{}\mathclose{{}\left\\{X\in\mathbb{S}^{n}\mid
p^{\star}=\langle
C,X\rangle,\mathcal{A}(X)=b,X\in\mathbb{S}^{n}_{+}}\right\\},$ (1a)
$\displaystyle\mathcal{D}^{\star}$
$\displaystyle=\mathopen{}\mathclose{{}\left\\{(y,Z)\in\mathbb{R}^{m}\times\mathbb{S}^{n}\mid
d^{\star}=b^{\mathsf{T}}y,Z+\mathcal{A}^{*}(y)=C,Z\in\mathbb{S}^{n}_{+}}\right\\}.$
(1b)

Assumption 2.2 ensures that $\mathcal{P}^{\star}\neq\emptyset$ and
$\mathcal{D}^{\star}\neq\emptyset$. In addition, if Assumption 2.1 holds, then
the solution sets $\mathcal{P}^{\star}$ and $\mathcal{D}^{\star}$ are nonempty
and compact [29, Section 2].

######  Proposition 2.1 ([29, Section 2]).

Under Assumptions 2.1 and 2.2, the mapping $\mathcal{A}$ is surjective and the
optimal solution sets $\mathcal{P}^{\star}$ and $\mathcal{D}^{\star}$ are
nonempty and compact.

It is clear that $\mathcal{P}^{\star}$ and $\mathcal{D}^{\star}$ are closed.
Given any $(y^{\star},Z^{\star})\in\mathcal{D}^{\star}$ and a strict primal
feasible point $\hat{X}$, we have a finite duality gap $\langle
C,\hat{X}\rangle-b^{\mathsf{T}}y^{\star}=\langle\hat{X},Z^{\star}\rangle\geq
0$. If $\|Z^{\star}\|\to\infty$, then
$\langle\hat{X},Z^{\star}\rangle\to\infty$ since
$Z^{\star}\in\mathbb{S}^{n}_{+}$ and $\hat{X}$ is positive definite. This is
impossible due to the finite duality gap, and thus, any optimal $Z^{\star}$ is
bounded. Assumption 2.1 ensures that $\mathcal{A}$ is a surjective mapping,
which means $\mathcal{A}^{*}$ is injective. Thus, any optimal $y^{\star}$ is
bounded, and $\mathcal{D}^{\star}$ is bounded. Similarly, the existence of a
strictly feasible point $(\hat{y},\hat{Z})$ ensures the compactness of
$\mathcal{P}^{\star}$.

The following result is a version of the KKT optimality condition for the SDPs
P and D.

######  Lemma 2.1 ([44, Lemma 3]).

Given a pair of primal and dual feasible solutions $X^{\star}$ and
$(y^{\star},Z^{\star})$, they are optimal if and only if there exists an
orthonormal matrix $Q\in\mathbb{R}^{n\times n}$ with $Q^{\mathsf{T}}Q=I$, such
that

$\displaystyle
X^{\star}=Q\cdot\mathrm{diag}(\lambda_{1},\ldots,\lambda_{n})\cdot
Q^{\mathsf{T}},\quad Z^{\star}=Q\cdot\mathrm{diag}(w_{1},\ldots,w_{n})\cdot
Q^{\mathsf{T}}$ (2)

and $\lambda_{i}w_{i}=0,i=1,\ldots,n$.

Given a pair of optimal solutions $X^{\star}$ and $(y^{\star},Z^{\star})$, the
complementary slackness condition (2) is equivalent to $Z^{\star}X^{\star}=0$
($X^{\star}$ and $Z^{\star}$ commutes so they share a common set of
eigenvectors as the columns of $Q$). This implies that

$\displaystyle\text{rank}(X^{\star})+\text{rank}(Z^{\star})\leq
n,\quad\text{range}(Z^{\star})\subset\text{null}(X^{\star}),\quad\text{and}\quad\text{range}(X^{\star})\subset\text{null}(Z^{\star}).$

We now introduce the notion of strict complementarity for a pair of optimal
solutions.

######  Definition 2.1 ([44, Definition 4]).

A pair of primal and dual optimal solutions $X^{\star}\in\mathcal{P}^{\star}$
and $(y^{\star},Z^{\star})\in\mathcal{D}^{\star}$ satisfies strict
complementarity if $\mathrm{rank}(X^{\star})+\mathrm{rank}(Z^{\star})=n$
holds, i.e., exactly one of the two conditions $\lambda_{i}=0$ and $w_{i}=0$
is true in (2).

If such a pair $X^{\star}$ and $(y^{\star},Z^{\star})$ exists, we also say
that the SDPs P and D satisfy strict complementarity. We note that strict
complementarity is not restrictive. It is a generic property of SDPs [44,
Theorem 15], and many structured SDPs from practical applications also satisfy
strict complementarity [45].

###  2.2 Exact penalization for constrained convex optimization

Consider a constrained convex optimization problem of the form

$\displaystyle\min$ $\displaystyle\quad f(x)$ (3) subject to
$\displaystyle\quad g_{i}(x)\leq 0,\quad i=1,\ldots,m,$ $\displaystyle\quad
x\in\mathcal{X}_{0},$

where $f:\mathbb{R}^{n}\to\mathbb{R}$ and
$g_{i}:\mathbb{R}^{n}\to\mathbb{R},i=1,\ldots,m$ are (possibly
nondifferentiable) convex functions, and
$\mathcal{X}_{0}\subseteq\mathbb{R}^{n}$ is a closed convex set (which are
defined by some simple constraints). The idea of exact penalty methods is to
reformulate the constrained optimization problem 3 by a problem with simple
constraints. In particular, upon defining an exact penalty function

$P(x)=\sum_{i=1}^{m}\;\max\\{0,g_{i}(x)\\},$

we consider a penalized problem

$\displaystyle\min$ $\displaystyle\quad\Phi_{\rho}(x):=f(x)+\rho P(x)$ (4)
subject to $\displaystyle\quad x\in\mathcal{X}_{0},$

where $\rho>0$ is a penalty parameter. When choosing $\rho$ large enough,
problems (3) and (4) are equivalent to each other in the sense that they have
the same optimal value and solution set.

######  Theorem 2.1 ([38, Theorem 7.21]).

Suppose that problem (3) satisfies Slater’s constraint qualification. There
exists a constant $\rho_{0}\geq 0$ such that for each $\rho>\rho_{0}$, a point
$\hat{x}$ is an optimal solution of (3) if and only if it is an optimal
solution of (4). In particular, we can choose
$\rho_{0}=\sup_{\lambda\in\Lambda}\|\lambda\|_{\infty}$, where
$\Lambda\subset\mathbb{R}^{m}$ is the set of optimal Lagrange multipliers
associated with $g_{i}(x)\leq 0,i=1,\ldots,m$.

Therefore, we can transform some nonsmooth constraints that are hard to handle
in (3) into the nonsmooth cost function of (4). Then, we can apply cutting
plane or bundle methods to solve the nonsmooth optimization (4). However, it
should be noted that the resulting problem (4) may be difficult to solve if
the penalty parameter $\rho$ is too large. As we will discuss in Section 3, in
some SDPs that arise from practical applications, the penalty parameter $\rho$
is known a priori [29, 19, 46].

###  2.3 The cutting plane and bundle methods

In this subsection, we briefly overview the generic bundle method; see [38,
Chapter 7] for details. Consider a generic nonsmooth constrained convex
optimization

$f^{\star}=\min_{x\in\mathcal{X}_{0}}\;\;f(x),$ (5)

where $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is convex but not necessarily
differentiable and $\mathcal{X}_{0}\subseteq\mathbb{R}^{n}$ is a closed convex
set. It is clear that problem 5 includes 4 as a special case.

The simplest method for solving 5 is arguably the subgradient method, which
constructs a sequence of points $x_{t}$ iteratively by updating

$\displaystyle x_{t+1}=\Pi_{\mathcal{X}_{0}}(x_{t}-\tau_{t}g_{t}),\quad
t=1,2,\ldots$

where $g_{t}\in\partial f(x_{t})$ is a subgradient of $f(\cdot)$ at the
current point $x_{t}$, and $\tau_{t}\in\mathbb{R}$ is a step size, and
$\Pi_{\mathcal{X}_{0}}(x)$ denotes the orthogonal projection of the point
$x\in\mathbb{R}^{n}$ onto $\mathcal{X}_{0}$. Recall that for a convex function
$f:\mathbb{R}^{n}\rightarrow\mathbb{R}$, a vector $g\in\mathbb{R}^{n}$ is
called a subgradient of $f$ at $x$ if

$f(y)\geq f(x)+\langle g,y-x\rangle,\;\;\forall y\in\mathbb{R}^{n}.$

The set of all subgradients of $f$ at $x$ is called the subdifferential,
denoted by $\partial f(x)$.

With mild assumptions and an appropriate choice of decreasing step sizes, the
subgradient method is guaranteed to generate a converging sequence
$\\{x_{t}\\}$ to an optimal solution of problem (5) (see [38, Theorem 7.4]).
Despite the simplicity of subgradient methods, it is generally challenging to
develop reliable and efficient step size rules for practical optimization
instances.

#### 2.3.1. The cutting plane method.

Another useful way to utilize subgradients is the idea of the cutting plane
method, which solves a lower approximation of the function $f(x)$ at every
iteration. Here, we assume that $\mathcal{X}_{0}$ is compact in 5; otherwise,
we consider minimizing $f(x)$ over $\mathcal{X}_{0}\cap\mathcal{C}$ where
$\mathcal{C}$ is a compact set containing an optimal solution.

The basic idea of the cutting plane method is to use the subgradient
inequality to construct lower approximations of $f(\cdot)$. In particular, at
iteration $t$, having points $x_{1},x_{2},\ldots,x_{t}$, function values
$f(x_{1}),f(x_{2}),\ldots,$ $f(x_{t})$, and the corresponding subgradients
$g_{1},g_{2},\ldots,g_{t}$, we construct a lower approximation using a piece-
wise affine function

$\displaystyle\hat{f}_{t}(x)=\max_{i=1,\ldots,t}\;\;f(x_{i})+\langle
g_{i},x-x_{i}\rangle.$ (6)

By definition of subgradients, it is clear that
$f(x)\geq\hat{f}_{t}(x),\forall x\in\mathbb{R}^{n}$. Starting from a triple
$\\{x_{1}\in\mathcal{X}_{0},f(x_{1}),g_{1}\in\partial f(x_{1})\\}$, the
cutting plane method solves the following master problem to generate the next
point,

$\displaystyle x_{t+1}$
$\displaystyle\in\operatorname*{\arg\\!\min}_{x\in\mathcal{X}_{0}}\;\;\hat{f}_{t}(x),\quad
t=1,2,\ldots.$

When $\mathcal{X}_{0}$ is a convex set defined by simple constraints (e.g., a
polyhedron), the above problem becomes a linear program (LP) for which very
efficient algorithms exist. The sequence generated by the cutting plane method
is guaranteed to satisfy that $\lim_{t\to\infty}f(x_{t})=f^{\star}$ [38,
Theorem 7.7]. However, the theoretical convergence rate is rather slow; it
generally takes $\mathcal{O}(1/\epsilon^{n})$ iterations to reach
$f(x_{t})-f^{\star}\leq\epsilon$ [47, Section 5.2]. The practical convergence
of the cutting plane method may be faster.

#### 2.3.2. The bundle method.

The bundle method improves the convergence rate and numerical behavior of the
cutting plane method by incorporating a regularization strategy. Unlike the
cutting plane method that only considers the lower approximation function, the
bundle method updates its iterates by solving a regularized master problem
(i.e., a proximal step to the lower approximation model $\hat{f}_{t}(x)$):

$\displaystyle x_{t+1}^{\star}$
$\displaystyle\in\operatorname*{\arg\\!\min}_{x\in\mathcal{X}_{0}}\;\;\hat{f}_{t}(x)+\frac{\alpha}{2}\|x-\omega_{t}\|_{2}^{2},$
(7)

where $\omega_{t}\in\mathbb{R}^{n}$ is the current reference point and
$\alpha>0$ penalizes the deviation from $\omega_{t}$. The bundle method only
updates the iterate $\omega_{t}$ when the decrease of the objective value
$f(\cdot)$ is at least a fraction of the decrease that the approximated model
$\hat{f}_{t}(\cdot)$ predicts. In particular, letting $0<\beta<1$, if

$\displaystyle\beta\mathopen{}\mathclose{{}\left(f(\omega_{t})-\hat{f}_{t}(x_{t+1}^{\star})}\right)\leq
f(\omega_{t})-f(x_{t+1}^{\star})$ (8)

then we set $\omega_{t+1}=x_{t+1}^{\star}$ (descent step); otherwise, we set
$\omega_{t+1}=\omega_{t}$ (null step). In any case, the subgradient
$g_{t+1}\in\partial f(x_{t+1}^{\star})$ at the new point $x_{t+1}^{\star}$ is
used to update the lower approximation $\hat{f}_{t+1}(x)$, e.g., using 6.

The seemingly subtle modifications above have a rather surprising consequence:
the cost value $f(\omega_{t})$ generated by the bundle method converges to
$f^{\star}$ for any constant $\alpha>0$ with a rate of
$\mathcal{O}(1/\epsilon^{3})$ when the objective function is Lipschitz
continuous [48]. Faster convergence rates appear under different assumptions
of $f(x)$; see [39, Table 1] for a detailed comparison. Note that subgradient
methods rely on very carefully controlled decreasing stepsizes which might be
inefficient and unreliable in practice, and the cutting plane method has a
slow convergence rate theoretically. On the contrary, the bundle method
appears more suitable to solve the nonsmooth problem 5.

#### 2.3.3. The bundle method with cut-aggregation.

The lower approximation model $\hat{f}_{t}(x)$ can be constructed using all
past subgradients as in 6, but this leads to a growing number of cuts or
constraints when solving the regularized master problem 7. Another useful cut-
aggregation idea [38, Chapter 7.4.4] allows us to simplify the collection of
$t$ lower bounds used by 6 into just two linear lower bounds. In particular,
the convergence of the bundle method is guaranteed as long as the lower
approximation model $\hat{f}_{t+1}$ satisfies the following three properties
[39]111The analysis in [39] focuses on unconstrained optimization where
$\mathcal{X}_{0}=\mathbb{R}^{n}$ in (5). We have extended the analysis [39]
for constrained optimization where $\mathcal{X}_{0}$ is a closed convex set;
see Appendix A.:

  *   * •

Minorant: the function $\hat{f}_{t+1}$ is a lower bound on $f$, i.e.

$\hat{f}_{t+1}(x)\leq f(x),\quad\forall x\in{\mathcal{X}_{0}}.$ (9a)

  * •

Subgradient lowerbound: $\hat{f}_{t+1}$ is lower bounded by the linearization
given by some subgradient $g_{t+1}\in\partial f(x_{t+1}^{\star})$ computed
after 7, i.e.

$\hat{f}_{t+1}(x)\geq f(x_{t+1}^{\star})+\langle
g_{t+1},x-x_{t+1}^{\star}\rangle,\quad\forall x\in{\mathcal{X}_{0}}.$ (9b)

  * •

Model subgradient lowerbound: $\hat{f}_{t+1}$ is lower bounded by the
linearization of the model $\hat{f}_{t}$ given by the subgradient
$s_{t+1}:=\alpha(\omega_{t}-x_{t+1}^{\star})\in\partial\hat{f}_{t}(x_{t+1}^{\star})+{\mathcal{N}_{\mathcal{X}_{0}}(x_{t+1}^{\star})}$,
i.e.

$\hat{f}_{t+1}(x)\geq\hat{f}_{t}(x_{t+1}^{\star})+\langle
s_{t+1},x-x_{t+1}^{\star}\rangle,\quad\forall x\in{\mathcal{X}_{0}},$ (9c)

where $\mathcal{N}_{\mathcal{X}_{0}}(x_{t+1}^{\star})$ denotes the normal cone
of $\mathcal{X}_{0}$ at point $x_{t+1}^{\star}$, i.e.
$\mathcal{N}_{\mathcal{X}_{0}}(x_{t+1}^{\star})=\\{v\in\mathbb{R}^{n}\mid\langle
v,x-x_{t+1}^{\star}\rangle\leq 0,\forall x\in\mathcal{X}_{0}\\}$. Note that
$s_{t+1}$ certifies the optimality of $x_{t+1}^{\star}$ for problem 7.

The lower bound 9c serves as an aggregation of all previous subgradient lower
bounds. Instead of (6), we can construct the lower approximation
$\hat{f}_{t+1}(x)$ as the maximum of two lower bounds as

$\displaystyle\hat{f}_{t+1}(x)=\max\mathopen{}\mathclose{{}\left\\{f(x_{t+1}^{\star})+\langle
g_{t+1},x-x_{t+1}^{\star}\rangle,\hat{f}_{t}(x_{t+1}^{\star})+\langle
s_{t+1},x-x_{t+1}^{\star}\rangle}\right\\}.$ (10)

Algorithm 1 General Bundle Method

1:Problem $f(x)$, $\mathcal{X}_{0}$; Parameters $\alpha>0$, $\beta\in(0,1)$,
$t_{\max}\geq 1$.

2:for $t=1,2,\ldots,t_{\max}$ do

3: Solve 7 to obtain a candidate iterate $x_{t+1}^{\star}$. $\triangleright$
master problem

4: if 8 is satisfied then

5: Set $\omega_{t+1}=x_{t+1}^{\star}$. $\triangleright$ descent step

6: else

7: Set $\omega_{t+1}=\omega_{t}$. $\triangleright$ null step

8: end if

9: Update approximation model $\hat{f}_{t}$ without violating 9a, 9b and 9c.

10: if stopping criterion then

11: Quit. $\triangleright$ termination

12: end if

13:end for

The overall process of the general bundle method is listed in Algorithm 1. The
convergence rates for Algorithm 1 have been recently revisited in [39, Theorem
2.1 and 2.3]. The big-O notation below suppresses some universal constants.

######  Lemma 2.2 ([39, Theorem 2.1 and 2.3]222The results in [39, Theorem 2.1
and 2.3] focus on unconstrained problems 5 with
$\mathcal{X}_{0}=\mathbb{R}^{n}$. Upon some adaptions, we can extend their
results for constrained convex optimization with a closed convex set
$\mathcal{X}_{0}\subseteq\mathbb{R}^{n}$. Details are presented in Appendix
A.).

Consider a convex and M-Lipschitz function $f(x)$ in 5. Let
$f^{\star}=\inf_{{x}\in\mathcal{X}_{0}}f(x)$ and
$\mathcal{C}=\\{x\in\mathcal{X}_{0}\mid f(x)=f^{\star}\\}$. If $\mathcal{C}$
is nonempty, the number of steps for Algorithm 1 before reaching an
$\epsilon>0$ optimality, i.e. $f(x)-f^{\star}\leq\epsilon$, is bounded by

$\displaystyle t\leq\mathcal{O}\biggl{(}\frac{12\alpha
M^{2}D^{4}}{\beta(1-\beta)^{2}\epsilon^{3}}\biggr{)},$

where $D=\sup_{k}\operatorname{dist}(x_{k},\mathcal{C})<\infty$. If $f(x)$
further satisfies the quadratic growth condition

$\displaystyle
f(x)-f^{\star}\geq\mu\cdot\operatorname{dist}^{2}(x,\mathcal{C}),\quad\forall
x\in\mathcal{X}_{0},$ (11)

where $\mu>0$ is a positive constant, then the number of steps for Algorithm 1
before reaching an $\epsilon$ optimality is bounded by

$\displaystyle
t\leq\mathcal{O}\biggl{(}\frac{16M^{2}}{\beta(1-\beta)^{2}\min\\{\alpha,\mu\\}\epsilon}\biggr{)}.$

When solving SDPs P-D, a specialized version called spectral bundle method
constructs a special lower approximation model satisfying 9a-9c. This idea was
first proposed in [19] for the dual SDP (D) with a constant trace property. In
this paper, we will show spectral bundle methods can be developed for both
general primal and dual SDPs (P) and (D). We will present the details in
Section 4 and Section 5.

##  3 Penalized nonsmooth formulations for SDPs

We here present an exact nonsmooth penalization of primal and dual SDPs P-D in
the form of 5, which allows us to apply the bundle method in Section 4 and
Section 5.

###  3.1 Exact penalization of primal and dual SDPs

The semidefinite constraints in P and D are nonsmooth and typically non-
trivial to deal with for numerical algorithms. A useful method proposed in
[19] is to move nonsmooth semidefinite constraints into the cost function. In
particular, for the primal SDP P, we consider a penalized nonsmooth
formulation

$\displaystyle\min_{X}$ $\displaystyle\langle
C,X\rangle+\rho\max\\{\lambda_{\max}(-X),0\\}$ (12)
$\displaystyle\mathrm{subject\leavevmode\nobreak\ to}$ $\displaystyle\langle
A_{i},X\rangle=b_{i},\quad i=1,\ldots,m,$

and for the dual SDP D, we consider the following penalized nonsmooth
formulation

$\displaystyle\min_{y}\;\;-b^{\mathsf{T}}y+\rho\max\mathopen{}\mathclose{{}\left\\{\lambda_{\max}\mathopen{}\mathclose{{}\left(\sum_{i=1}^{m}A_{i}y_{i}-C}\right),0}\right\\}.$
(13)

From Theorem 2.1, we expect that if the penalty parameter $\rho$ is large
enough, 12 and 13 are equivalent to the primal and dual SDPs P and D,
respectively. We have the following results (recall that $\mathcal{P}^{\star}$
and $\mathcal{D}^{\star}$ are the sets of primal and dual optimal solutions,
respectively; see 1).

######  Proposition 3.1.

Let
$\rho>\mathcal{D}_{Z_{\star}}:=\max_{(y^{\star},Z^{\star})\in\mathcal{D}^{\star}}\;{\mathop{\bf
tr}(Z^{\star})}.$ A point $\hat{X}$ is an optimal solution of the primal SDP P
if and only if it is an optimal solution of (12).

######  Proposition 3.2.

Let
$\rho>\mathcal{D}_{X_{\star}}:=\max_{X^{\star}\in\mathcal{P}^{\star}}\;{\mathop{\bf
tr}(X^{\star})}.$ A point $\hat{y}$ (with
$\hat{Z}=C-\sum_{i=1}^{m}A_{i}\hat{y}_{i}$) is an optimal solution of the dual
SDP D if and only if it is an optimal solution of (13).

Both Propositions 3.1 and 3.2 are direct consequences of Theorem 2.1. In
particular, a proof of Proposition 3.2 appeared in [29]. For completeness, we
provide a proof of Proposition 3.1 in Section B.1. In some applications, we
may have prior information on the bounds of $\mathcal{D}_{Z_{\star}}$ and
$\mathcal{D}_{X_{\star}}$ (for example, one may have explicit trace
constraints; see Section 3.2). In these cases, we can choose the penalty
parameter $\rho$ a priori.

######  Remark 3.1 (Exact penalization for primal SDPs).

The exact penalization for dual SDPs 13 is in the form of unconstrained
eigenvalue minimization; see [49] for excellent discussions on eigenvalue
optimization. To our best knowledge, all the existing results on the
application of bundle methods for solving SDPs focus on the dual formulation
13. One of the early results in [19] assumes a constant trace constraint
$\mathop{\bf tr}(X)=k>0$. This has been generalized to standard SDPs (i.e.
Proposition 3.2) in [29, 18]. However, the exact penalization for primal SDPs
12 has been less studied. We cannot find a formal statement of Proposition 3.1
in the literature. For completeness, we provide a proof of Proposition 3.1 in
Section B.1.

###  3.2 SDPs with trace constraints

Here, we show that the exact penalty formulations can be viewed as standard
SDPs with an explicit trace constraint $\mathop{\bf tr}(X)\leq\rho$ or
$\mathop{\bf tr}(Z)\leq\rho$. In particular, let us consider

$\displaystyle\max_{y,Z}$ $\displaystyle b^{\mathsf{T}}y$ (14)
$\displaystyle\mathrm{subject\leavevmode\nobreak\ to}$ $\displaystyle
Z+\sum_{i=1}^{m}A_{i}y_{i}=C,$ $\displaystyle\mathop{\bf
tr}(Z)\leq\rho,Z\in\mathbb{S}^{n}_{+},$

and

$\displaystyle-\min_{X}$ $\displaystyle\langle C,X\rangle$ (15)
$\displaystyle\mathrm{subject\leavevmode\nobreak\ to}$ $\displaystyle\langle
A_{i},X\rangle=b_{i},\quad i=1,\ldots,m,$ $\displaystyle\mathop{\bf
tr}(X)\leq\rho,X\in\mathbb{S}^{n}_{+}.$

######  Proposition 3.3.

The following statements hold:

  * •

For any $\rho>0$ such that 14 is strictly feasible, the exact penalization 12
and the modified SDP 14 have the same optimal cost value.

  * •

For any $\rho>0$ such that 15 is strictly feasible, the exact penalization 13
and the modified SDP 15 have the same optimal cost value.

###### Proof.

The equivalence comes from the strong duality. We present simple arguments
below. It is straightforward to verify that the Lagrange dual problem for 14
is

$\displaystyle\min_{X,t,Q}$ $\displaystyle\langle C,Q\rangle+\rho t$ (16)
$\displaystyle\mathrm{subject\leavevmode\nobreak\ to}$ $\displaystyle\langle
A_{i},Q\rangle=b_{i},\quad i=1,\ldots,m,$ $\displaystyle Q+tI=X,\;t\geq
0,\;X\in\mathbb{S}^{n}_{+}.$

Eliminating the variable $X$ leads to $Q+tI\in\mathbb{S}^{n}_{+}$, which is
equivalent to $t\geq\lambda_{\max}(-Q)$. Since $t\geq 0$, upon partially
minimizing over $t$, the problem 16 is equivalent to

$\displaystyle\min_{Q}$ $\displaystyle\langle
C,Q\rangle+\rho\max\\{\lambda_{\max}(-Q),0\\}$
$\displaystyle\mathrm{subject\leavevmode\nobreak\ to}$ $\displaystyle\langle
A_{i},Q\rangle=b_{i},\quad i=1,\ldots,m,$

which is clearly equivalent to 12. Since 14 is strictly feasible, strong
duality holds for 14 and 16, which confirms that 12 and 14 have the same
optimal cost value.

Similarly, the Lagrange dual problem for 15 is

$\displaystyle-\max_{y,t}$ $\displaystyle b^{\mathsf{T}}y-t\rho$ (17)
$\displaystyle\mathrm{subject\leavevmode\nobreak\ to}$ $\displaystyle
Z+\sum_{i=1}^{m}A_{i}y_{i}-tI=C,$ $\displaystyle t\geq
0,Z\in\mathbb{S}^{n}_{+}.$

Eliminating the variable $Z$ leads to
$C-\sum_{i=1}^{m}A_{i}y_{i}+tI\in\mathbb{S}^{n}_{+}$, which is equivalent to
$t\geq\lambda_{\max}(\sum_{i=1}^{m}A_{i}y_{i}-C)$. Combining this bound with
the constraint $t\geq 0$ leads to
$t\geq\max\\{\lambda_{\max}(\sum_{i=1}^{m}A_{i}y_{i}-C),0\\}$. Thus, the
problem 17 is equivalent to

$\displaystyle-\max_{y}\;\;$ $\displaystyle
b^{\mathsf{T}}y-\rho\max\Biggl{\\{}{\lambda_{\max}\mathopen{}\mathclose{{}\left(\sum_{i=1}^{m}A_{i}y_{i}-C}\right),0}\Biggr{\\}},$

which is equivalent to 13. Since the strongly duality holds for 15 and 17, we
know that 13 and 15 have the same optimal cost value. ∎

We note that Proposition 3.3 implies that when $\rho$ is large enough, i.e.,
satisfying the bounds in Propositions 3.1 and 3.2, 15 is equivalent to the
primal SDP P, and 14 is equivalent to the dual SDP D. This result becomes
obvious since the extra constraint trace constraint $\mathop{\bf
tr}(X)\leq\rho$ or $\mathop{\bf tr}(Z)\leq\rho$ does not affect the optimal
solutions.

If the constraints $\langle A_{i},X\rangle=b_{i},i=1,\ldots,m$ imply that
$\mathop{\bf tr}(X)=k$ for some $k>0$, i.e., the primal SDP P has an implicit
constant trace constraint, the exact dual penalization 13 can be simplified as

$\displaystyle\min_{y}\;\;-b^{\mathsf{T}}y+k\lambda_{\max}\mathopen{}\mathclose{{}\left(\sum_{i=1}^{m}A_{i}y_{i}-C}\right).$
(18)

This was first used to derive the original spectral bundle method in [19].
Similarly, if $Z+\sum_{i=1}^{m}A_{i}y_{i}=C$ imply that $\mathop{\bf
tr}(Z)=k$, i.e., $\mathop{\bf tr}(A_{i})=0$, the exact primal penalization 12
can be simplified as

$\displaystyle\min_{X}$ $\displaystyle\langle C,X\rangle+k\lambda_{\max}(-X)$
(19) $\displaystyle\mathrm{subject\leavevmode\nobreak\ to}$
$\displaystyle\langle A_{i},X\rangle=b_{i},\quad i=1,\ldots,m.$

For self-completeness, we provide short derivations for 18 and 19 in Section
B.2.

######  Remark 3.2 (Applications with constant trace).

We note that primal SDPs with a constant trace constraint are very common for
semidefinite relaxations of binary combinatorial optimization problems, such
as MaxCut [50] and Lov$\acute{a}$sz theta number [51]. Also, dual SDPs with a
constant trace constraint appears in certain matrix completion problem [52]
(see Example B.1) and the moment/sum-of-squares relaxation of polynomial
optimizations [46] (see Sections 6.4 and E.2). For these problems, the penalty
parameter $\rho$ is thus known a priori.

##  4 Spectral bundle methods for primal SDPs

We can apply the standard bundle method in Section 2.3 to solve the penalized
primal formulation 12 or dual formulation 13. This idea was first proposed in
[19], and further revised and developed in [20, 33, 29, 34]. To our best
knowledge, however, all previous studies [29, 19, 34, 20, 33] only consider
the penalized dual formulation 13. The dual formulation is in the form of
unconstrained eigenvalue optimization [49], for which it seems more convenient
to apply the bundle method.

In this section, we apply the bundle method to solve the penalized primal
formulation 12, which leads to a new family of spectral bundle algorithms.
Differences and connections between our new algorithms and the existing
spectral bundle algorithms will be clarified in Section 5.

###  4.1 A new family of spectral bundle algorithms for primal SDPs

For notational convenience, we denote the cost function in 12 as

$\displaystyle F(X):=\langle C,X\rangle+\rho\max\\{\lambda_{\max}(-X),0\\}.$
(20)

Directly applying the bundle method in Section 2.3 to solve 12 requires
computing a subgradient of $F(X)$ at every iteration $t$. It is known that for
every $X_{t}\in\mathbb{S}^{n}$, a subgradient $g_{t}\in\partial F(X_{t})$ is
given by [38, Example 2.89]

$\displaystyle g_{t}=\begin{cases}C-\rho v_{t}v_{t}^{\mathsf{T}},&\text{ if
}\lambda_{\max}(-X_{t})>0,\\\ C,&\text{ otherwise},\end{cases}$

where $v_{t}\in\mathbb{R}^{n}$ is a normalized eigenvector associated with
$\lambda_{\max}(-X_{t})$.

#### 4.1.1. Lower approximation models.

As discussed in Section 2.3.3, one key step in the bundle method is to
construct a valid lower approximation model of $F(X)$ at each iteration $t$.
Similar to (6), one natural choice is to use a piece-wise affine function,

$\displaystyle\hat{F}_{t}(X)$ $\displaystyle=\max_{i=1,\ldots,t}\;\;\langle
C,X_{i}\rangle+\rho\max\\{\lambda_{\max}(-X_{i}),0\\}+\langle
g_{i},X-X_{i}\rangle,$

where $X_{i},i=1,\ldots,t$ are the past iterates, and $g_{i}\in\partial
F(X_{i}),i=1,\ldots,t$ are subgradients. Via simple derivations, each affine
function corresponding to $X_{i}$ becomes

$\langle C,X_{i}\rangle+\rho\max\\{\lambda_{\max}(-X_{i}),0\\}+\langle
g_{i},X-X_{i}\rangle\\\ \\!=\\!\begin{cases}\langle
C,X_{i}\rangle+\rho\mathopen{}\mathclose{{}\left\langle
v_{i}v_{i}^{\mathsf{T}},-X}\right\rangle,&\text{if
}\lambda_{\max}(-X_{i})>0,\\\ \langle
C,X_{i}\rangle,&\text{otherwise},\end{cases}\\\ $ (21)

where $v_{i}\in\mathbb{R}^{n}$ is a normalized eigenvector corresponding to
$\lambda_{\max}(-X_{i})$.

For this special cost function, one key idea of the original spectral bundle
method in [19] is to improve the lower bound 21 using infinitely many affine
minorants. In particular, at iteration $t$, we compute a matrix
$P_{t}\in\mathbb{R}^{n\times r}$ with some small value $1\leq r<n$ and
orthonormal columns (i.e., $P_{t}^{\mathsf{T}}P_{t}=I\in\mathbb{S}^{r}$), and
define a lower approximation function

$\displaystyle\hat{F}_{P_{t}}(X)$ $\displaystyle=\langle
C,X\rangle+\rho\max_{S\in\mathbb{S}^{r}_{+},\mathop{\bf tr}(S)\leq
1}\mathopen{}\mathclose{{}\left\langle
P_{t}SP_{t}^{\mathsf{T}},-X}\right\rangle.$ (22)

It is clear that $F(X)\geq\hat{F}_{P_{t}}(X),\forall X\in\mathbb{S}^{n}$
thanks to the fact that

$\displaystyle\max\\{\lambda_{\max}(-X),0\\}=\max_{S\in\mathbb{S}^{n}_{+},\mathop{\bf
tr}(S)\leq 1}\langle S,-X\rangle,\qquad\forall X\in\mathbb{S}^{n},$

and

$\displaystyle\\{P_{t}SP_{t}^{\mathsf{T}}\in\mathbb{S}^{n}_{+}\mid
S\in\mathbb{S}^{r}_{+},\mathop{\bf tr}(S)\leq
1\\}\subset\\{S\in\mathbb{S}^{n}_{+}\mid\mathop{\bf tr}(S)\leq 1\\}.$

Meanwhile, it is not difficult to check that if $r=1$, and $P_{t}=v_{t}$ with
$v_{t}$ being the top eigenvector of $-X_{t}$, then $\hat{F}_{P_{t}}(X)$
defined in 22 is reduced to the approximation function in 21 with $i=t$. Thus,
when choosing $r>1$ and selecting $P_{t}$ spanning $v_{t}$, we have a strictly
better lower approximation using 22 than the simple linear function 21 based
on one subgradient. In principle, the columns of $P_{t}\in\mathbb{R}^{n\times
r}$ should consist of both top eigenvectors of the current iterate $X_{t}$ and
the accumulation of spectral information from past iterates
$X_{1},\ldots,X_{t-1}$ [19, 29].

Therefore, by construction, $\hat{F}_{P_{t}}(X)$ in 22 is naturally a minorant
satisfying 9a, and it also satisfies subgradient lowerbound 9b. However, a
further refinement is needed for 22 to fulfill the model subgradient lower
bound 9c. The spectral bundle method [19, 29] maintains a carefully selected
weight to capture past information. In particular, we introduce a matrix
$\bar{W}_{t}\in\mathbb{S}^{n}_{+}$ and $\mathop{\bf tr}(\bar{W}_{t})=1$, and
then build the lower approximated model using $\bar{W}_{t}$ and
$P_{t}\in\mathbb{R}^{n\times r}$,

$\displaystyle\hat{F}_{(\bar{W}_{t},P_{t})}(X)=\langle
C,X\rangle+\rho\max_{S\in\mathbb{S}^{r}_{+},\gamma\geq 0,\gamma+\mathop{\bf
tr}(S)\leq
1}\mathopen{}\mathclose{{}\left\langle\gamma\bar{W}_{t}+P_{t}SP_{t}^{\mathsf{T}},-X}\right\rangle.$
(23)

It is clear that this lower approximation model 23 is an improved
approximation than 22 (e.g., letting $\gamma=0$ in 23 recovers the
approximation model 22). Thus, it satisfies the inequalities 9a and 9b. Upon
careful construction of $\bar{W}_{t}$ at each iteration, we will show that 23
will also satisfy 9c. The construction details will be presented below.

#### 4.1.2. Spectral bundle algorithms.

We are ready to introduce a new family of spectral bundle algorithms for
primal SDPs, which we call $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP333This name
is consistent with [29] which focuses on solving the dual penalization
formulation 13.. As we will detail below, the family of spectral bundle
algorithms considers $r=r_{\mathrm{p}}+r_{\mathrm{c}}$ (where
$r_{\mathrm{p}}\geq 0,r_{\mathrm{c}}\geq 1$) normalized eigenvectors to form
the orthonormal matrix $P_{t}\in\mathbb{R}^{n\times r}$ that is used in 23.
The algorithms have the following steps.

Initialization: $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP starts with an initial
guess $\Omega_{0}\in\mathbb{S}^{n}$, and $P_{0}\in\mathbb{R}^{n\times{r}}$
formed by the $r=r_{\mathrm{p}}+r_{\mathrm{c}}$ top normalized eigenvectors of
$-\Omega_{0}$, any weight matrix $\bar{W}_{0}\in\mathbb{S}^{n}_{+}$ with
$\mathop{\bf tr}(\bar{W}_{0})=1$. Matrices $P_{0}$ and $\bar{W}_{0}$ are used
to construct an initial lower approximation model
$\hat{F}_{(\bar{W}_{0},P_{0})}(X)$ in 23.

Solve the master problem: Similar to 7, at iteration $t\geq 0$,
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP solves the following regularized master
problem

$\displaystyle(X_{t+1}^{\star},S_{t}^{\star},\gamma^{\star}_{t})=\operatorname*{\arg\\!\min}_{X\in\mathcal{X}_{0}}\;\;\hat{F}_{(\bar{W}_{t},P_{t})}(X)+\frac{\alpha}{2}\|X-\Omega_{t}\|^{2},$
(24)

where $\Omega_{t}\in\mathbb{S}^{n}$ is the current reference point (proximal
center), $\alpha>0$ serves as a penalty for the deviation from $\Omega_{t}$,
and $\mathcal{X}_{0}:=\\{X\in\mathbb{S}^{n}\;|\;\mathcal{A}(X)=b\\}.$ Solving
24 is the main computation in each iteration of
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP, and we provide its computational
details in Section 4.2.

Update reference point: Similar to 8, $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP
updates the next reference point $\Omega_{t+1}$ as follows: given
$\beta\in(0,1)$, if

$\displaystyle\beta\mathopen{}\mathclose{{}\left(F(\Omega_{t})-\hat{F}_{(\bar{W}_{t},P_{t})}(X_{t+1}^{\star})}\right)\leq
F(\Omega_{t})-F(X_{t+1}^{\star})$ (25)

holds true (i.e. at the candidate point $X_{t+1}^{\star}$, the decrease of the
objective value $F(\cdot)$ is at least $\beta$ fraction of the decrease in
objective value that the model $\hat{F}_{(\bar{W}_{t},P_{t})}(\cdot)$
predicts), we set $\Omega_{t+1}=X_{t+1}^{\star}$, which is called a descent
step. Otherwise, we let $\Omega_{t+1}=\Omega_{t}$, which is called a null
step.

Update the lower approximation model: $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP
updates the spectral matrices $\bar{W}_{t+1},P_{t+1}$ for the lower
approximation model $\hat{F}_{(\bar{W}_{t+1},P_{t+1})}(\cdot)$ using a
strategy similar to that in [29, 19]. We first compute the eigenvalue
decomposition of the small $r\times r$ matrix $S^{\star}_{t}$ as

$S^{\star}_{t}=\begin{bmatrix}Q_{1}&Q_{2}\end{bmatrix}\begin{bmatrix}\Sigma_{1}&0\\\
0&\Sigma_{2}\end{bmatrix}\begin{bmatrix}Q_{1}^{\mathsf{T}}\\\
Q_{2}^{\mathsf{T}}\end{bmatrix},$

where $Q_{1}\in\mathbb{R}^{r\times r_{\mathrm{p}}}$ consists of the top
$r_{\mathrm{p}}\geq 0$ orthonormal eigenvectors,
$\Sigma_{1}\in\mathbb{S}^{r_{\mathrm{p}}}$ is a diagonal matrix formed by the
top $r_{\mathrm{p}}$ eigenvalues of $S^{\star}_{t}$, and
$Q_{2}\in\mathbb{R}^{r\times(r-r_{\mathrm{p}})}$ and
$\Sigma_{2}\in\mathbb{S}^{(r-r_{\mathrm{p}})}$ captures the remaining
orthonormal eigenvectors and eigenvalues, respectively.

  * •

The orthonormal matrix $P_{t+1}$: we compute $V_{t}\in\mathbb{R}^{n\times
r_{\mathrm{c}}}$ with its columns being the top $r_{\mathrm{c}}\geq 1$
eigenvectors of $-X_{t+1}^{\star}$, which naturally contains the subgradient
information of the true objective function $F(X)$ at $X_{t+1}^{\star}$. Let
the range space of $P_{t+1}$ span $V_{t}$, which guarantees to improve the
lower approximation model. We also let $P_{t+1}$ contain the past important
information $P_{t}Q_{1}$ according to the top $r_{\mathrm{p}}$ eigenvalues of
$S^{\star}_{t}$. Therefore, we update $P_{t+1}$ as

$P_{t+1}=\mathrm{orth}\mathopen{}\mathclose{{}\left(\begin{bmatrix}V_{t}&P_{t}Q_{1}\end{bmatrix}}\right),$
(26)

where $\mathrm{orth}(\cdot)$ denotes an orthonormal process such that
$P_{t+1}^{\mathsf{T}}P_{t+1}=I_{r}$ with $r=r_{\mathrm{p}}+r_{\mathrm{c}}$.

  * •

The weight matrix $\bar{W}_{t+1}$: we keep the rest of the past information by
updating $\bar{W}_{t+1}$ as

$\bar{W}_{t+1}=\frac{1}{\gamma^{\star}_{t}+\mathop{\bf
tr}(\Sigma_{2})}\mathopen{}\mathclose{{}\left(\gamma^{\star}_{t}\bar{W}_{t}+P_{t}Q_{2}\Sigma_{2}Q_{2}^{\mathsf{T}}P_{t}^{\mathsf{T}}}\right).$
(27)

Note that $\bar{W}_{t+1}$ has been normalized such that $\mathop{\bf
tr}(\bar{W}_{t+1})=1$.

If $r_{\mathrm{p}}=0$, the updates in 26 and 27 become

$P_{t+1}=V_{t},\qquad\text{and}\qquad\bar{W}_{t+1}=\frac{1}{\mathop{\bf
tr}(W_{t}^{\star})}W_{t}^{\star},$

where we denote the optimal solution of
$\gamma\bar{W}_{t}+P_{t}SP_{t}^{\mathsf{T}}$ in 24 as

$W_{t}^{\star}=\gamma^{\star}_{t}\bar{W}_{t}+P_{t}S^{\star}_{t}P_{t}^{\mathsf{T}}.$
(28)

Overall, $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP generates a sequence of points
$\\{\Omega_{t},W_{t}^{\star},y_{t}^{\star}\\}$, where $y_{t}^{\star}$ is the
dual variable corresponding to the affine constraint in 24, and a sequence of
monotonically decreasing cost values $\\{F(\Omega_{t})\\}$. The detailed steps
of $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP are listed in Algorithm 2.

Algorithm 2 $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP– Spectral Bundle Method for
Primal SDPs

1:Problem data $A_{1},\ldots,A_{m},C\in\mathbb{S}^{n}$, $b\in\mathbb{R}^{n}$.
Parameters $r_{\mathrm{p}}\geq 0,r_{\mathrm{c}}\geq 1$, $\alpha>0$,
$\beta\in(0,1)$, $\epsilon\geq 0$, $t_{\max}\geq 1$. An initial point
$\Omega_{0}\in\mathbb{S}^{n}$.

2:Initialization: Let ${r}=r_{\mathrm{p}}+r_{\mathrm{c}}$. Initialize
$\bar{W}_{0}\in\mathbb{S}^{n}_{+}$ with $\mathop{\bf tr}(\bar{W}_{0})=1$.
Compute $P_{0}\in\mathbb{R}^{n\times{r}}$ with columns being the top ${r}$
orthonormal eigenvectors of $-\Omega_{0}$.

3:for $t=0,\ldots,t_{\max}$ do

4: Solve 24 to obtain $X_{t+1}^{\star},\gamma^{\star}_{t},S^{\star}_{t}$.
$\triangleright$ master problem

5: Form the iterate $W_{t}^{\star}$ in 28 and dual iterate $y_{t}^{\star}$ in
29.

6: if $t=0$ and $\mathcal{A}(\Omega_{t})\neq b$ then

7: Set primal iterate $\Omega_{t+1}=X_{t+1}^{\star}$. $\triangleright$ affine
constraint

8: else

9: if 25 holds then

10: Set primal iterate $\Omega_{t+1}=X_{t+1}^{\star}$. $\triangleright$
descent step

11: else

12: Set primal iterate $\Omega_{t+1}=\Omega_{t}$. $\triangleright$ null step

13: end if

14: end if

15: Compute $P_{t+1}$ as 26, and $\bar{W}_{t+1}$ as 27. $\triangleright$
update model

16: if stopping criterion then

17: Quit. $\triangleright$ termination

18: end if

19:end for

###  4.2 Computational details

At every iteration $t$, we need to solve the subproblem 24, which is the main
computation in $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP. Therefore, it is
crucial to solve the master problem 24 efficiently. We summarize the
computational details of solving 24 in the following proposition.

######  Proposition 4.1.

The master problem 24 is equivalent to the following problem

$\displaystyle\min_{W\in\hat{\mathcal{W}}_{t},\leavevmode\nobreak\
y\in\mathbb{R}^{m}}\;\;\langle W-C,\Omega_{t}\rangle-\langle
b-\mathcal{A}(\Omega_{t}),y\rangle+\frac{1}{2\alpha}\mathopen{}\mathclose{{}\left\|W-C+\sum_{i=1}^{m}A_{i}y_{i}}\right\|^{2},$
(29)

where the constraint set is defined as

$\hat{\mathcal{W}}_{t}:=\\{\gamma\bar{W}_{t}+P_{t}SP_{t}^{\mathsf{T}}\in\mathbb{S}^{n}\;|\;S\in\mathbb{S}^{r}_{+},\gamma\geq
0,\gamma+\mathop{\bf tr}(S)\leq\rho\\}.$ (30)

The optimal $X$ in 24 is recovered by

$\displaystyle
X_{t+1}^{\star}=\Omega_{t}+\frac{1}{\alpha}\mathopen{}\mathclose{{}\left(W_{t}^{\star}-C+\mathcal{A}^{*}(y_{t}^{\star})}\right),$
(31)

where $(W_{t}^{\star},y_{t}^{\star})$ is a minimizer of 29.

###### Proof.

Our proof relies on the strong duality of convex optimization. Upon applying
the definition of $\hat{F}_{(\bar{W}_{t},P_{t})}(\cdot)$ in (23), it is clear
that 24 becomes

$\displaystyle\min_{X\in\mathcal{X}_{0}}\langle
C,X\rangle+\rho\max_{S\in\mathbb{S}^{r}_{+},\gamma\geq 0,\gamma+\mathop{\bf
tr}(S)\leq
1}\;\;\langle\gamma\bar{W}_{t}+P_{t}SP_{t}^{\mathsf{T}},-X\rangle+\frac{\alpha}{2}\|X-\Omega_{t}\|^{2}$
$\displaystyle=$
$\displaystyle\min_{X\in\mathcal{X}_{0}}\max_{S\in\mathbb{S}^{r}_{+},\gamma\geq
0,\gamma+\mathop{\bf tr}(S)\leq\rho}\;\;\langle
C,X\rangle+\langle\gamma\bar{W}_{t}+P_{t}SP_{t}^{\mathsf{T}},-X\rangle+\frac{\alpha}{2}\|X-\Omega_{t}\|^{2}$
$\displaystyle=$
$\displaystyle\min_{X\in\mathcal{X}_{0}}\max_{W\in\hat{\mathcal{W}}_{t}}\;\;\langle
C-W,X\rangle+\frac{\alpha}{2}\|X-\Omega_{t}\|^{2},$

where the first equality brings the constant $\rho$ into the constraint, and
the second equality applies a change of variables
$W=\gamma\bar{W}_{t}+P_{t}SP_{t}^{\mathsf{T}}$ and uses the set
$\hat{\mathcal{W}}_{t}$ defined as 30. Since $\hat{\mathcal{W}}_{t}$ is
bounded, by strong duality [53, Corollary 37.3.2], we can switch the min-max
order and have the following equivalency

$\displaystyle\min_{X\in\mathcal{X}_{0}}\max_{W\in\hat{\mathcal{W}}_{t}}\;\;\langle
C-W,X\rangle+\frac{\alpha}{2}\|X-\Omega_{t}\|^{2}=\max_{W\in\hat{\mathcal{W}}_{t}}\min_{X\in\mathcal{X}_{0}}\;\;\langle
C-W,X\rangle+\frac{\alpha}{2}\|X-\Omega_{t}\|^{2}.$

Note that the inner minimization is an equality-constrained quadratic program,

$\min_{X\in\mathcal{X}_{0}}\;\;\langle
C-W,X\rangle+\frac{\alpha}{2}\|X-\Omega_{t}\|^{2},$ (32)

which can be simplified by considering its dual formulation. Specifically, we
introduce a dual variable $y\in\mathbb{R}^{m}$ and construct the Lagrangian
for 32 as follows

$\displaystyle L(X,y)$ $\displaystyle=\langle
C-W,X\rangle+\frac{\alpha}{2}\|X-\Omega_{t}\|^{2}+y^{\mathsf{T}}(b-\mathcal{A}(X)),$

which is strongly convex in $X$. The dual function for 32 is given by
$g(y):=\min_{X}L(X,y)$, where the unique minimizer $X$ is

$\displaystyle
X=\Omega_{t}+\frac{1}{\alpha}\mathopen{}\mathclose{{}\left(W-C+\sum_{i=1}^{m}A_{i}y_{i}}\right).$
(33)

Therefore, the dual function becomes

$g(y)=\langle C-W,\Omega_{t}\rangle+\langle
b-\mathcal{A}(\Omega_{t}),y\rangle-\frac{1}{2\alpha}\mathopen{}\mathclose{{}\left\|W-C+\sum_{i=1}^{m}A_{i}y_{i}}\right\|^{2}.$

By strong duality of 32, we have

$\displaystyle\max_{W\in\hat{\mathcal{W}}_{t}}\min_{X\in\mathcal{X}_{0}}\;\;\langle
C-W,X\rangle+\frac{\alpha}{2}\|X-\Omega_{t}\|^{2}$ $\displaystyle=$
$\displaystyle\max_{W\in\hat{\mathcal{W}}_{t}}\max_{y\in\mathbb{R}^{m}}\;\;\langle
C-W,\Omega_{t}\rangle+\langle
b-\mathcal{A}(\Omega_{t}),y\rangle-\frac{1}{2\alpha}\mathopen{}\mathclose{{}\left\|W-C+\sum_{i=1}^{m}A_{i}y_{i}}\right\|^{2}$
$\displaystyle=$
$\displaystyle\max_{W\in\hat{\mathcal{W}}_{t},y\in\mathbb{R}^{m}}\;\;\langle
C-W,\Omega_{t}\rangle+\langle
b-\mathcal{A}(\Omega_{t}),y\rangle-\frac{1}{2\alpha}\mathopen{}\mathclose{{}\left\|W-C+\sum_{i=1}^{m}A_{i}y_{i}}\right\|^{2},$

which is clearly equivalent to 29. Finally, the optimal $X$ in 31 is recovered
in 33 once we obtain the optimal dual variables $y_{t}^{\star}$ and
$W_{t}^{\star}$ from solving 29. This completes the proof. ∎

After the first iteration, if $\mathcal{A}(\Omega_{0})\neq b$, we update the
proximal center $\Omega_{1}=X^{\star}_{1}$. Then, the rest of the iterations
are naturally feasible to the affine constraint, i.e.,
$A(\Omega_{t})=b,\leavevmode\nobreak\ \forall t>0$. Therefore, the master
problem 29 can be further simplified as

$\displaystyle\min_{W\in\hat{\mathcal{W}}_{t},\leavevmode\nobreak\
y\in\mathbb{R}^{m}}\;\;\langle
W-C,\Omega_{t}\rangle+\frac{1}{2\alpha}\mathopen{}\mathclose{{}\left\|W-C+\sum_{i=1}^{m}A_{i}y_{i}}\right\|^{2}.$

######  Remark 4.1 (Complexity of solving 29 in
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP).

The subproblem 29 is a semidefinite program with a convex quadratic cost
function. The dimension of the PSD constraint is
$r=r_{\mathrm{p}}+r_{\mathrm{c}}$, which can be chosen to be very small (i.e.,
$r\ll n$). Thus, 29 can be efficiently solved using either standard conic
solvers (such as SeDuMi [54] and Mosek [41]) or customized interior-point
algorithms [19]. In addition, we note that the dual variable
$y\in\mathbb{R}^{m}$ in 29 admits an analytical closed-loop solution in terms
of $W$, and thus 29 can be further converted into the following form

$\displaystyle\min_{v\in\mathbb{R}^{1+r^{2}}}$ $\displaystyle
v^{\mathsf{T}}Qv+q^{\mathsf{T}}v+c$ (34)
$\displaystyle\mathrm{subject\leavevmode\nobreak\ to}$ $\displaystyle
v=\begin{bmatrix}\gamma&\mathrm{vec}\mathopen{}\mathclose{{}\left(S}\right)^{\mathsf{T}}\end{bmatrix}^{\mathsf{T}},$
$\displaystyle\gamma\geq 0,S\in\mathbb{S}^{r}_{+},\gamma+\mathop{\bf
tr}(S)\leq\rho,$

where $c\in\mathbb{R}$ is a constant, and $Q\in\mathbb{S}^{1+r^{2}}$ and
$q\in\mathbb{R}^{1+r^{2}}$ depend on the problem data and the matrices $P_{t}$
and $\bar{W}_{t}$ at step $t$. We present the details of transforming 29 into
34 in Section C.1. Therefore, at each iteration of
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP, we only need to solve (34), which
admits very efficient solutions since $r$ can be much smaller than the
original dimension $n$. As we shall see in Appendix C, the computational
complexity of solving the regularized sub-problem 24 in spectral bundle
methods for primal SDPs is very similar to that in spectral bundle methods for
dual SDPs. Furthermore, when $r=1$ (i.e.,
$r_{\mathrm{p}}=0,r_{\mathrm{c}}=1$), the problem 34 admits an analytical
closed-form solution, and no other solver is required; see Section C.3.

###  4.3 Convergence results

We here present two convergence guarantees Theorem 4.1 and Theorem 4.2 for
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP. Consistent with Lemma 2.2, when strong
duality holds for P and D, Theorem 4.1 below provides a convergence rate of
$\mathcal{O}(1/\epsilon^{3})$ in terms of cost value gap, approximate primal
feasibility, approximate dual feasibility, and approximate primal-dual
optimality gap. The convergence rate improves to $\mathcal{O}(1/\epsilon)$
under the condition of strict complementarity (see Definition 2.1).

######  Theorem 4.1.

Suppose Assumptions 2.1 and 2.2 are satisfied. Given any $\beta\in(0,1)$,
$r_{\mathrm{c}}\geq 1$, $r_{\mathrm{p}}\geq 0$, $\alpha>0$,
$r=r_{\mathrm{c}}+r_{\mathrm{p}}$, $\rho>2\mathcal{D}_{Z_{\star}}+1$, and
target accuracy $\epsilon>0$, the $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP in
Algorithm 2 produces iterates $\Omega_{t}$, $W_{t}^{\star}$, and
$y_{t}^{\star}$ with $F(\Omega_{t})-F(X^{\star})\leq\epsilon$ and

approximate primal feasibility:
$\displaystyle\quad\mathcal{A}(\Omega_{t})-b=0,\quad\lambda_{\min}(\Omega_{t})\geq-\epsilon,$
(35a) approximate dual feasibility:
$\displaystyle\quad\mathopen{}\mathclose{{}\left\|W_{t}^{\star}+\mathcal{A}^{*}(y_{t}^{\star})-C}\right\|^{2}\leq\epsilon,\quad
W_{t}^{\star}\succeq 0,$ (35b) approximate primal-dual optimality:
$\displaystyle\quad|\langle C,\Omega_{t}\rangle-\langle
b,y_{t}^{\star}\rangle|\leq\sqrt{\epsilon}$ (35c)

by some iteration $t\leq\mathcal{O}(1/\epsilon^{3})$. If additionally, strict
complementarity holds for P and D, then these conditions 35 are reached by
some iteration $t\leq\mathcal{O}(1/\epsilon)$.

In addition to strict complementarity, an improved convergence rate can be
established if the number of current eigenvectors at every iteration satisfies

$r_{\mathrm{c}}\geq\sigma_{\mathrm{p}}:=\max_{X^{\star}\in\mathcal{P}^{\star}}\dim(\text{null}(X^{\star}))$
(36)

with proper choice of $\alpha$ and $\beta$. Under these conditions, Theorem
4.2 ensures that $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP converges linearly
once the iterate $\Omega_{t}$ is close enough to the set of the primal optimal
solutions $\mathcal{P}^{\star}$.

######  Theorem 4.2.

Suppose Assumptions 2.1 and 2.2 are satisfied and strict complementarity holds
for P and D. There exist constants $T_{0}>0$ and $\eta>0$. Under a proper
selection of $\alpha\geq\eta$, $\beta\in(0,\frac{1}{2}]$,
$r_{\mathrm{c}}\geq\sigma_{\mathrm{p}}$, $r_{\mathrm{p}}\geq 0$,
$r=r_{\mathrm{c}}+r_{\mathrm{p}}$, $\rho>2\mathcal{D}_{Z_{\star}}+1$, and
target accuracy $\epsilon>0$, after at most ${T}_{0}$ iterations, the
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP in Algorithm 2 only takes descent steps
and converges linearly to an optimal solution. Consequently,
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP produces iterates $\Omega_{t}$,
$W_{t}^{\star}$, and $y_{t}^{\star}$ satisfying
$F(\Omega_{t})-F(X^{\star})\leq\epsilon$ and 35 by at most
${T}_{0}+\mathcal{O}(\log(1/\epsilon))$ iterations.

The proof sketches are provided in Section 4.4. The constants $T_{0}$ and
$\eta$ only depend on problem data and are independent of the sub-optimality
$\epsilon$, and we provide some discussions in Section 4.4.2 (see Lemma D.5
and Section D.4 in the appendix for further details).

######  Remark 4.2.

The convergence results in Theorems 4.1 and 4.2 can be viewed as the
counterparts of [29, Theorems 3.1 and 3.2] when applying the spectral bundle
method for solving primal SDPs. Note that [29, Theorems 3.1 and 3.2] focus on
solving dual SDPs only (we will review some details in Section 5.1). As
highlighted earlier, all existing studies [19, 20, 33, 29] only consider the
penalized dual formulation 13. Here, we establish a class of spectral bundle
methods, i.e., $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP in Algorithm 2, to
directly solve primal SDPs with similar computational complexity and
convergence behavior. However, we remark that the value of $r_{\mathrm{c}}$ in
(36) can be drastically different from that in [29] for linear convergence. A
detailed comparison will be presented in Section 5.2.

###  4.4 Proof sketches

Here, we provide some proof sketches for Theorems 4.1 and 4.2. The proofs
largely follow the strategies in [29, Section 3], [37, Section 4] and [38,
Section 7]. We complete some detailed calculations, handle the constrained
case for the penalized primal SDPs (12) (note that the penalized dual case is
unconstrained), and fix minor typos in [29, 37, 38]. We have also provided an
extension in Lemma 4.3 compared to [29, Section 3]. We do not claim main
contributions for establishing the proofs, and we provide them for self-
completeness and for the convenience of interested readers.

#### 4.4.1. Proof of Theorem 4.1

One main step in the proof of Theorem 4.1 is to characterize the improvements
in terms of primal feasibility, dual feasibility, and primal-dual optimality
at every descent step, which is summarized in the following lemma. Its proof
is provided in Appendix D.

######  Lemma 4.1.

In $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP, let $\beta\in(0,1)$,
$r_{\mathrm{c}}\geq 1$, $r_{\mathrm{p}}\geq 0$, $\alpha>0$,
$r=r_{\mathrm{c}}+r_{\mathrm{p}}$, $\rho>2\mathcal{D}_{Z_{\star}}+1$. Then, at
every descent step $t>0$, the following results hold.

  1. 1.

The approximate primal feasibility for $\Omega_{t+1}$ satisfies

$\displaystyle\lambda_{\min}(\Omega_{t+1})\geq\frac{-(F(\Omega_{t})-F(X^{\star}))}{\mathcal{D}_{Z_{\star}}+1},\quad
and\quad\mathcal{A}(\Omega_{t+1})=b.$

  2. 2.

The approximate dual feasibility for $(W_{t}^{\star},y_{t}^{\star})$ satisfies

$\displaystyle W_{t}^{\star}\succeq 0,\quad
and\quad\|W_{t}^{\star}-C+\mathcal{A}^{*}(y_{t}^{\star})\|^{2}\leq\frac{2\alpha}{\beta}(F(\Omega_{t})-F(X^{\star})).$

  3. 3.

The approximate primal-dual optimality for
$(\Omega_{t+1},W_{t}^{\star},y_{t}^{\star})$ satisfies

$\displaystyle\langle C,\Omega_{t+1}\rangle-\langle b,y_{t}^{\star}\rangle$
$\displaystyle\geq-\frac{\rho(F(\Omega_{t})-F(X^{\star}))}{\mathcal{D}_{Z_{\star}}+1}-\mathcal{D}_{\Omega_{0}}\sqrt{\frac{2\alpha}{\beta}(F(\Omega_{t})-F(X^{\star}))},$
$\displaystyle\langle C,\Omega_{t+1}\rangle-\langle b,y_{t}^{\star}\rangle$
$\displaystyle\leq\frac{1-\beta}{\beta}(F(\Omega_{t})-F(X^{\star}))+\mathcal{D}_{\Omega_{0}}\sqrt{\frac{2\alpha}{\beta}(F(\Omega_{t})-F(X^{\star}))},$

where $\mathcal{D}_{\Omega_{0}}:=\max_{F(\Omega_{t})\leq
F(\Omega_{0})}\|\Omega_{t}\|$ is bounded due to the compactness of
$\mathcal{P}^{\star}$(see Lemma D.1).

Theorem 4.1 is a direct consequence by combining Lemma 4.1 with the
convergence results of the generic bundle method in Lemma 2.2.

Proof of Theorem 4.1: Consider the generic convex optimization 5. For any
lower approximation function $\hat{f}_{t+1}(\cdot)$ satisfying 9a-9c, Lemma
2.2 guarantees that the generic bundle method in Algorithm 1 generates
iterates $\omega_{t}$ with the gap $f(\omega_{t})-f(x^{\star})$ converging to
zero at a rate of $\mathcal{O}(1/\epsilon^{3})$. This rate improves to
$\mathcal{O}(1/\epsilon)$ whenever the quadratic growth condition 11 occurs.

In our case of SDPs, the quadratic growth of the primal penalized cost
function $F(X)$ in 20 holds whenever strict complementarity holds (see Lemma
D.3): fix $\epsilon>0$ and define a $\epsilon$-sublevel set
$\mathcal{P}_{\epsilon}=\\{X\in\mathbb{S}^{n}\;|\;\mathcal{A}(X)=b,F(X)-F(X^{\star})\leq\epsilon\\},$
then there exist a constant $\mu>0$ such that

$F(X)-F(X^{\star})\geq\mu\cdot\mathrm{dist}^{2}(X,\mathcal{P}^{\star}),\,\;\forall
X\in\mathcal{P}_{\epsilon}.$ (37)

In other words, the square of the distance between a point $X$ and the set of
primal optimal solutions $\mathcal{P}^{\star}$ is bounded by the optimality of
the objective value when strict complementarity holds. Therefore, if the lower
approximation model $\hat{F}_{(\bar{W}_{t},P_{t})}(X)$ in 23 satisfies 9a, 9b
and 9c, the gap $F(\Omega_{t})-F(X^{\star})$ from
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP converges to zero at the rate of
$\mathcal{O}\mathopen{}\mathclose{{}\left(1/\epsilon^{3}}\right)$. This rate
improves to $\mathcal{O}(1/\epsilon)$ whenever strict complementarity holds.

By Lemma 4.1, the convergence results for approximate primal feasibility 35a,
approximate dual feasibility 35b, and primal and dual gap 35c are naturally
established. It remains to verify 9a, 9b and 9c for the lower approximation
$\hat{F}_{(\bar{W}_{t},P_{t})}(X)$, which is provided in Section D.1. This
completes the proof. $\blacksquare$

#### 4.4.2. Proof of Theorem 4.2

Thanks to Theorem 4.1, the iterate $\Omega_{t}$ will be sufficiently close to
the set of optimal solutions $\mathcal{P}^{\star}$ after some finite number of
iterations (which is independent of $\epsilon$). In particular, let ${T}_{0}$
be the number of iterations that ensures

$\min_{X^{\star}\in\mathcal{P}^{\star}}\|\Omega_{t}-X^{\star}\|_{\mathrm{op}}\leq\frac{\delta}{3},$

where $\delta$ is a constant eigenvalue gap parameter (see Lemma D.5 in
Section D.3). Then, the lower approximation model
$\hat{F}_{(\bar{W}_{t},P_{t})}(X)$ in 23 becomes quadratically close to the
true penalized cost function $F(X)$, as summarized in Lemma 4.2. Thanks to the
quadratic closeness of $\hat{F}_{(\bar{W}_{t},P_{t})}(X)$, Algorithm 2 will
take only descent steps after $T_{0}$ iterations, and we have contractions in
terms of the cost value gap and the distance to the set of primal optimal
solutions (Lemma 4.3).

######  Lemma 4.2.

Suppose Assumptions 2.1 and 2.2 are satisfied and strict complementarity holds
for P and D. Let $r_{\mathrm{c}}\geq\sigma_{\mathrm{p}}$. After ${T}_{0}$
iterations, there exists a constant $\eta>0$ (independent of $\epsilon$) such
that

$\hat{F}_{(\bar{W}_{t},P_{t})}(X)\leq
F(X)\leq\hat{F}_{(\bar{W}_{t},P_{t})}(X)+\frac{\eta}{2}\|X-\Omega_{t}\|^{2},\quad\forall
X\in\mathcal{X}_{0}.$

######  Lemma 4.3.

Under the conditions in Lemma 4.2, for any $\alpha\geq\eta$ and
$t\geq{T}_{0}$, Algorithm 2 with $\beta\in(0,\frac{1}{2}]$ takes only descent
steps and guarantees two contractions444Note that [29, Lemma 3.8] only shows
the linear convergence in terms of the distance to the set of optimal
solutions. We here further establish the linear convergence in terms of cost
value gap $F(\Omega_{t})-F(X^{\star})$, which can be directly used in Lemma
4.1.

$\displaystyle\operatorname{dist}(\Omega_{t+1},\mathcal{P}^{\star})$
$\displaystyle\leq\sqrt{\frac{{\alpha}/{2}}{\mu+{\alpha}/{2}}}\operatorname{dist}(\Omega_{t},\mathcal{P}^{\star}),$
(38a) $\displaystyle F(\Omega_{t+1})-F(X^{\star})$
$\displaystyle\leq\mathopen{}\mathclose{{}\left(1-\min\mathopen{}\mathclose{{}\left\\{\frac{\mu}{2\alpha},\frac{1}{2}}\right\\}\beta}\right)\mathopen{}\mathclose{{}\left(F(\Omega_{t})-F(X^{\star})}\right),$
(38b)

where $\mu$ is the quadratic growth constant in 37 for the initial sublevel
set.

Proof of Theorem 4.2: Theorem 4.2 is a direct consequence by combining Lemma
4.3 with Lemma 4.1. After $T_{0}$ iterations, (38b) guarantees the linear
convergence of $F(\Omega_{t})-F(X^{\star})$. Thus, together with Lemma 4.1,
this ensures the linear convergence of approximate primal feasibility,
approximate dual feasibility, and approximate primal-dual optimality after
$T_{0}$ iterations. $\blacksquare$

We finally note that $T_{0}$ is a constant depending on problem data only. The
proofs of Lemma 4.2 and 4.3 are provided in Section D.3.

##  5 Spectral bundle methods for dual SDPs

In this section, we first review the existing spectral bundle method for dual
SDPs that was originally proposed in [19] and further developed in [20, 33,
29]. The (sub)linear convergence results have been recently established in
[29]. We then compare the spectral bundle method for dual SDPs with that for
primal SDPs developed in Section 4.

###  5.1 Spectral bundle algorithms for dual SDPs

The idea in [19] applies the standard bundle method in Section 2.3 to the
penalized dual formulation 13. One key step is to construct an appropriate
lower approximation model. For notational simplicity, let us denote the
objective function in 13 as

$\displaystyle
F_{\mathrm{d}}(y):=-b^{\mathsf{T}}y+\rho\max\mathopen{}\mathclose{{}\left\\{\lambda_{\max}\mathopen{}\mathclose{{}\left(\mathcal{A}^{*}(y)-C}\right),0}\right\\}.$

Similar to 23, the family of spectral bundle methods for dual SDPs uses a
positive semidefinite definite matrix $\bar{W}_{t}\in\mathbb{S}^{n}_{+}$ with
$\mathop{\bf tr}(\bar{W}_{t})=1$ and an orthonormal
$P_{t}\in\mathbb{R}^{n\times r}$ matrix at iteration $t$ to lower approximate
$F_{d}(y)$. Specifically, the lower approximation model is constructed as

$\displaystyle\hat{F}_{\mathrm{d},(\bar{W}_{t},V_{t})}(y)=-b^{\mathsf{T}}y+\rho\max_{S\in\mathbb{S}^{r}_{+},\gamma\geq
0,\gamma+\mathop{\bf tr}(S)\leq
1}\mathopen{}\mathclose{{}\left\langle\gamma\bar{W}_{t}+P_{t}SP_{t}^{\mathsf{T}},\mathcal{A}^{*}(y)-C}\right\rangle.$
(39)

It is clear to see that $\hat{F}_{\mathrm{d},(\bar{W}_{t},V_{t})}(y)\leq
F_{\mathrm{d}}(y),\forall y\in\mathbb{R}^{m}$.

Following [29], we present a family of spectral bundle algorithms for dual
SDPs, called $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD, which has the following
steps:

  * •

Initialization: $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD starts with an initial
guess $\omega_{0}\in\mathbb{R}^{m}$, $P_{0}\in\mathbb{R}^{n\times r}$ formed
by top $r$ normalized eigenvectors of $\mathcal{A}^{*}(\omega_{0})-C$, and
$\bar{W}_{0}\in\mathbb{S}^{n}$ with $\mathop{\bf tr}(\bar{W}_{0})=1$. The
matrices $P_{0}$ and $\bar{W}_{0}$ are used to construct an initial lower
approximation model $\hat{F}_{\mathrm{d},(\bar{W}_{t},V_{t})}(y)$ in 39.

  * •

Solve the master problem: Similar to 7, at iteration $t\geq 0$,
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD solves the following regularized master
problem

$\displaystyle\mathopen{}\mathclose{{}\left(y_{t+1}^{\star},S_{t}^{\star},\gamma^{\star}_{t}}\right)=$
$\displaystyle\operatorname*{\arg\\!\min}_{y\in\mathbb{R}^{n}}\;\;\hat{F}_{\mathrm{d},(\bar{W}_{t},P_{t})}(y)+\frac{\alpha}{2}\|y-\omega_{t}\|^{2},$
(40)

where $\omega_{t}\in\mathbb{R}^{m}$ is the current reference position, and
$\alpha>0$ is as a penalty parameter.

  * •

Update reference point: Similar to 8, $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD
updates the next reference point $\omega_{t+1}$ as follows: given
$\beta\in(0,1)$, if

$\displaystyle\beta\mathopen{}\mathclose{{}\left(F_{\mathrm{d}}(\omega_{t})-\hat{F}_{\mathrm{d},(\bar{W}_{t},P_{t})}(y_{t+1}^{\star})}\right)\leq
F_{\mathrm{d}}(\omega_{t})-F_{\mathrm{d}}(y_{t+1}^{\star})$ (41)

holds true, we let $\omega_{t+1}=y_{t+1}^{\star}$ (descent step), otherwise
$\omega_{t+1}=\omega_{t}$ (null step).

  * •

Update the lower approximation model: $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD
updates the matrices $\bar{W}_{t+1},P_{t+1}$ for the lower approximation model
$\hat{F}_{\mathrm{d},(\bar{W}_{t},V_{t})}(y)$ using formulas in 26 and 27 with
one difference that the matrix $V_{t}\in\mathbb{R}^{n\times r_{\mathrm{c}}}$
is formed by the top $r_{\mathrm{c}}\geq 1$ eigenvectors of
$\mathcal{A}^{*}(y_{t+1}^{\star})-C$.

The overall process for $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD is listed in
Algorithm 3. $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD generates a sequence of
iterates $\\{\omega_{t},W_{t}^{\star}\\}$ with monotonically decreasing cost
values $\\{F_{d}(\omega_{t})\\}$. Solving the subproblem 40 is the main
computation in each iteration of $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD.
Similar to Proposition 4.1, we have the following result.

######  Proposition 5.1 ([29, Section 2.3]).

The master problem 40 is equivalent to

$\displaystyle\min_{W\in\hat{\mathcal{W}}_{t}}\;\;\langle
b,\omega_{t}\rangle+\langle
W,C-\mathcal{A}^{*}(\omega_{t})\rangle+\frac{1}{2\alpha}\mathopen{}\mathclose{{}\left\|b-\mathcal{A}W}\right\|^{2},$
(42)

where the constraint set is defined as

$\hat{\mathcal{W}}_{t}:=\\{\gamma\bar{W}_{t}+P_{t}SP_{t}^{\mathsf{T}}\in\mathbb{S}^{n}\;|\;S\in\mathbb{S}^{r}_{+},\gamma\geq
0,\gamma+\mathop{\bf tr}(S)\leq\rho\\}.$ (43)

The optimal $y$ in 40 is recovered as
$y_{t+1}^{\star}=\omega_{t}+\frac{1}{\alpha}\bigl{(}b-\mathcal{A}(W_{t}^{\star})\bigr{)},$
where $W_{t}^{\star}$ is a minimizer of 42.

Similar to Remark 4.1, problem 42 is a quadratic problem with a small
semidefinite constraint, which can be reformulated into a problem of the form
34; see Section C.2. The convergence results of
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD are summarized in Theorems 5.1 and 5.2.

Algorithm 3 $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD– Spectral Bundle Method for
Dual SDPs

1:Problem data $A_{1},\ldots,A_{m},C\in\mathbb{S}^{n}$, $b\in\mathbb{R}^{n}$.

2:Parameters $r_{\mathrm{p}}\geq 0,r_{\mathrm{c}}\geq 1$, $\alpha>0$,
$\beta\in(0,1)$, $\epsilon\geq 0$. An initial point
$\omega_{0}\in\mathbb{R}^{m}$.

3:Initialization:  Let ${r}=r_{\mathrm{p}}+r_{\mathrm{c}}$. Initialize
$\bar{W}_{0}\in\mathbb{S}^{n}_{+}$ with $\mathop{\bf tr}(\bar{W}_{0})=1$.
Compute $P_{0}\in\mathbb{R}^{n\times{r}}$ with columns being the top ${r}$
orthonormal eigenvectors of $\mathcal{A}^{*}(\omega_{0})-C$.

4:for $t=0,\ldots,t_{\max}$ do

5: Solve 40 to obtain $y_{t+1}^{\star},\gamma^{\star}_{t},S^{\star}_{t}$.
$\triangleright$ master problem

6: Form the iterate $W_{t}^{\star}$ as in 43.

7: if 41 holds then

8: Set primal iterate $\omega_{t+1}=y_{t+1}^{\star}$. $\triangleright$ descent
step

9: else

10: Set primal iterate $\omega_{t+1}=\omega_{t}$. $\triangleright$ null step

11: end if

12: Compute $P_{t+1}$ as 26, and $\bar{W}_{t+1}$ as 27. $\triangleright$
update model

13: if stopping criterion are met then

14: Quit. $\triangleright$ termination

15: end if

16:end for

######  Theorem 5.1 ([29, Theorem 3.1]).

Suppose Assumptions 2.1 and 2.2 are satisfied. Given any $\beta\in(0,1)$,
$r_{\mathrm{c}}\geq 1$, $r_{\mathrm{p}}\geq 0$, $\alpha>0$,
$r=r_{\mathrm{c}}+r_{\mathrm{p}}$, $\rho>2\mathcal{D}_{X_{\star}}$,
$P_{0}\in\mathbb{R}^{n\times r}$, $\omega_{0}\in\mathbb{R}^{m}$, and accuracy
$\epsilon>0$, the $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD in Algorithm 3
produces iterates $\omega_{t}$ and $W_{t}^{\star}$ with
$F_{\mathrm{d}}(\omega_{t})-F_{\mathrm{d}}(y^{\star})\leq\epsilon$ and

approximate primal feasibility:
$\displaystyle\quad\|b-\mathcal{A}(W_{t}^{\star})\|^{2}\leq\epsilon,\quad
W_{t}^{\star}\succeq 0,$ (44a) approximate dual feasibility:
$\displaystyle\quad\lambda_{\min}(C-\mathcal{A}^{*}(\omega_{t}))\geq-\epsilon,$
(44b) approximate primal-dual optimality: $\displaystyle\quad|\langle
C,W_{t}^{\star}\rangle-\langle b,\omega_{t}\rangle|\leq\sqrt{\epsilon}$ (44c)

by some iteration $t\leq\mathcal{O}(1/\epsilon^{3})$. If additionally, strict
complementarity holds, then these conditions are reached by some iteration
$t\leq\mathcal{O}(1/\epsilon)$.

Along with strict complementarity, a further improvement on the convergence
rate can be shown if the number of selected current eigenvectors at every
iteration satisfies

$r_{\mathrm{c}}\geq\sigma_{\mathrm{d}}:=\max_{(y^{\star},Z^{\star})\in\mathcal{D}^{\star}}\dim(\text{null}(Z^{\star})).$
(45)

######  Theorem 5.2 ([29, Theorem 3.2]).

Suppose Assumptions 2.1 and 2.2 are satisfied and strict complementarity holds
for P and D. There exist constants $T_{0}>0$ and $\eta>0$. Under a proper
selection of $\alpha$ and any $\beta\in(0,\frac{1}{2}]$,
$r_{\mathrm{c}}\geq\sigma_{\mathrm{d}}$, $r_{\mathrm{p}}\geq 0$,
$r=r_{\mathrm{c}}+r_{\mathrm{p}}$, $\rho>2\mathcal{D}_{X_{\star}}$,
$P_{0}\in\mathbb{R}^{n\times r}$, $\omega_{0}\in\mathbb{R}^{m}$, and accuracy
$\epsilon>0$, after at most ${T}_{0}$ iterations, the
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD in Algorithm 3 only takes descent steps
and converges linearly to an optimal solution. Consequently,
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD produces iterates $\omega_{t}$ and
$W_{t}^{\star}$ with
$F_{\mathrm{d}}(\omega_{t})-F_{\mathrm{d}}(y^{\star})\leq\epsilon$ and 44 by
at most $T_{0}+\mathcal{O}(\log(1/\epsilon))$ iterations.

The constants $T_{0}>0$ and $\eta>0$ only depend on problem data and are
independent of $\epsilon$. We refer the interested reader to [29, Section
3.3.1] for details.

######  Remark 5.1.

Historically, the spectral bundle method was first introduced in [19] to
tackle SDP relaxations for large-scale combinatorial problems. The method in
[19] works for dual SDPs with an explicit constant trace constraint; see 18.
The algorithm [19] requires one current eigenvalue, i.e. $r_{\mathrm{c}}=1$,
and allows different values of $r_{\mathrm{p}}$. This selection of parameters
is different from the selection of parameters reviewed in Algorithm 3 (we
follow the choice in [29]). The convergence results in Theorem 5.2 show that
the parameter $r_{\mathrm{p}}$ can be chosen as zero which still guarantees
the linear convergence whenever $r_{\mathrm{c}}\geq\sigma_{\mathrm{d}}$,
$\alpha$ is chosen correctly, and the SDP satisfies strict complementarity.
The result is based on one critical observation that
$\gamma\bar{W}_{t}+P_{t}SP_{t}^{\mathsf{T}}$ in 39 is an approximation of the
optimal primal variable $X^{\star}$ and $P_{t}$ is an approximation of the
null space of $Z^{\star}$. A key result in the analysis is based on a novel
eigenvalue approximation in [29, Lemma 3.9]. One immediate benefit of
selecting $r_{\mathrm{c}}$ instead of $r_{\mathrm{p}}$ is that
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD becomes particularly efficient when the
primal SDP admits low-rank solutions as $r_{\mathrm{c}}$ can be chosen small
to fulfill the rank condition 45. In this case, the master problem 40 can be
solved efficiently. The low-rank property indeed holds for many SDPs from
combinatorial problems and phase retrieval [55].

###  5.2 Comparison and connections

In this subsection, we compare the differences and draw the connections
between primal and dual formulations of the spectral bundle method.

Table 1: Comparison between $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP in
Algorithm 2 and $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD in Algorithm 3.
$W_{t}^{\star}$ denotes the optimizer of the master problem in constraint set
$\hat{\mathcal{W}}_{t}$ in 30 for $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP (or
43 for $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD). The first block of rows
(Comp.) denotes computational details, and the second block of rows (Conv.)
shows convergence details.

| Descriptions |  $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP |  $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD
---|---|---|---
Comp. | Domain | $\mathcal{X}_{0}=\\{X\in\mathbb{S}^{n}\;|\;\mathcal{A}(X)=b\\}$ | $\mathbb{R}^{m}$
Reference, candidate | $\Omega_{t}\in\mathbb{S}^{n},X_{t+1}^{\star}\in\mathbb{S}^{n}$ | $\omega_{t}\in\mathbb{R}^{m},y_{t+1}^{\star}\in\mathbb{R}^{m}$
Penalized parameter | $\rho>\sup_{(y^{\star},Z^{\star})\in\mathcal{D}^{\star}}\;{\mathop{\bf tr}(Z^{\star})}$ | $\rho>\sup_{X^{\star}\in\mathcal{P}^{\star}}\;{\mathop{\bf tr}(X^{\star})}$
Nonsmooth objective | $\langle C,X\rangle+\rho\max\\{\lambda_{\max}(-X),0\\}$ | $-b^{\mathsf{T}}y+\rho\max\mathopen{}\mathclose{{}\left\\{\lambda_{\max}\mathopen{}\mathclose{{}\left(\mathcal{A}^{*}(y)-C}\right),0}\right\\}$
Approximation $(\hat{f})$ | $\langle C,X\rangle+\max_{W\in\hat{\mathcal{W}}_{t}}\langle W,-X\rangle$ | $-b^{\mathsf{T}}y+\max_{W\in\hat{\mathcal{W}}_{t}}\mathopen{}\mathclose{{}\left\langle W,\mathcal{A}^{*}(y)-C}\right\rangle$
Master problem | $\min_{X\in\mathcal{X}_{0}}\;\hat{f}(X)+\frac{\alpha}{2}\|X-\Omega_{t}\|^{2}$ | $\min_{y\in\mathbb{R}^{m}}\;\hat{f}(y)+\frac{\alpha}{2}\|y-\omega_{t}\|^{2}$
Current information | Top $r_{\mathrm{c}}$ eigenvectors of $-X_{t+1}^{\star}$ | Top $r_{\mathrm{c}}$ eigenvectors of $\mathcal{A}^{*}\mathopen{}\mathclose{{}\left(y_{t+1}^{\star}}\right)-C$
Conv. | Primal feasibility | $\mathcal{A}(\Omega_{t})-b=0,\lambda_{\min}(\Omega_{t})\geq-\epsilon$ | $\|\mathcal{A}(W_{t}^{\star})-b\|^{2}\leq\epsilon,W_{t}^{\star}\succeq 0$
Dual feasibility | $\mathopen{}\mathclose{{}\left\|W_{t}^{\star}+\mathcal{A}^{*}(y_{t}^{\star})-C}\right\|^{2}\leq\epsilon,W_{t}^{\star}\succeq 0$ | $\lambda_{\min}\mathopen{}\mathclose{{}\left(\mathcal{A}^{*}\mathopen{}\mathclose{{}\left(\omega_{t}}\right)-C}\right)\geq-\epsilon$
Duality gap | $\lvert\langle C,\Omega_{t}\rangle-\langle b,y_{t}^{\star}\rangle\rvert\leq\sqrt{\epsilon}$ | $\lvert\langle C,W_{t}^{\star}\rangle-\langle b,\omega_{t}\rangle\rvert\leq\sqrt{\epsilon}$
Bound on $r_{\mathrm{c}}$ | $\max_{X^{\star}\in\mathcal{P}^{\star}}\dim(\text{null}(X^{\star}))$ | $\max_{(y^{\star},Z^{\star})\in\mathcal{D}^{\star}}\dim(\text{null}(Z^{\star}))$

It is clear that both $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP in Algorithm 2
(primal SDPs) and $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD in Algorithm 3 (dual
SDPs) follow the same framework of the generic bundle method (see Algorithm
1). One main difference is that $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP solves
a constrained nonsmooth problem 12 while
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD deals with an unconstrained nonsmooth
problem 13. The unconstrained nonsmooth problem 13 is in the form of
eigenvalue minimization which has extensive literature [49]. Table 1 presents
a comparison of the computational details and convergence results between
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP and
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD. In particular, we have the following
observations:

  * •

$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP generates iterates
$(\Omega_{t+1},W_{t}^{\star},y_{t}^{\star})$ that 1) exactly satisfy the
primal affine constraint $\mathcal{A}(\Omega_{t+1})=b$ and dual PSD constraint
$\lambda_{\min}(W_{t}^{\star})\geq 0$, but 2) do not exactly satisfy the
primal PSD constraint $\lambda_{\min}(\Omega_{t+1})\geq-\epsilon$ and the dual
affine constraint
$\mathopen{}\mathclose{{}\left\|W_{t}^{\star}+\mathcal{A}^{*}(y_{t}^{\star})-C}\right\|^{2}\leq\epsilon$.

  * •

$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD outputs iterates
$(W_{t}^{\star},y_{t}^{\star})$ that 1) exactly satisfy the dual affine
constraint (by construction) and the primal PSD constraint
$(\lambda_{\min}(W_{t}^{\star})\geq 0)$, but 2) do not exactly satisfy the
primal affine constraint $\|\mathcal{A}(W_{t}^{\star})-b\|^{2}\leq\epsilon$
and the dual PSD constraint
$\lambda_{\min}(C-\mathcal{A}^{*}(\omega_{t+1}))\geq-\epsilon$.

  * •

Each iteration of $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP and
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD requires solving a small quadratic SDP
in the form of 34. When choosing the same values for $r_{\mathrm{p}}$ and
$r_{\mathrm{c}}$, the subproblems from $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP
and $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD have the same dimension, and thus
the computational complexity is very similar (although their constructions of
34 are slightly different; see Appendix C).

Both $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP and
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD have very similar sublinear and linear
convergence rates, as shown in Theorems 4.1 and 4.2 and Theorems 5.1 and 5.2,
respectively. However, we here highlight a key difference in the condition for
linear convergence. When constructing the lower approximation model at each
iteration, $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP selects the top
$r_{\mathrm{c}}$ eigenvectors of $-X^{\star}_{t+1}$ to approximate the null
space of the primal optimal solutions $X^{\star}$, while
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD uses $r_{\mathrm{c}}$ eigenvectors of
$\mathcal{A}^{*}(\omega_{t+1})-C$ to approximate the null space of the dual
optimal solutions $Z^{\star}$. If the null space of the primal optimal
solutions in $\mathcal{P}^{\star}$ has a smaller dimension than that of the
dual optimal solutions in $\mathcal{D}^{\star}$, i.e.,

$\displaystyle\max_{X^{\star}\in\mathcal{P}^{\star}}\dim(\text{null}(X^{\star}))\ll\max_{(y^{\star},Z^{\star})\in\mathcal{D}^{\star}}\dim(\text{null}(Z^{\star})),$
(46a) then the master problem in $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP can
select a smaller number of eigenvectors, leading to a smaller PSD constraint
in 34. This is reflected in the linear convergence condition on $r_{c}$; see
36 and 45. In this case, it will be more numerically beneficial to apply
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP that solves the primal SDPs directly as
it is easier to solve the subproblem 34 at each iteration. If, on the other
hand, we have the following relationship
$\displaystyle\max_{X^{\star}\in\mathcal{P}^{\star}}\dim(\text{null}(X^{\star}))\gg\max_{(y^{\star},Z^{\star})\in\mathcal{D}^{\star}}\dim(\text{null}(Z^{\star})),$
(46b) then it will be more numerically beneficial to apply
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD that solves the dual SDPs directly.

######  Remark 5.2 (Low-rank primal solutions versus low-rank dual solutions).

Note that 46a implies that the dual SDP D admits low-rank optimal solutions
(i.e., $\mathrm{rank}(Z^{\star})$ is small), while 46b indicates that the
primal SDP P has low-rank optimal solutions (i.e., $\mathrm{rank}(X^{\star})$
is small). Therefore, it is crucial to choose the appropriate algorithm to
solve the problem if we have prior knowledge of the rank property of SDPs. For
instance, SDP-based optimization problems from sum-of-square relaxation and
their applications [35, 36] are likely to admit low-rank dual optimal
solutions in D. On the other hand, SDP relaxation of combinatorial problems
such as Max-Cut [17], matrix completion [52], and phase retrieval [55] are
likely to admit low-rank primal solutions in P. We provide numerical
experiments of SOS optimization and Max-Cut in Section 6, which indeed
validate the effect of the rank property on the convergence behavior of
Algorithms 2 and 3. We finally note that the low-rank properties also depend
on how SDPs are formulated in different applications; see the conversion
between primal and dual SDPs in Section C.4.

##  6 Implementation and numerical experiments

We have implemented $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP and
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD in Algorithms 2 and 3 in an open-source
MATLAB package, which is available at

https://github.com/soc-ucsd/specBM.

In this section, we discuss some implementation details and present three sets
of numerical experiments to validate the performance of
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP and
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD, especially their linear convergence
behaviors. The numerical results confirm the discussions in Section 5.2:
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP works better when the dual SDP D admits
low-rank solutions (i.e., 46a holds). Similarly,
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD is more beneficial when the primal SDP
P admits low-rank solutions (i.e., 46b holds).

In Section 6.2, similar to [29, Section 5.1], we consider SDPs with randomly
generated problem data, which demonstrates that
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP and
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD admit sublinear convergence under
different configurations of $(r_{\mathrm{p}},r_{\mathrm{c}})$ and that the
algorithms have a linear convergence when the rank conditions in 36 and 45
hold. In Section 6.3, we consider a benchmark combinatorial problem – Max-Cut
and its standard SDP relaxation, which is likely to admit low-rank primal
solutions. In this case, it is more desirable to apply
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD than
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP. In Section 6.4, we show that
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP is well-suited for SDP relaxation
rising from sum-of-squares (SOS) optimization, which appears to admit low-rank
dual solutions. In this case, a small number of current eigenvectors
$r_{\mathrm{c}}$ is sufficient to fulfill the rank condition 36 that ensures
fast linear convergence. To highlight the efficiency of
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP, we further compare it with the state-
of-the-art interior-point and first-order SDP solvers: we choose SDPT3 [40]
and MOSEK [41] as the interior-point solvers, and we select CDCS [42] and
SDPNAL+ [16] as the first-order solvers.

All numerical experiments are conducted on a PC with a 12-core Intel i7-12700K
<EMAIL_ADDRESS>and 32GB RAM.

###  6.1 Implementation details

One major computation in each iteration of
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP and
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD is to solve the master problems 29 and
42. As discussed in Remark 4.1, we have implemented an automatic
transformation from 29 and 42 into standard conic form 34. Then, a subroutine
is required to solve 34. In our current implementation, we used MOSEK [41] to
get an exact solution of 34 when $r=r_{\mathrm{p}}+r_{\mathrm{c}}>1$ and
implemented the analytical solution in Section C.3 when
$r_{\mathrm{p}}=0,r_{\mathrm{c}}=1$. We note that customized algorithms can be
developed to 34 with varying accuracy which will further improve numerical
efficiency at each iteration.

To form the lower approximation models 23 and 39, we computed eigenvalues and
eigenvectors using the routine eig in MATLAB. Faster eigenvalue/eigenvector
computations can be implemented using dsyevx in LAPACK, similar to [16,
Section 5]. The orthogonalization process to update the matrix $P_{t+1}$ in 26
was implemented using the routine orth in MATLAB.

#### 6.1.1. Adaptive strategy on the regularization parameter.

Algorithms 2 and 3 are guaranteed to converge with any regularization
parameter $\alpha>0$ in the subproblems 24 and 40. Yet, the value of $\alpha$
largely influences the practical convergence performance, as highlighted in
[56] and [38, Page 380]. In our implementation of Algorithm 2, the parameter
$\alpha$ uses the adaptive updating rule below555The implementation of
Algorithm 3 uses the same updating rule and replaces $\Omega_{t}$,
$X_{t+1}^{\star}$, $F$, and $\hat{F}_{(\bar{W}_{t},P_{t})}$ with $\omega_{t}$,
$y_{t+1}^{\star}$, $F_{\mathrm{d}}$ and
$\hat{F}_{\mathrm{d},(\bar{W}_{t},P_{t})}$ respectively.:

$\displaystyle\alpha=\begin{cases}\min\\{2\alpha,\alpha_{\max}\\},\leavevmode\nobreak\
\text{if
}m_{\mathrm{l}}\mathopen{}\mathclose{{}\left(F(\Omega_{t})-\hat{F}_{(\bar{W}_{t},P_{t})}(X_{t+1}^{\star})}\right)\geq
F(\Omega_{t})-F(X_{t+1}^{\star})\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ N_{\mathrm{c}}\geq N_{\min}\\\
\max\\{\alpha/2,\alpha_{\min}\\},\leavevmode\nobreak\ \text{if
}m_{\mathrm{r}}\mathopen{}\mathclose{{}\left(F(\Omega_{t})-\hat{F}_{(\bar{W}_{t},P_{t})}(X_{t+1}^{\star})}\right)\leq
F(\Omega_{t})-F(X_{t+1}^{\star}),\end{cases}$

where $m_{\mathrm{r}}>\beta$ and $0<m_{\mathrm{l}}<\beta$ are two nonnegative
parameters that indicate the effectiveness of the current approximation model
$\hat{F}_{(\bar{W}_{t},P_{t})}$ and the candidate point $X_{t+1}^{\star}$,
$\alpha_{\min}$ and $\alpha_{\max}$ are two nonnegative parameters that keep
$\alpha$ staying in the interval $[\alpha_{\min},\alpha_{\max}]$,
$N_{\mathrm{c}}$ counts the number of consecutive null steps, and $N_{\min}$
is the threshold that controls the frequency of increasing $\alpha$.

In our implementation for Algorithms 2 and 3, the default parameters are
chosen as
$m_{\mathrm{l}}=0.001,\alpha_{\min}=10^{-5},\alpha_{\max}=100,\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ N_{\min}=10.$ The parameter $\beta$ and
$m_{\mathrm{r}}$ were tuned slightly for different classes of instances in our
experiments. The initial points are chosen as $\Omega_{0}=I$ and
$\omega_{0}=0$ for Algorithm 2 and Algorithm 3 respectively.

#### 6.1.2. Suboptimality measures.

For the pair of primal and dual SDPs P and D, we measure the feasibility and
optimality of a candidate solution
$(X,y,Z)\in\mathbb{S}^{n}\times\mathbb{R}^{m}\times\mathbb{S}^{n}$ using

$\displaystyle\eta_{1}$
$\displaystyle=\frac{\|\mathcal{A}(X)-b\|}{1+\|b\|},\qquad\quad\eta_{2}=\min\\{0,\lambda_{\min}(X)\\},$
(47a) $\displaystyle\eta_{3}$
$\displaystyle=\frac{\|\mathcal{A}^{*}(y)+Z-C\|}{1+\|C\|},\quad\eta_{4}=\min\\{0,\lambda_{\min}(Z)\\},$
(47b) $\displaystyle\eta_{5}$ $\displaystyle=\frac{|\langle
C,X\rangle-b^{\mathsf{T}}y|}{1+|\langle C,X\rangle|+|b^{\mathsf{T}}y|}.$ (47c)

In 47a, $\eta_{1}$ and $\eta_{2}$ measure the violation of the affine
constraint and the conic constraint in the primal SDP, respectively. The
measures $\eta_{3}$/$\eta_{4}$ in 47b quantify the violation of the
affine/conic constraints in the dual SDP. The last index $\eta_{5}$ in 47c
measures the duality gap.

As stated in Theorem 4.1, all iterates $\Omega_{t}$ (primal variable),
$y_{t}^{\star}$ and $W_{t}^{\star}$ (dual variables) from
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP ensure $\eta_{1}=0$ and $\eta_{4}=0$
(up to machine accuracy). On the other hand, the iterates $W_{t}^{\star}$
(primal variable) and $\omega_{t}$ (dual variable) from
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD guarantee $\eta_{2}=0$ and $\eta_{3}=0$
(see Theorem 5.1). In Sections 6.2 and 6.3, we ran
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP and
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD for a fixed number of iterations and
report the measures $\eta_{1},\ldots,\eta_{5}$ for the final iterate. In
Section 6.4, for a given tolerance $\textsf{tol}\geq 0$, we terminate the
algorithms when

$\max\\{\eta_{1},\ldots,\eta_{5}\\}\leq\textsf{tol}.$

The performance of $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP was compared with
the baseline solvers SDPT3 [40], MOSEK [41], CDCS [42] and SDPNAL+ [16] in our
last experiment. In Sections 6.2 and 6.3, $f(\omega_{t})$ denotes the value of
the generic cost function 5, and it refers to the cost function in 12 for
primal SDPs and in 13 for dual SDPs.

###  6.2 SDPs with randomly generated problem data

Our first experiment demonstrates the (sub)linear convergence of
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP and
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD under different configurations of
$(r_{\mathrm{p}},r_{\mathrm{c}})$. Similar to [29, Section 5.1], we randomly
generated two SDPs, satisfying strict complementarity, in the form of P and D.
Both SDP instances have a PSD constraint of dimension $n=1000$ and affine
constraints of size $m=200$. The first SDP admits a low-rank dual solution
($\text{rank}(Z^{\star})=3$) and the second SDP admits a low-rank primal
solution ($\text{rank}(X^{\star})=3$). Details of generating these SDP
instances are discussed in Section E.1.

We consider two different configurations of the parameters $r_{\mathrm{p}}$
and $r_{\mathrm{c}}$:

  * •

$r_{\mathrm{p}}=0$, while changing
$r_{\mathrm{c}}=r^{\star}-1,\leavevmode\nobreak\
r_{\mathrm{c}}=r^{\star},\leavevmode\nobreak\ \text{and
}r_{\mathrm{c}}=r^{\star}+1,$

  * •

$r_{\mathrm{c}}=1$, while changing
$r_{\mathrm{p}}=r^{\star}-2,\leavevmode\nobreak\
r_{\mathrm{p}}=r^{\star}-1,\leavevmode\nobreak\ \text{and
}r_{\mathrm{p}}=r^{\star},$

where we set $r^{\star}=3$ for both SDP instances (thus we have
$r^{\star}=\text{rank}(Z^{\star})$ or $r^{\star}=\text{rank}(X^{\star})$). In
the first setting, we do not consider any past information but only rely on
the current information $r_{\mathrm{c}}$, while in the second setting, we keep
the minimum amount of current information $r_{\mathrm{c}}$ and rely on the
accumulated past information $r_{\mathrm{p}}$. As the discussion in Section
5.2, $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP and
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD iteratively approximate two different
spaces: the null space of the primal solutions $X^{\star}$ and the null space
of the dual solutions $Z^{\star}$, respectively. When an SDP admits low-rank
dual solutions (i.e., the null space of the primal solutions $X^{\star}$ has a
low dimension), it is computationally more beneficial to choose
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP to solve the SDP. On the other hand,
when an SDP admits low-rank primal solutions (i.e., the null space of the dual
solutions $Z^{\star}$ has a low dimension), it is computationally more
beneficial to choose $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD.

We use the first SDP instance with a low-rank dual solution to demonstrate the
benefits of $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP and its fast convergence
guarantees in Theorem 4.2. For this SDP instance, we ran
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP for both settings and
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD for only the first setting. Then, we
use the second SDP instance with a low-rank primal solution to highlight the
benefits of $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD and validate its fast
convergence guarantees in Theorem 5.2. For the second SDP instance with a low-
rank primal solution, we ran $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD for both
settings and $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP for only the first
setting. The penalty parameter $\rho$ is set $2\mathop{\bf tr}(Z^{\star})+2$
and $2\mathop{\bf tr}(X^{\star})+2$ for $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP
and $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD respectively. The step-size
parameters $\beta$ and $m_{\mathrm{r}}$ are chosen as $0.4$ and $0.7$
respectively.

(a) An SDP instance with $\mathrm{rank}(Z^{\star})=3$

(b) An SDP instance with $\mathrm{rank}(X^{\star})=3$

Figure 1: The relative optimality gap of different choices of $r_{\mathrm{p}}$ and $r_{\mathrm{c}}$ in $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP and $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD for solving two SDPs with $X\in\mathbb{S}^{1000}_{+}$: the left SDP instance admits a low-rank dual solution $\mathrm{rank}(Z^{\star})=3$ which suits well for $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP, while the right SDP instance admits a low-rank primal solution $\mathrm{rank}(X^{\star})=3$ that suits well for $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD.  Table 2: Computational results of $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP and$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD for solving two random SDPs ($n=1000$ and $m=200$) under different values of $r_{\mathrm{p}}$ and $r_{\mathrm{c}}$. “Semi Feasi.” denotes the violation of PSD constraints $\max\\{\eta_{2},\eta_{4}\\}$, “Affine Feasi.” denotes the violation of affine constraints $\max\\{\eta_{1},\eta_{3}\\}$, “Dual Gap” denotes the duality gap $\eta_{5}$ in 47, and “Cost Opt.” denotes the cost-value gap $\eta_{\mathrm{c}}={(f(\omega_{t})-f^{\star})}/{f^{\star}}$. SDP instance | Algorithm | $(r_{\mathrm{p}},r_{\mathrm{c}})$ | Semi Feasi. | Affine Feasi. | Dual Gap | Cost Opt.
---|---|---|---|---|---|---
$\begin{aligned} &\text{Instance 1}\\\ &\text{rank}(Z^{\star})=3\end{aligned}$ |  $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD | $(0,2)$ | $-7.53\mathrm{e}{-1}$ | $1.55\mathrm{e}{-4}$ | $4.28\mathrm{e}{-5}$ | $13.1$
$(0,3)$ | $-7.59\mathrm{e}{-1}$ | $1.54\mathrm{e}{-4}$ | $3.62\mathrm{e}{-4}$ | $13.1$
$(0,4)$ | $-9.96\mathrm{e}{-1}$ | $3.31\mathrm{e}{-4}$ | $1.10\mathrm{e}{-4}$ | $28.6$
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP | $(1,1)$ | $-3.67\mathrm{e}{-3}$ | $3.27\mathrm{e}{-4}$ | $1.24\mathrm{e}{-5}$ | $1.38\mathrm{e}{-4}$
$(2,1)$ | $-2.92\mathrm{e}{-3}$ | $3.18\mathrm{e}{-5}$ | $1.03\mathrm{e}{-5}$ | $1.02\mathrm{e}{-4}$
$(3,1)$ | $-7.22\mathrm{e}{-3}$ | $9.00\mathrm{e}{-6}$ | $9.74\mathrm{e}{-7}$ | $2.17\mathrm{e}{-4}$
$(0,2)$ | $-3.58\mathrm{e}{-4}$ | $2.46\mathrm{e}{-4}$ | $5.06\mathrm{e}{-5}$ | $3.22\mathrm{e}{-5}$
$(0,3)$ | $-1.17\mathrm{e}{-8}$ | $1.57\mathrm{e}{-6}$ | $3.38\mathrm{e}{-7}$ | $2.14\mathrm{e}{-9}$
$(0,4)$ | $-1.40\mathrm{e}{-7}$ | $3.22\mathrm{e}{-6}$ | $6.91\mathrm{e}{-7}$ | $1.05\mathrm{e}{-8}$
$\begin{aligned} &\text{Instance 2}\\\ &\text{rank}(X^{\star})=3\end{aligned}$ |  $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP | $(0,2)$ | $-144$ | $6.82\mathrm{e}{-1}$ | $9.99\mathrm{e}{-1}$ | $4801$
$(0,3)$ | $-121$ | $5.61\mathrm{e}{-1}$ | $9.99\mathrm{e}{-1}$ | $4041$
$(0,4)$ | $-106$ | $4.58\mathrm{e}{-1}$ | $9.99\mathrm{e}{-1}$ | $3369$
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD | $(1,1)$ | $-1.80$ | $3.11\mathrm{e}{-2}$ | $8.21\mathrm{e}{-3}$ | $6.34\mathrm{e}{-1}$
$(2,1)$ | $-5.11\mathrm{e}{-1}$ | $2.00\mathrm{e}{-2}$ | $4.78\mathrm{e}{-3}$ | $2.06\mathrm{e}{-1}$
$(3,1)$ | $-4.95\mathrm{e}{-1}$ | $4.25\mathrm{e}{-3}$ | $4.62\mathrm{e}{-3}$ | $17.7$
$(0,2)$ | $-3.04\mathrm{e}{-1}$ | $5.47\mathrm{e}{-3}$ | $6.74\mathrm{e}{-4}$ | $1.03\mathrm{e}{-1}$
$(0,3)$ | $-4.66\mathrm{e}{-6}$ | $6.65\mathrm{e}{-5}$ | $4.91\mathrm{e}{-5}$ | $4.74\mathrm{e}{-6}$
$(0,4)$ | $-9.30\mathrm{e}{-6}$ | $6.65\mathrm{e}{-5}$ | $1.32\mathrm{e}{-5}$ | $3.94\mathrm{e}{-6}$

In all cases, we ran $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP and
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD for 300 iterations. In this experiment,
we also computed the cost value gap as
$\eta_{\mathrm{c}}={(f(\omega_{t})-f^{\star})}/{f^{\star}}.$ The convergence
behaviors of the cost value gap are illustrated in Figure 1. In Table 2, we
list the suboptimality measures for the final iterates, where “Semi Feasi.”
denotes the violation of PSD constraints $\max\\{\eta_{2},\eta_{4}\\}$,
“Affine Feasi.” denotes the violation of affine constraints
$\max\\{\eta_{1},\eta_{3}\\}$, “Dual Gap” denotes the duality gap $\eta_{5}$
in 47, and “Cost Opt.” denotes the cost-value gap $\eta_{\mathrm{c}}$. As
expected, the value of $r_{\mathrm{c}}$ greatly affects the convergence
performance for both $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP and
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD. In the first SDP with a low-rank dual
solution, we observe that $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP has a fast
convergence behavior when choosing $r_{\mathrm{c}}\geq
3=\text{dim}(\text{null}(X^{\star}))$, while
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD has a slow performance in all settings.
On the other hand, in the second SDP with a low-rank primal solution,
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD enjoys the fast convergence when
$r_{\mathrm{c}}\geq 3=\text{dim}(\text{null}(Z^{\star}))$, while
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP converges poorly in all different
configurations. The numerical results confirm the theoretical convergence
results in Theorems 4.2 and 5.2 and our discussions in Section 5.2.

###  6.3 Max-Cut

(a) Max-Cut problem with Graph G1

(b) Max-Cut problem with Graph G25

Figure 2: The relative optimality gap of different choices of $r_{\mathrm{p}}$
and $r_{\mathrm{c}}$ in $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP and
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD for solving SDP relaxations of two Max-
Cut problems.

In this experiment, we consider the maximum cut problem, which is a benchmark
combinatorial optimization problem. The SDP relaxation is likely to have low-
rank primal solutions, for which $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD is
better suited. Consider an undirected graph
$\mathcal{G}(\mathcal{V},\mathcal{E})$ defined by a set of vertices
$\mathcal{V}=\\{1,2,\ldots,n\\}$ and a set of edges
$\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$, and each edge
$\\{i,j\\}\in\mathcal{E}$ has a weight $w_{ij}=w_{ji}$. The max-cut problem
aims to find a maximum cut that separates the vertices into two different
groups. This can be formulated as a binary quadratic program [50]

$\min_{x_{i}^{2}=1,i=1,\ldots,n}\;\;\frac{1}{4}x^{\mathsf{T}}Lx,\\\ $ (48)

where $L\in\mathbb{S}^{n}$ is the Laplacian matrix of
$\mathcal{G}(\mathcal{V},\mathcal{E})$ defined as $L_{ij}=\sum_{j\neq
i}w_{ij}$ if $i=j$, and $L_{ij}=-w_{ij}$ otherwise. A well-known semidefinite
relaxation [50] for the Max-Cut problem 48 is

$\displaystyle\min_{X}$ $\displaystyle\frac{1}{4}\langle L,X\rangle$ (49)
$\displaystyle\mathrm{subject\leavevmode\nobreak\ to}$ $\displaystyle
X_{ii}=1,\ldots,n,$ $\displaystyle X\in\mathbb{S}^{n}_{+}.$

If the optimal solution of SDP relaxation 49 satisfies
$\mathrm{rank}(X^{\star})=1$, then the SDP relaxation 49 is exact and one can
recover a globally optimal solution to 48. However, the rank-one solution may
not exist. Instead, many max-cut instances admit low-rank optimal solutions:
$1<\mathrm{rank}(X^{\star})\ll n$, as observed in [17, Section 1.3.2], [29,
Section 5.1]. For these SDP instances, we expect that
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD exhibits faster linear convergence when
choosing a small value of $r_{\mathrm{c}}$, while
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP only has slower sublinear convergence
for the same choice of $r_{\mathrm{c}}$.

We run both $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP and
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD for a fixed number of 300 iterations
for two Max-Cut instances G1 and G25 from [57]. There are 800 nodes in graph
G1, and 2000 nodes in graph G25, and so the PSD constraints in 49 have a
dimension of 800 and 2000, respectively. Despite the large value of $n$, we
observed a low rank primal solution $\mathrm{rank}(X^{\star})=12$ for G1 and
$\mathrm{rank}(X^{\star})=19$ for G25. Note that the SDP 49 has a constant
trace property $\mathop{\bf tr}(X)=n$. Thus, we chose the penalty parameter
$\rho=2n+2$ for $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD. We can also estimate
the penalty parameter for $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP as666This
estimate is due to the structure of Max-Cut problems 49 that $\mathop{\bf
tr}\mathopen{}\mathclose{{}\left(\sum_{i=1}^{m}A_{i}y_{i}}\right)=b^{\mathsf{T}}y,\forall
y\in\mathbb{R}^{m}$, and the fact that
$\hat{y}=\mathopen{}\mathclose{{}\left(\frac{n}{4}\min\\{0,\lambda_{\min}(L)\\}}\right)\mathbf{1}$
is a dual feasible solution, where $\mathbf{1}\in\mathbb{R}^{m}$ is an all one
vector. Hence, we have $\mathop{\bf tr}(Z^{\star})=\frac{1}{4}\mathop{\bf
tr}(L)-\sum_{i=1}^{m}\mathop{\bf
tr}(A_{i})y^{\star}_{i}=\frac{1}{4}\mathop{\bf
tr}(L)-b^{\mathsf{T}}y^{\star}\leq\frac{1}{4}\mathop{\bf
tr}(L)-b^{\mathsf{T}}\hat{y},\forall(y^{\star},Z^{\star})\in\mathcal{D}^{\star}.$
$\rho=2\mathopen{}\mathclose{{}\left(\frac{1}{4}\mathop{\bf
tr}(L)-\frac{n}{4}\min\\{0,\lambda_{\min}(L)\\}}\right)+2$. We set
$r_{\mathrm{p}}=0$ since this parameter has little impact on convergence in
this case. The other parameters were chosen as those in Section 6.2.

The numerical results are illustrated in Figure 2 and Table 3. In all cases,
compared with $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP,
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD returns solutions with much higher
accuracy within the same number of iterations. In particular, since the rank
condition 45 is expected to hold, $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD shows
a faster linear convergence rate, while $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP
converges to an optimal solution slowly. Again, this is consistent with the
theoretical expectation in Section 5.2.

Table 3: Computational results of $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP and $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD for solving max-cut SDPs 49 under different values of $r_{\mathrm{p}}$ and $r_{\mathrm{c}}$. “Semi Feasi.” denotes the violation of PSD constraints $\max\\{\eta_{2},\eta_{4}\\}$, “Affine Feasi.” denotes the violation of affine constraints $\max\\{\eta_{1},\eta_{3}\\}$, “Dual Gap” denotes the duality gap $\eta_{5}$ in 47, and “Cost Opt.” denotes the cost-value gap $\eta_{\mathrm{c}}={(f(\omega_{t})-f^{\star})}/{f^{\star}}$. Max-cut instance | Algorithm | $(r_{\mathrm{p}},r_{\mathrm{c}})$ | Semi Feasi. | Affine Feasi. | Dual Gap | Cost Opt.
---|---|---|---|---|---|---
$\begin{aligned} &\text{G1}\\\ &\text{rank}(X^{\star})=12\\\ &n=800\\\ &m=800\end{aligned}$ |  $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP | $(0,12)$ | $-1.39\mathrm{e}{-2}$ | $7.28\mathrm{e}{-2}$ | $1.68\mathrm{e}{-2}$ | $6.85\mathrm{e}{-2}$
$(0,13)$ | $-1.39\mathrm{e}{-2}$ | $6.94\mathrm{e}{-2}$ | $1.77\mathrm{e}{-2}$ | $5.49\mathrm{e}{-2}$
$(0,14)$ | $-2.36\mathrm{e}{-2}$ | $6.63\mathrm{e}{-2}$ | $1.47\mathrm{e}{-2}$ | $5.97\mathrm{e}{-2}$
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD | $(0,12)$ | $-3.80\mathrm{e}{-6}$ | $7.16\mathrm{e}{-4}$ | $1.23\mathrm{e}{-5}$ | $6.86\mathrm{e}{-7}$
$(0,13)$ | $-2.52\mathrm{e}{-9}$ | $2.69\mathrm{e}{-5}$ | $4.58\mathrm{e}{-8}$ | $2.31\mathrm{e}{-9}$
$(0,14)$ | $-1.69\mathrm{e}{-9}$ | $4.42\mathrm{e}{-5}$ | $1.83\mathrm{e}{-7}$ | $1.97\mathrm{e}{-9}$
$\begin{aligned} &\text{G25}\\\ &\text{rank}(X^{\star})=19\\\ &n=2000\\\ &m=2000\end{aligned}$ |  $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP | $(0,18)$ | $-1.27\mathrm{e}{-1}$ | $4.73\mathrm{e}{-1}$ | $3.88\mathrm{e}{-2}$ | $3.07\mathrm{e}{-1}$
$(0,19)$ | $-9.24\mathrm{e}{-2}$ | $1.18$ | $2.33\mathrm{e}{-1}$ | $2.51\mathrm{e}{-1}$
$(0,20)$ | $-8.27\mathrm{e}{-2}$ | $6.62\mathrm{e}{-1}$ | $1.30\mathrm{e}{-1}$ | $2.24\mathrm{e}{-1}$
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD | $(0,18)$ | $-2.17\mathrm{e}{-7}$ | $6.34\mathrm{e}{-4}$ | $3.63\mathrm{e}{-6}$ | $9.37\mathrm{e}{-8}$
$(0,19)$ | $-1.72\mathrm{e}{-9}$ | $3.04\mathrm{e}{-6}$ | $1.77\mathrm{e}{-8}$ | $2.76\mathrm{e}{-9}$
$(0,20)$ | $-2.83\mathrm{e}{-9}$ | $2.76\mathrm{e}{-5}$ | $2.99\mathrm{e}{-8}$ | $3.01\mathrm{e}{-9}$

###  6.4 Quartic polynomial optimization on a sphere

In our last numerical experiments, we consider SOS relaxations for polynomial
optimization, which are likely to admit low-rank dual solutions. We expect
that $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP is better suited than
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD. Our numerical results further show
that $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP outperforms a set of baseline
solvers, including interior-point solvers SDPT3 [40], MOSEK [41], and first-
order solvers CDCS [42], SDPNAL+ [16].

Consider a constrained polynomial optimization problem over a sphere
$\min_{x\in\mathcal{S}^{n-1}}\;p_{0}(x),$ where
$p_{0}(x):\mathbb{R}^{n}\rightarrow\mathbb{R}$ is a polynomial and
$\mathcal{S}^{n-1}=\\{x\in\mathbb{R}^{n}\mid\|x\|^{2}=1\\}$ is the
$n$-dimensional unit sphere. This problem is in general NP-hard, but it can be
approximated well using the moment/SOS relaxation [35, 36]. In particular, the
$k$th-order SOS relaxation is

$\displaystyle\max_{\gamma,\sigma_{0},\psi_{1}}$ $\displaystyle\quad\gamma$
(50) $\displaystyle\mathrm{subject\leavevmode\nobreak\ to}$
$\displaystyle\quad p_{0}(x)-\gamma-\psi_{1}h(x)=\sigma_{0},$
$\displaystyle\quad\sigma_{0}\in\Sigma[x]_{n,2k},\;\psi_{1}\in\mathbb{R}[x]_{n,2(k-\lceil
h\rceil)},$

where $h(x)=\|x\|^{2}-1$, $\lceil h\rceil=\lceil\mathrm{deg}(h)/2\rceil$,
$\mathbb{R}[x]_{n,2(k-\lceil h\rceil)}$ denotes the real polynomial in $n$
variables and degree at most $2(k-\lceil h\rceil)$, and $\Sigma[x]_{n,2k}$
denote the cone of SOS polynomials in $\mathbb{R}[x]_{n,2k}$. It is well-known
that 50 can be equivalently reformulated into the standard primal SDP in the
form P with some extra free variables; see e.g., [35] for details. We observe
that these SDPs are likely to admit low-rank dual solutions (this observation
is consistent with the flat extension theory on the moment side [36, Theorem
3.7] when it is formulated as a dual SDP).

Motivated by the benchmark problems [58], we consider three instances of 50 in
our numerical experiments:

  1. 1.

Modified Broyden tridiagonal polynomial

$\displaystyle q_{1}(x)=$
$\displaystyle\mathopen{}\mathclose{{}\left(\mathopen{}\mathclose{{}\left(3-2x_{1}}\right)x_{1}-2x_{2}+1}\right)^{2}+$
$\displaystyle\sum_{i=2}^{d-1}\mathopen{}\mathclose{{}\left(\mathopen{}\mathclose{{}\left(3-2x_{i}}\right)x_{i}-x_{i-1}-2x_{i+1}+1}\right)^{2}+\mathopen{}\mathclose{{}\left(\mathopen{}\mathclose{{}\left(3-2x_{d}}\right)x_{d}-x_{d-1}+1}\right)^{2}+\mathopen{}\mathclose{{}\left(\sum_{i=1}^{d}x_{i}}\right)^{2}.$

  2. 2.

Modified Rosenbrock polynomial
$q_{2}(x)=1+\sum_{i=2}^{d}100\mathopen{}\mathclose{{}\left(x_{i}-x_{i-1}^{2}}\right)^{2}+\mathopen{}\mathclose{{}\left(1-x_{i}}\right)^{2}+\mathopen{}\mathclose{{}\left(\sum_{i=1}^{d}x_{i}}\right)^{2}.$

  3. 3.

Random quartic polynomial $q_{3}(x)=\langle c_{d,4},[x]_{d,4}\rangle,$ where
$[x]_{d,4}$ is the standard monomial bases with $d$ variables and degree at
most $4$, and $c_{d,4}$ is a randomly generated coefficient vector.

We used the package SOSTOOLS [59] to recast the SOS relaxation 50 with $k=2$
into a standard SDP in the form of P. We tested
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP for the above three polynomials with
different dimensions (the performance of
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMD was very poor in our experiments, and
we omitted it here). The dimension of the SDP relaxations ranges from $n=496$
to $861$ and $m=45,897$ to $134,889$. Recall in Section 3.1 that we need to
choose the penalty parameter $\rho>\mathcal{D}_{Z_{\star}}$ for
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP. For SDPs from 50, we can show that any
$\rho\geq(1+1)^{2}=4$ is a valid exact penalty parameter (see Lemma E.1)
thanks to the unit sphere constraint. We thus chose $\rho=10$ for all cases.
The parameters $(r_{\mathrm{p}},r_{\mathrm{c}})$ are chosen as $(0,3),(0,5),$
and $(0,7)$ for three different problems. In our experiments, we ran
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP until it reached tolerance
$1\mathrm{e}{-4}$.

Table 4: Computational results on the SOS relaxation 50 of quartic polynomial
optimization on a sphere: $\epsilon_{\mathrm{p}}$ denotes the primal
feasibility, $\epsilon_{\mathrm{d}}$ denotes the primal feasibility, and
$\epsilon_{\mathrm{g}}$ denotes the duality gap. Entries marked by oom
indicate “out of memory” runtime errors.

Dimension Case Metric SDPT3 [40] MOSEK [41] CDCS [42] SDPNAL+ [16]
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP $\begin{subarray}{c}d:\ 30\\\ \\\ n:\
496\ \\\ \\\ m:\ 45,879\ \end{subarray}$ $q_{1}$ $\epsilon_{\mathrm{p}}$ oom
$9.2\mathrm{e}{-11}$ $9.9\mathrm{e}{-5}$ $2.1\mathrm{e}{-11}$
$5.9\mathrm{e}{-14}$ $\epsilon_{\mathrm{d}}$ $2.0\mathrm{e}{-9}$
$4.1\mathrm{e}{-5}$ $9.5\mathrm{e}{-5}$ $9.2\mathrm{e}{-5}$
$\epsilon_{\mathrm{g}}$ $5.9\mathrm{e}{-12}$ $9.0\mathrm{e}{-6}$
$1.0\mathrm{e}{-2}$ $6.8\mathrm{e}{-6}$ cost $-25.074$ $-25.073$ $-42.850$
$-25.074$ time $887.1$ $100.9$ $1209$ 17.5 $q_{2}$ $\epsilon_{\mathrm{p}}$ oom
$8.2\mathrm{e}{-10}$ $1.5\mathrm{e}{-5}$ $4.7\mathrm{e}{-5}$
$5.4\mathrm{e}{-14}$ $\epsilon_{\mathrm{d}}$ $1.9\mathrm{e}{-8}$
$9.3\mathrm{e}{-5}$ $8.1\mathrm{e}{-5}$ $8.9\mathrm{e}{-5}$
$\epsilon_{\mathrm{g}}$ $2.9\mathrm{e}{-10}$ $3.9\mathrm{e}{-5}$
$6.7\mathrm{e}{-3}$ $5.3\mathrm{e}{-5}$ cost $-32.136$ $-32.137$ $-32.137$
$-32.135$ time $583$ $45.1$ $17.3$ 12.8 $q_{3}$ $\epsilon_{\mathrm{p}}$ oom
$9.1\mathrm{e}-{9}$ $1.0\mathrm{e}{-4}$ $3.2\mathrm{e}{-16}$
$1.9\mathrm{e}{-14}$ $\epsilon_{\mathrm{d}}$ $2.0\mathrm{e}{-7}$
$1.5\mathrm{e}{-6}$ $5.2\mathrm{e}{-4}$ $8.7\mathrm{e}{-5}$
$\epsilon_{\mathrm{g}}$ $4.9\mathrm{e}{-9}$ $2.4\mathrm{e}{-6}$
$3.5\mathrm{e}{-3}$ $2.2\mathrm{e}{-5}$ cost $2.023$ $2.023$ $2.021$ $2.023$
time $645.2$ $55.3$ $548.6$ 16.9 $\begin{subarray}{c}d:\ 35\\\ \\\ n:\ 666\\\
\\\ m:\ 81,584\end{subarray}$ $q_{1}$ $\epsilon_{\mathrm{p}}$ oom oom
$9.6\mathrm{e}{-5}$ $1.1\mathrm{e}{-15}$ $9.9\mathrm{e}{-14}$
$\epsilon_{\mathrm{d}}$ $1.9\mathrm{e}{-5}$ $6.8\mathrm{e}{-5}$
$9.9\mathrm{e}{-5}$ $\epsilon_{\mathrm{g}}$ $2.1\mathrm{e}{-5}$
$2.7\mathrm{e}{-4}$ $9.4\mathrm{e}{-6}$ cost $-30.053$ $-30.070$ $-30.050$
time $289.6$ $3449.4$ 41.4 $q_{2}$ $\epsilon_{\mathrm{p}}$ oom oom
$2.2\mathrm{e}{-5}$ $4.2\mathrm{e}{-5}$ $1.3\mathrm{e}{-13}$
$\epsilon_{\mathrm{d}}$ $9.4\mathrm{e}{-5}$ $9.5\mathrm{e}{-5}$
$8.8\mathrm{e}{-5}$ $\epsilon_{\mathrm{g}}$ $7.7\mathrm{e}{-5}$
$2.4\mathrm{e}{-3}$ $2.0\mathrm{e}{-5}$ cost $-37.100$ $-37.221$ $-37.091$
time $112.3$ $100.5$ 42.7 $q_{3}$ $\epsilon_{\mathrm{p}}$ oom oom
$1.0\mathrm{e}{-4}$ $9.3\mathrm{e}{-13}$ $1.9\mathrm{e}{-14}$
$\epsilon_{\mathrm{d}}$ $6.9\mathrm{e}{-5}$ $6.5\mathrm{e}{-7}$
$9.4\mathrm{e}{-5}$ $\epsilon_{\mathrm{g}}$ $3.5\mathrm{e}{-5}$
$6.2\mathrm{e}{-3}$ $2.5\mathrm{e}{-6}$ cost $2.121$ $2.039$ $2.121$ time
$170.1$ $1208.4$ 44.5 $\begin{subarray}{c}d:\ 40\\\ \\\ n:\ 861\\\ \\\ m:\
134,889\end{subarray}$ $q_{1}$ $\epsilon_{\mathrm{p}}$ oom oom
$1.0\mathrm{e}{-4}$ $3.6\mathrm{e}{-12}$ $7.4\mathrm{e}{-14}$
$\epsilon_{\mathrm{d}}$ $3.2\mathrm{e}{-6}$ $1.0\mathrm{e}{-4}$
$9.6\mathrm{e}{-5}$ $\epsilon_{\mathrm{g}}$ $1.4\mathrm{e}{-5}$
$4.9\mathrm{e}{-2}$ $5.2\mathrm{e}{-6}$ cost $-35.037$ $-47.792$ $-35.034$
time $738.9$ $1881.1$ 84.8 $q_{2}$ $\epsilon_{\mathrm{p}}$ oom oom
$1.7\mathrm{e}{-5}$ $3.3\mathrm{e}{-14}$ $1.3\mathrm{e}{-13}$
$\epsilon_{\mathrm{d}}$ $9.3\mathrm{e}{-5}$ $9.9\mathrm{e}{-5}$
$5.5\mathrm{e}{-5}$ $\epsilon_{\mathrm{g}}$ $5.8\mathrm{e}{-5}$
$2.8\mathrm{e}{-5}$ $1.6\mathrm{e}{-5}$ cost $-42.041$ $-41.997$ $-42.050$
time $254.9$ $187.6$ 82.1 $q_{3}$ $\epsilon_{\mathrm{p}}$ oom oom
$1.0\mathrm{e}{-4}$ $9.1\mathrm{e}{-16}$ $2.0\mathrm{e}{-14}$
$\epsilon_{\mathrm{d}}$ $1.6\mathrm{e}{-6}$ $7.5\mathrm{e}{-5}$
$8.2\mathrm{e}{-5}$ $\epsilon_{\mathrm{g}}$ $4.9\mathrm{e}{-5}$
$2.4\mathrm{e}{-3}$ $1.0\mathrm{e}{-4}$ cost $2.681$ $2.695$ $2.681$ time
$144.5$ $1339.6$ 93.1

To further demonstrate the performance of
$(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP, we compare it with SDPT3 [40], MOSEK
[41], CDCS [42], SDPNAL+ [16]. For SDPT3 and MOSEK, we used their default
parameters. For CDCS, we use the sos solver with a maximum 10,000 iterations
with tolerance $1\mathrm{e}{-4}$, and this sos option is customized for
solving SDPs from SOS relaxations (by exploiting a property called partial
orthogonality [60]). For SDPNAL+, we used $1\mathrm{e}{-4}$ as the tolerance,
turn off their stagnation detection, and run it using their default parameter
with a maximum of 20,000 iterations and a maximum of 10,000 seconds runtime.

The computational results are listed in Table 4. To be consistent with other
solvers, we report the cost value, time consumption, primal feasibility, dual
feasibility, and duality gap (see 47) of the final outcome, i.e.,

$\epsilon_{\mathrm{p}}=\eta_{1},\;\epsilon_{\mathrm{d}}=\eta_{3},\;\epsilon_{\mathrm{g}}=\eta_{5}.$

As we can see in Table 4, our algorithm $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP
solves all SDP instances to the desired accuracy within a reasonable time, and
it consistently outperforms the baseline solvers. For the interior-point
solvers, $\sf{SDPT3}$ runs out of memory in all cases on our computer.
$\sf{MOSEK}$ was able to return solutions of high accuracy for medium-size
problems ($d=30$) but took more time consumption. $\sf{MOSEK}$ also
encountered memory issues for larger instances ($d\geq 35$). For the first-
order solvers, CDCS solved all tested problems with medium accuracy, but the
runtime was worse than $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP (indeed, our
algorithm $(r_{\mathrm{p}},r_{\mathrm{c}})$-SBMP was one order of magnitude
faster than CDCS in some cases); the solver SDPNAL+ solved all SDPs to the
desired accuracy for the measures $\epsilon_{\mathrm{p}}$ and
$\epsilon_{\mathrm{d}}$, while the duality gap $\epsilon_{\mathrm{g}}$ remains
unsatisfactory. We note that the design of SDPNAL+ does not consider the
duality gap $\epsilon_{\mathrm{g}}$ as a stopping criterion. This partially
explains the poor performance of the duality gap $\epsilon_{\mathrm{g}}$ in
the final iterate.

##  7 Conclusion

In this paper, we have presented an overview and comparison of spectral bundle
methods for solving primal and dual SDPs. All the existing results focus on
solving dual SDPs. We have established a family of spectral bundle methods for
solving primal SDPs directly. The algorithm developments mirror the elegant
duality between primal and dual SDPs. We have presented the sublinear
convergence rates for this family of spectral bundle methods and shown that
the algorithm enjoys linear convergence with proper parameter choice and low-
rank dual solutions. The convergence behaviors and computational complexity of
spectral bundle methods for both primal and dual SDPs are in general similar,
but they have different features. It is clear that the existing spectral
bundle methods are well-suited for SDPs with low-rank primal solutions, and
our new spectral bundle method works well for SDPs with low-rank dual
solutions. These theoretical findings are supported by a range of large-scale
numerical experiments. We have further demonstrated that our new spectral
bundle method achieves state-of-the-art efficiency and scalability when
solving the SDP relaxations from polynomial optimization.

Potential future directions include incorporating other types of constraints
(such as nonnegative, second-order cone constraints, etc.), considering
second-order information [33] for the lower approximation, and analyzing the
algorithm performance when the subproblem 24 is solved inexactly [61].
Finally, we remark that our current prototype implementation shows promising
numerical performance, and it would also be very interesting to further
develop reliable and efficient open-source implementations of these spectral
bundle methods.

## References

  * [1] Lieven Vandenberghe and Stephen Boyd.  Semidefinite programming.  SIAM review, 38(1):49–95, 1996.
  * [2] Grigoriy Blekherman, Pablo A Parrilo, and Rekha R Thomas.  Semidefinite optimization and convex algebraic geometry.  SIAM, 2012.
  * [3] Yang Zheng, Giovanni Fantuzzi, and Antonis Papachristodoulou.  Chordal and factor-width decompositions for scalable semidefinite and polynomial optimization.  Annual Reviews in Control, 52:243–279, 2021.
  * [4] Lieven Vandenberghe, Martin S Andersen, et al.  Chordal graphs and semidefinite optimization.  Foundations and Trends® in Optimization, 1(4):241–433, 2015.
  * [5] Anirudha Majumdar, Georgina Hall, and Amir Ali Ahmadi.  Recent scalability improvements for semidefinite programming with applications in machine learning, control, and robotics.  Annual Review of Control, Robotics, and Autonomous Systems, 3:331–360, 2020.
  * [6] Stephen Boyd, Laurent El Ghaoui, Eric Feron, and Venkataramanan Balakrishnan.  Linear matrix inequalities in system and control theory.  SIAM, 1994.
  * [7] Renata Sotirov.  SDP relaxations for some combinatorial optimization problems.  In Handbook on Semidefinite, Conic and Polynomial Optimization, pages 795–819. Springer, 2012.
  * [8] Gert RG Lanckriet, Nello Cristianini, Peter Bartlett, Laurent El Ghaoui, and Michael I Jordan.  Learning the kernel matrix with semidefinite programming.  Journal of Machine learning research, 5(Jan):27–72, 2004.
  * [9] Lieven Vandenberghe and Stephen Boyd.  Applications of semidefinite programming.  Applied Numerical Mathematics, 29(3):283–299, 1999.
  * [10] Henry Wolkowicz, Romesh Saigal, and Lieven Vandenberghe.  Handbook of semidefinite programming: theory, algorithms, and applications, volume 27.  Springer Science & Business Media, 2012.
  * [11] Amir Ali Ahmadi, Georgina Hall, Antonis Papachristodoulou, James Saunderson, and Yang Zheng.  Improving efficiency and scalability of sum of squares optimization: Recent advances and limitations.  In 2017 IEEE 56th annual conference on decision and control (CDC), pages 453–462. IEEE, 2017.
  * [12] Zaiwen Wen, Donald Goldfarb, and Wotao Yin.  Alternating direction augmented lagrangian methods for semidefinite programming.  Mathematical Programming Computation, 2(3):203–230, 2010.
  * [13] Brendan O’Donoghue, Eric Chu, Neal Parikh, and Stephen Boyd.  Conic optimization via operator splitting and homogeneous self-dual embedding.  Journal of Optimization Theory and Applications, 169(3):1042–1068, June 2016.
  * [14] Michael Garstka, Mark Cannon, and Paul Goulart.  Cosmo: A conic operator splitting method for convex conic problems.  Journal of Optimization Theory and Applications, 190(3):779–810, 2021.
  * [15] Xin-Yuan Zhao, Defeng Sun, and Kim-Chuan Toh.  A Newton-CG augmented Lagrangian method for semidefinite programming.  SIAM Journal on Optimization, 20(4):1737–1765, 2010.
  * [16] Liuqin Yang, Defeng Sun, and Kim-Chuan Toh.  SDPNAL+: a majorized semismooth Newton-CG augmented Lagrangian method for semidefinite programming with nonnegative constraints.  Mathematical Programming Computation, 7(3):331–366, 2015.
  * [17] Alp Yurtsever, Joel A Tropp, Olivier Fercoq, Madeleine Udell, and Volkan Cevher.  Scalable semidefinite programming.  SIAM Journal on Mathematics of Data Science, 3(1):171–200, 2021\.
  * [18] Lijun Ding, Alp Yurtsever, Volkan Cevher, Joel A Tropp, and Madeleine Udell.  An optimal-storage approach to semidefinite programming using approximate complementarity.  SIAM Journal on Optimization, 31(4):2695–2725, 2021.
  * [19] Christoph Helmberg and Franz Rendl.  A spectral bundle method for semidefinite programming.  SIAM Journal on Optimization, 10(3):673–696, 2000.
  * [20] Christoph Helmberg and Krzysztof C Kiwiel.  A spectral bundle method with bounds.  Mathematical Programming, 93(2):173–194, 2002.
  * [21] Yang Zheng and Giovanni Fantuzzi.  Sum-of-squares chordal decomposition of polynomial matrix inequalities.  Mathematical Programming, 197(1):71–108, 2023.
  * [22] Mituhiro Fukuda, Masakazu Kojima, Kazuo Murota, and Kazuhide Nakata.  Exploiting sparsity in semidefinite programming via matrix completion i: General framework.  SIAM Journal on optimization, 11(3):647–674, 2001.
  * [23] Yang Zheng, Giovanni Fantuzzi, Antonis Papachristodoulou, Paul Goulart, and Andrew Wynn.  Chordal decomposition in operator-splitting methods for sparse semidefinite programs.  Mathematical Programming, 180(1):489–532, 2020.
  * [24] Yifan Sun, Martin S Andersen, and Lieven Vandenberghe.  Decomposition in conic optimization with partially separable structure.  SIAM Journal on Optimization, 24(2):873–897, 2014.