# Fast, Scalable, Warm-Start Semidefinite Programming
with Spectral Bundling and Sketching

Rico Angell Manning College of Information and Computer Sciences, University
of Massachusetts Amherst, Amherst, MA, USA. Correspondence to Rico Angell
(rangell@cs.umass.edu). Andrew McCallum11footnotemark: 1

###### Abstract

While semidefinite programming (SDP) has traditionally been limited to
moderate-sized problems, recent algorithms augmented with matrix sketching
techniques have enabled solving larger SDPs. However, these methods achieve
scalability at the cost of an increase in the number of necessary iterations,
resulting in slower convergence as the problem size grows. Furthermore, they
require iteration-dependent parameter schedules that prohibit effective
utilization of warm-start initializations important in practical applications
with incrementally-arriving data or mixed-integer programming. We present
SpecBM, a provably correct, fast and scalable algorithm for solving massive
SDPs that can leverage a warm-start initialization to further accelerate
convergence. Our proposed algorithm is a spectral bundle method for solving
general SDPs containing both equality and inequality constraints. Moveover,
when augmented with an optional matrix sketching technique, our algorithm
achieves the dramatically improved scalability of previous work while
sustaining convergence speed. We empirically demonstrate the effectiveness of
our method, both with and without warm-starting, across multiple applications
with large instances. For example, on a problem with 600 million decision
variables, SpecBM achieved a solution of standard accuracy in less than 7
minutes, where the previous state-of-the-art scalable SDP solver requires more
than 16 hours. Our method solves an SDP with more than $10^{13}$ decision
variables on a single machine with 16 cores and no more than 128GB RAM; the
previous state-of-the-art method had not achieved an accurate solution after
72 hours on the same instance. We make our implementation in pure JAX publicly
available111https://github.com/rangell/specbm.

##  1 Introduction

Semidefinite programming (SDP) is a convex optimization paradigm capable of
modeling or approximating many practical problems in combinatorial
optimization [3], neural network verification [71, 26], robotics [73], optimal
experiment design [92], VLSI [91], and systems and control theory [15].
Despite their widespread applicability, practitioners often dismiss the use of
SDPs under the assumption that optimization is intractable at real-world
scale. This assumption is grounded in the prohibitively high computational
complexity of standard SDP solvers (e.g. interior point methods [5, 90]).

The main challenge when solving SDPs with standard constrained optimization
approaches is the high cost of projection onto the feasible region. Projecting
a symmetric matrix onto the semidefinite cone requires a full
eigendecomposition, an operation that scales cubicly with the problem
dimension. The high computational cost of this single operation severely
limits the applicability of projected gradient and ADMM-based approaches [18,
69]. Interior point methods, the de facto approach for solving small-to-
moderate sized SDPs, require even more computational resources.

Projection-free and conditional gradient-based methods have been proposed to
avoid the computational burden of projecting iterates onto the semidefinite
cone [106, 104, 10]. These methods avoid projection onto the feasible region
with various techniques including subgradient descent in the dual space and a
conditional gradient method on an unconstrained augmented objective function.
The bulk of the computational burden for these methods is typically an extreme
eigenvector calculation. The scalablity of these methods is still limited,
despite the improved per-iteration complexity, since they require storing the
entire primal matrix which scales quadratically with the problem dimension.

The most well-known method for scaling semidefinte programming without storing
the entire primal matrix in memory is the Burer–Monteiro (BM) factorization
heuristic [21]. The core idea of the BM method is to explicitly control the
memory usage of the optimization algorithm by restricting the primal matrix to
a low-rank factorization. However, this method sacrifices the convexity of the
problem for scalability. In some cases, the BM method either requires having a
high rank factorization, leading to burdensome memory requirements, or risks
optimization getting stuck in local minima [97]. Yurtsever et al. [107] extend
the conditional gradient augmented Lagrangian approach [106] with a matrix
sketching technique. This approach avoids storing the entire primal matrix
while maintaining the convexity of the problem leading to provable convergence
for general SDPs. A significant drawback of this approach is the slow
empirical convergence particularly as the SDP grows in size. In addition, this
method uses iteration-dependent parameters prohibiting it from reliably
leveraging a warm-start initialization.

Spectral bundle methods, first proposed by Helmberg and Rendl [42], are an
appealing framework for solving SDPs due to their low per-iteration
computational complexity and fast empirical convergence. Recently, Ding and
Grimmer [28] proved non-asymptotic convergence rates for a spectral bundle
method designed to solve equality-constrained SDPs. These theoretical results
quantified the computational benefits of such methods, corroborating previous
empirical results. Moreover, spectral bundle methods are amenable to the same
matrix sketching techniques used by Yurtsever et al. [107] and have no
iteration-dependent parameters.

In this work, we show that a family of spectral bundle methods designed for a
broader class of SDPs is a practical approach for solving large problem
instances quickly. The method is powerful in that it allows the user to
control the tradeoff between per-iteration complexity and the empirical speed
at which the algorithm converges. In addition, the method can be augmented
with a matrix sketching technique that can dramatically improve the
scalability of the algorithm while maintaining its fast convergence.

We demonstrate the empirical efficacy of our proposed spectral bundle method
on three types of practical large SDPs seeing staggering performance
improvements in comparison with the previous state-of-the-art scalable
semidefinte programming algorithm for general SDPs. SpecBM is able to obtain a
good quality solution to an SDP with over $10^{13}$ decision variables while
the previous state-of-the-art fails to reach an accurate solution within 72
hours. In addition, we see a more than 1000x speedup of the previous state-of-
the-art scalable SDP solver on an instance with 600 million decision
variables. Finally, we find that warm-starting can lead to more than a 100x
speedup in the convergence of SpecBM as compared to cold-starting while the
previous state-of-the-art often is not able to reliably take advantage of a
warm-start initialization.

###  1.1 Contributions.

Our contributions can be summarized as follows:

  * •

We present a spectral bundle method (SpecBM) integrated with an optional
matrix sketching technique for solving SDPs with both equality and inequality
constraints. (Sections 3 & 4)

  * •

We prove correctness and show non-asymptotic convergence rates for SpecBM
under several standard assumptions. (Section 3.3)

  * •

We empirically show that SpecBM is fast, scalable, and can more reliably
leverage a warm-start initialization on three types of practical SDPs.
(Section 5)

  * •

We provide a standalone implementation of the algorithm in pure JAX [19, 36],
enabling efficient execution on various hardware (CPU, GPU, and TPU) and wide-
spread use. We make our implementation available at:
https://github.com/rangell/specbm.

##  2 Preliminaries

In this section, we present the necessary notation, assumptions, and
background material necessary to understanding the spectral bundle method.

###  2.1 Semidefinite Programming

We begin with the notation necessary to define the semidefinite programming
problem. Let $\SS^{n}$ be the set of all real symmetric $n\times n$ matrices
and $\SS^{n}_{+}:=\\{X\in\SS^{n}:X\succeq 0\\}$ be the positive semidefinite
cone. Let $\langle\cdot,\cdot\rangle$ denote the standard inner product for
vectors and the Frobenius inner product for matrices. For any matrix
$M\in\SS^{n}_{+}$ let $\lambda_{\max}(M)$ denote the maximum eigenvalue of $M$
and let $\mathrm{tr}(M)$ denote the trace of $M$. Let
$\mathcal{A}:\SS^{n}\to\mathbb{R}^{m}$ be a given linear operator and
$\mathcal{A}\mathclose{\vphantom{)}}^{*}:\mathbb{R}^{m}\to\SS^{n}$ be its
adjoint, i.e. $\langle\mathcal{A}X,y\rangle=\langle
X,\mathcal{A}\mathclose{\vphantom{)}}^{*}y\rangle$ for any $X\in\SS^{n}$ and
$y\in\mathbb{R}^{m}$, which take the following form,

$\mathcal{A}X=\begin{bmatrix}\langle A_{1},X\rangle\\\ \vdots\\\ \langle
A_{m},X\rangle\\\
\end{bmatrix}\quad\textrm{and}\quad\mathcal{A}\mathclose{\vphantom{)}}^{*}y=\sum_{i=1}^{m}y_{i}A_{i},$

for given matrices $A_{i}\in\SS^{n}$. For any given
$\mathcal{I}\subset\\{1,\ldots,m\\}$, let $\mathcal{A}_{\mathcal{I}}$ denote
the linear operator restricted to matrices $A_{i}$ for $i\in\mathcal{I}$ and
$y_{\mathcal{I}}$ to be the corresponding subset of rows for any
$y\in\mathbb{R}^{m}$. In general, we utilize $(\cdot)_{\mathcal{I}}$ to
extract the subset of components corresponding to $\mathcal{I}$ from an
arbitrary vector in $\mathbb{R}^{m}$. Finally, let
$\mathcal{I}^{\prime}:=\\{1,\cdots,m\\}\setminus\mathcal{I}$ be the compliment
of $\mathcal{I}$.

Given fixed $C\in\SS^{n}$, $\mathcal{A}$, $b\in\mathbb{R}^{m}$, and
$\mathcal{I}$, the primal (P) and dual (D) semidefinite programs can be
written as follows:

$\displaystyle\operatorname*{maximize}_{X\in\SS^{n}_{+}}\quad\langle
C,X\rangle$ (P)
$\displaystyle\operatorname*{subject\;to}\;\;\;\mathcal{A}_{\mathcal{I}^{{}^{\prime}}}X=b_{\mathcal{I}^{{}^{\prime}}}$
$\displaystyle\qquad\qquad\quad\;\;\mathcal{A}_{\mathcal{I}}X\leq
b_{\mathcal{I}}$
$\displaystyle\operatorname*{minimize}_{y\in\mathbb{R}^{m}}\quad\langle
b,y\rangle$ (D)
$\displaystyle\operatorname*{subject\;to}\;\;\;C-\mathcal{A}\mathclose{\vphantom{)}}^{*}y\preceq
0$ $\displaystyle\qquad\qquad\quad\;\;y_{\mathcal{I}}\geq 0$

Note that while all SDPs of this form can technically be written as equality
constrained SDPs (i.e. _standard-form_ SDPs), doing so in practice is ill-
advised since the problem dimension, $n$, grows by one for every inequality
constraint. In addition, we also represent the primal feasibility set as the
convex set $\mathcal{K}:=\\{z\in\mathbb{R}^{m}:z_{\mathcal{I}}\leq
b_{\mathcal{I}},z_{\mathcal{I}^{\prime}}=b_{\mathcal{I}^{\prime}}\\}$, and
measure the primal infeasibility $\mathrm{dist}(\mathcal{A}X,\mathcal{K})$. We
use $\mathrm{dist}(z,\mathcal{Z})$ to denote the Euclidean distance from a
point $z$ to a closed set $\mathcal{Z}$ and $\mathrm{proj}_{\mathcal{Z}}(z)$
to denote the Euclidean projection of $z$ onto $\mathcal{Z}$ (i.e. the closest
point in $\mathcal{Z}$ to $z$). Equivalently,

$\mathrm{dist}(z,\mathcal{Z})=\inf_{z^{\prime}\in\mathcal{Z}}\|z-z^{\prime}\|=\|z-\mathrm{proj}_{\mathcal{Z}}(z)\|.$

We denote the solution sets of (P) and (D) by $\mathcal{X}_{\star}$ and
$\mathcal{Y}_{\star}$, respectively. We make the standard assumption that (P)
and (D) satisfy _strong duality_ , namely that the solution sets
$\mathcal{X}_{\star}$ and $\mathcal{Y}_{\star}$ are nonempty, compact, and
every pair of solutions
$(X_{\star},y_{\star})\in\mathcal{X}_{\star}\times\mathcal{Y}_{\star}$ has
zero duality gap: $p_{\star}:=\langle C,X_{\star}\rangle=\langle
b,y_{\star}\rangle=:d_{\star}$. Strong duality holds whenever Slater’s
condition holds and $\mathcal{A}$ is a surjective linear operator.
Additionally, we say that the SDP satisfies _strict complementarity_ if there
exists $(X_{\star},y_{\star})$ such that for induced dual slack matrix
$Z_{\star}:=C-\mathcal{A}\mathclose{\vphantom{)}}^{*}y_{\star}$,

$\mathrm{rank}(X_{\star})+\mathrm{rank}(Z_{\star})=n.$

Strict complementarity is satisfied by generic SDPs [4] and many well-
structured SDPs [29]. We denote the maximum nuclear norm of the primal
solution set as
$\mathbb{N}(\mathcal{X}_{\star}):=\sup_{X_{\star}\in\mathcal{X}_{\star}}\|X_{\star}\|\mathclose{\vphantom{)}}_{*}$.

In many important applications, the SDPs take a highly structured and sparse
forms that admit low-rank solutions. The cost matrix $C$ is often sparse,
containing many fewer than $n^{2}$ non-zero entries. Additionally, the number
of primal constraints $m$ is also much less $n^{2}$, often being proportional
to $n$ or the number of non-zero entries in $C$. We would like to exploit
these features with a provably fast and scalable optimization algorithm that
can solve any SDP.

We also use $\langle\cdot,\cdot\rangle$ to denote the standard inner product
on vectors and $\|\cdot\|$ to denote the induced $\ell^{2}$-norm. We use
$\|\cdot\|_{\infty}$ to represent the $\ell^{\infty}$-norm. We use
$\|\cdot\|\mathclose{\vphantom{)}}_{\textrm{F}}$ to denote the Frobenius norm,
$\|\cdot\|\mathclose{\vphantom{)}}_{*}$ to denote the nuclear norm, and
$\|\cdot\|\mathclose{\vphantom{)}}_{\textrm{op}}$ to denote the operator norm.
We use $\mathcal{O}(\cdot)$ to represent standard big-O asymptotic notation.
We use $\textrm{nnz}(M)$ to denote the number of non-zero elements in a
(sparse) matrix $M$.

###  2.2 Proximal Bundle Method

Proximal bundle methods [55, 63, 33, 49] are a class of optimization
algorithms for solving unconstrained convex minimization problems of the form

$\operatorname*{minimize}_{y\in\mathbb{R}^{m}}f(y)$ (1)

where $f:\mathbb{R}^{m}\to(-\infty,+\infty]$ is a proper closed convex
function that attains its minimum value, $\inf f$, on some nonempty set. At
each iteration, the proximal bundle method proposes an update to the current
iterate $y_{t}$ by applying a proximal step to an approximation of the
objective function, $\hat{f}_{t}$:

$\tilde{y}_{t+1}\leftarrow\operatorname*{arg\,min}_{y\in\mathbb{R}^{m}}\hat{f}_{t}(y)+\frac{\rho}{2}\|y-y_{t}\|^{2}$
(2)

where $\rho>0$. The next iterate $y_{t+1}$ is set equal to $\tilde{y}_{t+1}$
only when the decrease in the objective value is at least a predetermined
fraction of the decrease predicted by the model $\hat{f}_{t}$, i.e.

$\beta(f(y_{t})-\hat{f}_{t}(\tilde{y}_{t+1}))\leq f(y_{t})-f(\tilde{y}_{t+1})$
(3)

for some fixed $\beta\in(0,1)$. The iterations where (3) is satisfied (and
thus, $y_{t+1}\leftarrow\tilde{y}_{t+1}$) are referred to as _descent steps_.
Otherwise, the algorithm takes a _null step_ and sets $y_{t+1}\leftarrow
y_{t}$. Regardless of whether (3) is satisfied or not, $\tilde{y}_{t+1}$ is
used to construct the next model $\hat{f}_{t+1}$.

##### Model Requirements.

The model $\hat{f}_{t}$ can take many forms, but is usually constructed using
subgradients of $f$ at past and current iterates. Let $\partial
f(y):=\\{g:f(y^{\prime})\geq f(y)+\langle g,y^{\prime}-y\rangle,\,\forall
y^{\prime}\in\mathbb{R}^{m}\\}$ denote the subdifferential of $f$ at $y$
(i.e. the set of subgradients of $f$ evaluated at a point $y$). Following
prior work, we require the next model $\hat{f}_{t+1}$ to satisfy some mild
assumptions:

  1. 1.

Minorant.

$\hat{f}_{t+1}(y)\leq f(y)\quad\textrm{ for all }y\in\mathbb{R}^{m}$ (4)

  2. 2.

Objective subgradient lower bound. For any $g_{t+1}\in\partial
f(\tilde{y}_{t+1})$,

$\hat{f}_{t+1}(y)\geq f(\tilde{y}_{t+1})+\langle
g_{t+1},y-\tilde{y}_{t+1}\rangle\quad\textrm{ for all }y\in\mathbb{R}^{m}$ (5)

  3. 3.

Model subgradient lower bound. The first order optimality conditions for (2)
gives the subgradient
$s_{t+1}:=\rho(y_{t}-\tilde{y}_{t+1})\in\partial\hat{f}_{t}(\tilde{y}_{t+1})$.
After a null step $t$,

$\hat{f}_{t+1}(y)\geq\hat{f}_{t}(\tilde{y}_{t+1})+\langle
s_{t+1},y-\tilde{y}_{t+1}\rangle\quad\textrm{ for all }y\in\mathbb{R}^{m}$ (6)

The first two conditions are almost inherent, serving to guarantee that a new
model integrates first-order information from the objective at
$\tilde{y}_{t+1}$. The third condition is workable and, in essence, demands
the new model to preserve a degree of the approximation accuracy exhibited by
the preceding model.

##### Convergence Guarantees.

Recently, Díaz and Grimmer [27] proved non-asymptotic convergence rates for
the proximal bundle method under a wide variety of conditions. The following
theorem summarizes the results relevant to this work.

######  Theorem 2.1.

Let $\beta>0$ and $\rho>0$ be constant, $f:\mathbb{R}^{m}\to(-\infty,+\infty]$
be a proper closed convex function with nonempty set of minimizers
$\mathcal{Y}_{\star}$, and the sequence of models produced by the proximal
bundle method
$\left\\{\hat{f}_{t+1}:\mathbb{R}^{m}\to(-\infty,+\infty]\right\\}$ satisfy
the conditions (4), (5), and (6). If the iterates $y_{t}$ and subgradients
$g_{t+1}$ are bounded, then the iterates $y_{t}$ have
$f(y_{t})-f(y_{\star})\leq\varepsilon$ for all
$t\geq\mathcal{O}(1/\varepsilon^{3})$. Additionally, if
$f(y)-f(y_{\star})\geq\mathcal{O}(\mathrm{dist}^{2}(y,\mathcal{Y}_{\star}))$
for all $y$, then the bound improves to $f(y_{t})-f(y_{\star})\leq\varepsilon$
for all $t\geq\mathcal{O}(1/\varepsilon)$.

We will use this result to prove non-asymptotic convergence rates for the
spectral bundle method presented in the next section.

##  3 Spectral Bundle Method

In this section, we will present our proposed algorithm for applying the
proximal bundle method to the SDP defined in (P) and (D), called the spectral
bundle method (SpecBM). Our proposed spectral bundle method is a powerful
algorithm for solving SDPs with both equality and inequality constraints and
can also be integrated with a matrix sketching technique for dramatically
improved scalability as demonstrated in Section 4. To apply the proximal
bundle method, we first need to convert the constrained convex optimization
problem into an equivalent unconstrained problem. Begin by considering the
following Lagrangian formulation of the model problem

$p_{\star}=\max_{X\in\SS^{n}_{+}}\;\min_{y\in\mathbb{R}^{m}:y_{\mathcal{I}}\geq
0}\;\mathcal{L}(X,y)=\min_{y\in\mathbb{R}^{m}:y_{\mathcal{I}}\geq
0}\;\max_{X\in\SS^{n}_{+}}\;\mathcal{L}(X,y)=d_{\star},$ (7)

where the Lagrangian is defined as

$\mathcal{L}(X,y):=\langle
C-\mathcal{A}\mathclose{\vphantom{)}}^{*}y,X\rangle+\langle b,y\rangle.$ (8)

For any $\alpha>\sup_{X_{\star}\in\mathcal{X}_{\star}}\mathrm{tr}(X_{\star})$,
let $\mathcal{X}:=\\{X\in\SS^{n}_{+}:\mathrm{tr}(X)\leq\alpha\\}$ be the trace
constrained subset of the primal domain. Additionally, define
$\mathcal{Y}:=\\{y\in\mathbb{R}^{m}:y_{\mathcal{I}}\geq 0\\}$. Since
$\mathcal{X}\subset\mathcal{X}_{\star}$, we know

$\max_{X\in\SS^{n}_{+}}\;\min_{y\in\mathcal{Y}}\;\mathcal{L}(X,y)=\max_{X\in\mathcal{X}}\;\min_{y\in\mathcal{Y}}\;\mathcal{L}(X,y)=\min_{y\in\mathcal{Y}}\;\max_{X\in\mathcal{X}}\;\mathcal{L}(X,y)$
(9)

Furthermore, it is easy to show the following fact [70],

$\max_{X\succeq 0\,,\,\mathrm{tr}(X)\leq 1}\langle
Z,X\rangle=\max\\{\lambda_{\textrm{max}}(Z),0\\}.$ (10)

Incorporating (10) with (9) allows us to write an equivalent formulation of
the original dual problem as follows

$\operatorname*{minimize}_{y\in\mathbb{R}^{m}}f(y):=\alpha\max\\{\lambda_{\textrm{max}}(C-\mathcal{A}\mathclose{\vphantom{)}}^{*}y),0\\}+\langle
b,y\rangle+\iota_{\mathcal{Y}}(y),$ (pen-D)

where $\iota_{\mathcal{Y}}(\cdot)$ is the indicator function defined such that
$\iota_{\mathcal{Y}}(y)=0$ if $y\in\mathcal{Y}$ and
$\iota_{\mathcal{Y}}(y)=+\infty$ otherwise. The problem (pen-D) is an
unconstrained convex optimization problem as desired in order to apply the
proximal bundle method. In the following subsections, we will detail how the
model is constructed and how to compute the candidate iterates.

###  3.1 Spectral Bundle Model

The form of (pen-D) is not quite conducive to defining the model and solving
for the candidate iterate $\tilde{y}_{t+1}$. To address this, we will consider
an equivalent formulation to (pen-D) that is more amenable to this effort. We
begin by considering the dual form of the indicator function
$\iota_{\mathcal{Y}}(\cdot)$

$\iota_{\mathcal{Y}}(y)=\sup\\{-\langle\nu,y\rangle:\nu_{\mathcal{I}}\leq
0,\nu_{\mathcal{I}^{\prime}}=0\\}.$ (11)

Define $\mathrm{N}:=\\{\nu\in\mathbb{R}^{m}:\nu_{\mathcal{I}}\leq
0,\nu_{\mathcal{I}^{{}^{\prime}}}=0\\}$. Substituting (11) into (pen-D) and
utilizing (10) gives

$f(y)=\sup_{(X,\nu)\in\mathcal{X}\times\mathrm{N}}\langle
C-\mathcal{A}\mathclose{\vphantom{)}}^{*}y,X\rangle+\langle b-\nu,y\rangle.$
(12)

##### Defining the model.

This equivalent formulation is clearly just as difficult to optimize as
(pen-D), but is helpful in defining the model we will utilize in the spectral
bundle method. The main idea is to consider (12) over a low-dimensional
subspace of $\mathcal{X}$ such that (4), (5), and (6) are satisfied. The
subspace we will consider is parameterized by matrices
$\bar{X}_{t}\in\SS^{n}_{+}$ and $V\in\mathbb{R}^{n\times k}$ and is defined as
follows

$\widehat{\mathcal{X}}_{t}:=\left\\{\eta\bar{X}_{t}+V_{t}SV_{t}^{\top}:\eta\,\mathrm{tr}(\bar{X}_{t})+\mathrm{tr}(S)\leq\alpha,\eta\geq
0,S\in\SS^{k}_{+}\right\\}.$ (13)

The matrix $V_{t}$ has $k$ orthonormal column vectors, where $k$ is a small
user defined parameter. When $k=1$, $V_{t}$ is just a maximum eigenvector of
$C-\mathcal{A}\mathclose{\vphantom{)}}^{*}\tilde{y}_{t}$, which also means
$b-\alpha\mathcal{A}V_{t}V_{t}^{\top}\in\partial f(\tilde{y}_{t})$. When
$k>1$, the columns of $V_{t}$ include a maximum eigenvector as well as other
eigenvectors based on the candidate dual iterate $\tilde{y}_{t}$ and
(potentially) the accumulation of past spectral information. The columns of
$V_{t}$ are further partitioned into $k_{c}\geq 1$ current eigenvectors and
$k_{p}\geq 0$ orthonormal vectors representing past spectral information. For
scalability reasons, we expect $k=k_{c}+k_{p}$ to be relatively small. The
matrix $\bar{X}_{t}$ is a carefully selected weighted sum of past spectral
bounds such that $\mathrm{tr}(\bar{X}_{t})\leq\alpha$ and enables the last
model condition (6) to be satisfied. Hence, the model is

$\hat{f}_{t}(y):=\sup_{(X,\nu)\in\widehat{\mathcal{X}}_{t}\times\mathrm{N}}\langle
C-\mathcal{A}\mathclose{\vphantom{)}}^{*}y,X\rangle+\langle b-\nu,y\rangle.$
(14)

We show that this model satisfies the conditions (4), (5), and (6) in Appendix
B.6.

##### Updating the model.

In every step, independent of whether a descent step or null step is taken,
the model needs to be updated. At every step, we solve the following minimax
problem

$\min_{y\in\mathbb{R}^{m}}\sup_{(X,\nu)\in\widehat{\mathcal{X}}_{t}\times\mathrm{N}}\langle
C-\mathcal{A}\mathclose{\vphantom{)}}^{*}y,X\rangle+\langle
b-\nu,y\rangle+\frac{\rho}{2}\|y-y_{t}\|^{2},$ (15)

where $(\tilde{y}_{t+1},X_{t+1},\nu_{t+1})$ is the minimax solution. The
structure of $\widehat{\mathcal{X}}_{t}$ allows us to rewrite $X_{t+1}$ as

$X_{t+1}=\eta_{t+1}\bar{X}_{t}+V_{t}S_{t+1}V_{t}^{\top}.$ (16)

To compute $\bar{X}_{t+1}$ and $V_{t+1}$ we first compute an
eigendecomposition of $S_{t+1}$ to separate current and past spectral
information as follows

$S_{t+1}=Q_{\overline{p}}\Lambda_{\overline{p}}Q_{\overline{p}}^{\top}+Q_{\underline{c}}\Lambda_{\underline{c}}Q_{\underline{c}}^{\top},$
(17)

where $\Lambda_{\overline{p}}$ is a diagonal matrix containing the largest
$k_{p}$ eigenvalues and $Q_{\overline{p}}\in\mathbb{R}^{k\times k_{p}}$
contains the corresponding eigenvectors. The diagonal matrix
$\Lambda_{\underline{c}}$ contains the remaining $k_{c}$ eigenvalues and
$Q_{\underline{c}}\in\mathbb{R}^{k\times k_{c}}$ contains the corresponding
eigenvectors. It is important to note that we can compute this
eigendecomposition efficiently Given this decomposition, we set the next
model’s $\bar{X}_{t+1}$ as follows

$\bar{X}_{t+1}\leftarrow\eta_{t+1}\bar{X}_{t}+V_{t}Q_{\underline{c}}\Lambda_{\underline{c}}Q_{\underline{c}}^{\top}V_{t}^{\top}.$
(18)

In the case that $k_{p}=0$, observe that $\bar{X}_{t}=X_{t}$ for all $t$. To
update $V_{t+1}$, we first compute the top $k_{c}$ orthonormal eigenvectors
$v_{1},\ldots,v_{k_{c}}$ of
$C-\mathcal{A}\mathclose{\vphantom{)}}^{*}\tilde{y}_{t+1}$. Then, we set
$V_{t+1}$ to $k$ orthonormal vectors that span the columns of
$[V_{t}Q_{\overline{p}};v_{1},\ldots,v_{k_{c}}]$ as follows

$V_{t+1}\leftarrow\texttt{orthonormalize}\Big{(}[V_{t}Q_{\overline{p}};v_{1},\ldots,v_{k_{c}}]\Big{)},$
(19)

where we can use QR decompostion to orthonormalize the vectors. If $k_{p}=0$,
this process is unnecessary and we just set
$V_{t+1}\leftarrow[v_{1},\ldots,v_{k_{c}}]$.

###  3.2 Computing the Candidate Iterate

We will now discuss how to compute the candidate iterate $\tilde{y}_{t+1}$ by
applying a proximal step to the current approximation of the objective
function (i.e. solve (2)). Notice that the candidate iterate can be written as
follows

$\tilde{y}_{t+1}=\operatorname*{arg\,min}_{y\in\mathbb{R}^{m}}\sup_{(X,\nu)\in\widehat{\mathcal{X}}_{t}\times\mathrm{N}}F(X,\nu,y):=\langle
C-\mathcal{A}\mathclose{\vphantom{)}}^{*}y,X\rangle+\langle
b-\nu,y\rangle+\frac{\rho}{2}\|y-y_{t}\|^{2}.$ (20)

To solve for $\tilde{y}_{t+1}$, start by fixing
$(X,\nu)\in\widehat{\mathcal{X}}_{t}\times\mathrm{N}$ arbitrarily. Then,
completing the square gives

$\displaystyle\operatorname*{arg\,min}_{y\in\mathbb{R}^{m}}F(X,\nu,y)$
$\displaystyle=\operatorname*{arg\,min}_{y\in\mathbb{R}^{m}}\langle
b-\nu-\mathcal{A}X,y\rangle+\frac{\rho}{2}\|y-y_{t}\|^{2}$ (21)
$\displaystyle=\operatorname*{arg\,min}_{y\in\mathbb{R}^{m}}\frac{\rho}{2}\left\|y-y_{t}+\frac{1}{\rho}(b-\nu-\mathcal{A}X)\right\|^{2}$
$\displaystyle=y_{t}-\frac{1}{\rho}(b-\nu-\mathcal{A}X).$

This implies that the candidate iterate can be computed as follows

$\tilde{y}_{t+1}=\operatorname*{arg\,min}_{y\in\mathbb{R}^{m}}\hat{f}_{t}(y)+\|y-y_{t}\|^{2}=y_{t}-\frac{1}{\rho}(b-\nu_{t+1}-\mathcal{A}X_{t+1}),$
(22)

where $(X_{t+1},\nu_{t+1})\in\widehat{\mathcal{X}}_{t}\times\mathrm{N}$ is the
solution to the following optimization problem

$(X_{t+1},\nu_{t+1})\in\operatorname*{arg\,max}_{(X,\nu)\in\widehat{\mathcal{X}}_{t}\times\mathrm{N}}\psi_{t}(X,\nu):=\langle
C,X\rangle+\langle
b-\nu-\mathcal{A}X,y_{t}\rangle-\frac{1}{2\rho}\|b-\nu-\mathcal{A}X\|^{2}$
(23)

To solve for
$(X_{t+1},\nu_{t+1})\in\widehat{\mathcal{X}}_{t}\times\mathrm{N}$, we propose
using the alternating maximization algorithm (also referred to as block
coordinate ascent). After initializing $\tilde{\nu}=0$, then the alternating
maximization algorithm has the following update steps.

$\displaystyle\tilde{X}$
$\displaystyle\leftarrow\operatorname*{arg\,max}_{X\in\widehat{\mathcal{X}}_{t}}\;\psi_{t}(X,\tilde{\nu})$
(24) $\displaystyle\tilde{\nu}$
$\displaystyle\leftarrow\operatorname*{arg\,max}_{\nu\in\mathrm{N}}\;\psi_{t}(\tilde{X},\nu)=\mathrm{proj}_{\mathrm{N}}(b-\mathcal{A}\tilde{X}-\rho
y_{t})$ (25)

The steps (24) and (25) are repeated until convergence. Clearly $\psi_{t}$ is
smooth, and in this case the alternating minimization algorithm is known to
converge at a $\mathcal{O}(1/\varepsilon)$ rate [14]. If $\psi_{t}$ happens to
also be strongly concave, the rate improves to
$\mathcal{O}(\log(1/\varepsilon))$ [61]. We solve (24) using a primal-interior
point method derived in Appendix A.2222Ding and Grimmer [28] claim that (24)
could be solved using (accelerated) projected gradient descent and describe
the method necessary for projection and proper scaling, but they use Mosek
[9], an off-the-shelf solver, in their experiments. We find that while
theoretically possible, projected gradient descent does not work well in
practice since choosing the correct step size to obtain a high quality
solution to 24 is difficult and varies between time steps $t$. We instead use
(and advocate for) the primal-dual interior point method derived in Appendix
A.2. This interior point method has time complexity $\mathcal{O}(k^{6})$ and
requires $\mathcal{O}(k^{4})$ memory.

Algorithm 1 Spectral Bundle Method (SpecBM)

1: Input: Problem data $(C,\mathcal{A},b,\mathcal{I})$, parameters
$(\alpha,\rho,\beta,k_{c},k_{p})$, and initialization of $X_{0}=\bar{X}_{0}$,
$y_{0}$, and orthonormal $V_{0}$

2: Output: $X_{T},y_{T}$

3: for $t=0,1,\ldots,T$ do

4:
$(X_{t+1},\nu_{t+1})\leftarrow{\operatorname*{arg\,max}\limits_{(X,\nu)\in\widehat{\mathcal{X}}_{t}\times\mathrm{N}}}\psi_{t}(X,\nu)$
// where $X_{t+1}=\eta_{t+1}\bar{X}_{t}+V_{t}S_{t+1}V_{t}^{\top}$

5: $\tilde{y}_{t+1}\leftarrow
y_{t}-\frac{1}{\rho}(b-\nu_{t+1}-\mathcal{A}X_{t+1})$

6: if $\beta(f(y_{t})-\hat{f}_{t}(\tilde{y}_{t+1}))\leq
f(y_{t})-f(\tilde{y}_{t+1})$ then

7: $y_{t+1}\leftarrow\tilde{y}_{t+1}$ // descent step

8: else

9: $y_{t+1}\leftarrow y_{t}$ // null step

10: Compute the eigendecomposition
$S_{t+1}=Q_{\overline{p}}\Lambda_{\overline{p}}Q_{\overline{p}}^{\top}+Q_{\underline{c}}\Lambda_{\underline{c}}Q_{\underline{c}}^{\top}$
// see (17)

11:
$\bar{X}_{t+1}\leftarrow\eta_{t+1}\bar{X}_{t}+V_{t}Q_{\underline{c}}\Lambda_{\underline{c}}Q_{\underline{c}}^{\top}V_{t}^{\top}$

12: Compute $k_{c}$ maximum eigenvectors $v_{1},\ldots,v_{k_{c}}$ of
$C-\mathcal{A}^{*}\tilde{y}_{t+1}$

13:
$V_{t+1}\leftarrow\texttt{orthonormalize}\Big{(}[V_{t}Q_{\overline{p}};v_{1},\ldots,v_{k_{c}}]\Big{)}$

If we substitute
$\nu_{t+1}=\mathrm{proj}_{\mathrm{N}}(b-\mathcal{A}X_{t+1}-\rho y_{t})$ into
(22), it can be seen that the candidate $\tilde{y}_{t+1}$ is always feasible
and is also complementary to the dual slack variable $\nu_{t+1}$,

$\big{(}\tilde{y}_{t+1}\big{)}_{\mathcal{I}}=\max\left\\{\big{(}y_{t}\big{)}_{\mathcal{I}}-\frac{1}{\rho}\Big{(}b_{\mathcal{I}}-\mathcal{A}_{\mathcal{I}}X_{t+1}\Big{)},0\right\\}\geq
0\quad\textrm{and}\quad\langle\tilde{y}_{t+1},\nu_{t+1}\rangle=0.$ (26)

This view allows us to see that we can equivalently express the candidate
iterate as
$\tilde{y}_{t+1}=\mathrm{proj}_{\mathcal{Y}}\left(y_{t}-\frac{1}{\rho}(b-\mathcal{A}X_{t+1})\right)$.

The complete Spectral Bundle Method (SpecBM) is detailed in Algorithm 1. It is
important to note that the main cost of computing $f(y_{t})$ and
$f(\tilde{y}_{t+1})$ is a maximum eigenvalue computation of
$C-\mathcal{A}\mathclose{\vphantom{)}}^{*}y_{t}$ and
$C-\mathcal{A}\mathclose{\vphantom{)}}^{*}\tilde{y}_{t+1}$, respectively.
Lastly, we compute $\hat{f}_{t}(\tilde{y}_{t+1})$ using a primal-dual interior
point method derived in Appendix A.1. This interior point method has time
complexity $\mathcal{O}(k^{6})$ and requires $O(k^{4})$ memory.

###  3.3 Convergence Analysis

In this section, we present the non-asymptotic convergence results for SpecBM.
Under any setting of the algorithms parameters, Theorem 3.1 guarantees
sublinear convergence for both primal and dual problems, showing that the
iterates $X_{t}$ and $y_{t}$ converge in terms of objective gap and
feasilibity. The results are summarized in the following theorem.

######  Theorem 3.1.

Suppose strong duality holds. For any fixed $\rho>0,\beta\in(0,1),k_{c}\geq
1,k_{p}\geq 0$, and $\alpha\geq 2\mathbb{N}(\mathcal{X}_{\star})$, SpecBM
produces iterates $X_{t}\succeq 0$ and $y_{t}\in\mathcal{Y}$ such that for any
$\varepsilon\in(0,1]$,

approximate penalized dual optimality: $\displaystyle
f(y_{t})-f(y_{\star})\leq\varepsilon,$ approximate primal feasiblity:
$\displaystyle\mathrm{dist}(\mathcal{A}X_{t},\mathcal{K})\leq\sqrt{\varepsilon},$
approximate dual feasiblity:
$\displaystyle\lambda_{\max}(C-\mathcal{A}\mathclose{\vphantom{)}}^{*}y_{t})\leq\varepsilon,$
approximate primal-dual optimality: $\displaystyle|\langle
b,y_{t}\rangle-\langle C,X_{t}\rangle|\leq\sqrt{\varepsilon},$

by some iteration $\mathcal{O}(1/\varepsilon^{3})$. And, if strict
complementarity holds, then these conditions are achieved by some iteration
$\mathcal{O}(1/\varepsilon)$. Additionally, if Slater’s condition holds and
$y_{\star}$ is unique, then the approximate primal feasibility and approximate
primal-dual optimality respectively improve to

$dist(\mathcal{A}X_{t},\mathcal{K})\leq\varepsilon\quad\textrm{and}\quad|\langle
b,y_{t}\rangle-\langle C,X_{t}\rangle|\leq\varepsilon.$

The proof of this theorem is given in Appendix B.6 and the general strategy is
to show that the assumptions made in Theorem 2.1 hold (including the model
conditions (4), (5), and (6)) and then use the following lemmas to obtain the
stated result. Note the following fact necessary to use Theorem 2.1.

######  Fact 3.2.

The set $\\{y\in\mathcal{Y}:f(y_{\star})\leq f(y)\leq f(y_{0})\\}$ is compact
and contains all iterates $y_{t}$. Moreover, the iterates $y_{t}$ and
subgradients $g_{t+1}$ are bounded, i.e. $\sup_{\,t\geq
0}\\{\mathrm{dist}(y_{t},\mathcal{Y}_{*})\\}<\infty$ and $\sup_{\,t\geq
0}\left\\{\|g_{t+1}\|:g_{t+1}\in\partial f(\tilde{y}_{t+1})\right\\}<\infty$.

The following lemma guarantees quadratic growth whenever strict
complementarity holds.

######  Lemma 3.3.

Let $y\in\mathcal{Y}$ be arbitrary. Then,

$\mathrm{dist}^{2^{d}}(y,\mathcal{Y}_{\star})\leq\mathcal{O}(f(y)-f(y_{\star})),$
(27)

where $d\in\\{0,1,\ldots,m\\}$ is the singularity degree _[82, 31, 80]_. If
strict complementarity holds, then $d\leq 1$, and if Slater’s condition holds,
then $d=0$.

The following three lemmas guarantee primal feasiblity, dual feasibility, and
primal-dual objective optimality given convergence of the penalized dual
objective, i.e. $f(y_{t})-f(y_{\star})\leq\varepsilon$.

######  Lemma 3.4.

At every descent step $t$,

$\mathrm{dist}(\mathcal{A}X_{t+1},\mathcal{K})\leq\mathcal{O}\\!\left(\sqrt{f(y_{t})-f(y_{\star})}\right).$
(28)

Additionally, if Slater’s condition holds and $y_{\star}$ is unique, then

$\mathrm{dist}(\mathcal{A}X_{t+1},\mathcal{K})\leq\mathcal{O}(f(y_{t})-f(y_{\star})).$
(29)

######  Lemma 3.5.

Suppose strong duality holds and $\alpha\geq
2\mathbb{N}(\mathcal{X}_{\star})$. Then, at every descent step $t$,

$\lambda_{\max}(C-\mathcal{A}\mathclose{\vphantom{)}}^{*}y_{t})\leq\mathcal{O}\\!\left(f(y_{t})-f(y_{\star})\right).$
(30)

######  Lemma 3.6.

At every descent step $t$,

$\left|\langle b,y_{t+1}\rangle-\langle
C,X_{t+1}\rangle\right|\leq\mathcal{O}\\!\left(f(y_{t})-f(y_{\star})\right)+\mathcal{O}\\!\left(\sqrt{f(y_{t})-f(y_{\star})}\right).$
(31)

Additionally, if Slater’s condition holds and $y_{\star}$ is unique, then

$\left|\langle b,y_{t+1}\rangle-\langle
C,X_{t+1}\rangle\right|\leq\mathcal{O}(f(y_{t})-f(y_{\star})).$ (32)

The proofs of these results are given in Appendix B.3-B.5, respectively. It is
important to note that Ding and Grimmer [28] show for a similar spectral
bundle method designed for equality-constrained SDPs that under stronger
conditions, namely $k_{c}\geq\mathrm{rank}(X_{\star})$ (among others), that
$f(y_{t})-f(y_{\star})$ converges at a $\mathcal{O}(\log(1/\varepsilon))$
rate. It should be possible to extend that result to our algorithm, but we
leave this for future work.

##  4 Scaling with Matrix Sketching

The memory required to store the primal iterates $\bar{X}_{t}$ (similarly
$X_{t}$) is prohibitive to scaling SDPs to large problem instances. We will
utilize a _Nyström sketch_ [87, 37, 39, 56] to track a compressed low-rank
projection of the primal iterate as it evolves which at any iteration can be
used to compute a provably accurate low-rank approximation of $\bar{X}_{t}$.

First, notice that the operations carried out by SpecBM do not require
explicitly storing $\bar{X}_{t}$. In fact, if we are not interested in the
primal iterates and we only want to solve the dual problem (e.g. in the case
where we want to test the feasbility of the primal problem), we only need to
store $\langle C,\bar{X}_{t}\rangle$, $\mathrm{tr}(\bar{X}_{t})$, and
$\mathcal{A}\bar{X}_{t}$. This means that our representation of $\bar{X}_{t}$
is independent from the operations taken by SpecBM. In addition, these values
can be efficiently maintained given low-rank updates to $\bar{X}_{t}$ due to
the linearity of the operations,

$\displaystyle\langle C,\bar{X}_{t+1}$ $\displaystyle\rangle\leftarrow\langle
C,\eta_{t+1}\bar{X}_{t}+V_{t}Q_{\underline{c}}\Lambda_{\underline{c}}Q_{\underline{c}}^{\top}V_{t}^{\top}\rangle=\eta_{t+1}\langle
C,\bar{X}_{t}\rangle+\mathrm{tr}\Big{(}V_{t}^{\top}Q_{\underline{c}}^{\top}\big{(}CV_{t}Q_{\underline{c}}\Lambda_{\underline{c}}\big{)}\Big{)},$
(33) $\displaystyle\mathrm{tr}(\bar{X}_{t+1})$
$\displaystyle\leftarrow\mathrm{tr}\Big{(}\eta_{t+1}\bar{X}_{t}+V_{t}Q_{\underline{c}}\Lambda_{\underline{c}}Q_{\underline{c}}^{\top}V_{t}^{\top}\Big{)}=\eta_{t+1}\mathrm{tr}\big{(}\bar{X}_{t}\big{)}+\mathrm{tr}\\!\left(\Lambda_{\underline{c}}\right),$
(34) $\displaystyle\mathcal{A}\bar{X}_{t+1}$
$\displaystyle\leftarrow\mathcal{A}\Big{(}\eta_{t+1}\bar{X}_{t}+V_{t}Q_{\underline{c}}\Lambda_{\underline{c}}Q_{\underline{c}}^{\top}V_{t}^{\top}\Big{)}=\eta_{t+1}\mathcal{A}\bar{X}_{t}+\mathcal{A}\big{(}V_{t}Q_{\underline{c}}\Lambda_{\underline{c}}Q_{\underline{c}}^{\top}V_{t}^{\top}\big{)},$
(35)

where
$\mathcal{A}\big{(}V_{t}Q_{\underline{c}}\Lambda_{\underline{c}}Q_{\underline{c}}^{\top}V_{t}^{\top}\big{)}$
can be computed efficiently since

$\Big{(}\mathcal{A}\big{(}V_{t}Q_{\underline{c}}\Lambda_{\underline{c}}Q_{\underline{c}}^{\top}V_{t}^{\top}\big{)}\Big{)}_{i}=\mathrm{tr}\Big{(}V_{t}^{\top}Q_{\underline{c}}^{\top}\big{(}A_{i}V_{t}Q_{\underline{c}}\Lambda_{\underline{c}}\big{)}\Big{)}.$
(36)

###  4.1 Nyström Sketch

Consider any iterate $\bar{X}_{t}\in\SS^{n}_{+}$. Let $r\in\mathbb{N}$ be a
parameter that trades of accuracy for scalability. To construct the Nyström
sketch, we sample a random projection matrix $\Omega\in\mathbb{R}^{n\times r}$
such that $\Omega_{ij}\sim\mathcal{N}(0,1)$ are sampled i.i.d. The projection
of $\bar{X}_{t}$, or _sketch_ , can be computed as follows

$P_{t}=\bar{X}_{t}\Omega\in\mathbb{R}^{n\times r}.$ (37)

We can efficiently maintain $P_{t}$ as $\bar{X}_{t}$ evolves with low-rank
updates,

$P_{t+1}\leftarrow\Big{(}\eta_{t+1}\bar{X}_{t}+V_{t}Q_{\underline{c}}\Lambda_{\underline{c}}Q_{\underline{c}}^{\top}V_{t}^{\top}\Big{)}\Omega=\eta_{t+1}P_{t}+V_{t}Q_{\underline{c}}\Lambda_{\underline{c}}\Big{(}Q_{\underline{c}}^{\top}V_{t}^{\top}\Omega\Big{)}.$
(38)

Note that computing $P_{t+1}$ in this manner is equivalent to applying the
projection matrix to $X_{t+1}$. Given the sketch $P_{t}$ and the projection
matrix $\Omega$, we can compute a rank-$r$ approximation to $\bar{X}_{t}$. The
approximation is as follows

$\widehat{\bar{X}_{t}}:=P_{t}(\Omega^{\top}P_{t})^{+}P_{t}^{\top}=(\bar{X}_{t}\Omega)(\Omega^{\top}\bar{X}_{t}\Omega)^{+}(\bar{X}_{t}\Omega)^{\top},$
(39)

where $\vphantom{\cdot}{}^{+}$ is the Moore-Penrose inverse. The
reconstruction is called a _Nyström approximation_. We will almost always
compute the best rank-$r$ approximation of $\widehat{\bar{X}_{t}}$ for memory
efficiency. In practice, we use the numerically stable implementation of the
Nyström approximation provided by Yurtsever et al. [107]. The Nyström
approximation yields a provably good approximation for any sketched matrix.
See the following fact from [87],

######  Fact 4.1.

Fix $X\in\SS^{n}_{+}$ arbitrarily. Let $P:=X\Omega$ where
$\Omega\in\mathbb{R}^{n\times r}$ has independently sampled standard normal
entries. For each $r^{\prime}<r-1$, the Nyström approximation $\widehat{X}$ as
defined in (39) satisfies

$\mathbb{E}_{\Omega}\big{\|}X-\widehat{X}\big{\|}\mathclose{\vphantom{)}}_{*}\leq\left(1+\frac{r^{\prime}}{r-r^{\prime}-1}\right)\big{\|}X-[\mkern-2.5mu[X]\mkern-2.5mu]_{r}\big{\|}\mathclose{\vphantom{)}}_{*},$
(40)

where $\mathbb{E}_{\Omega}$ is the expectation with respect to $\Omega$ and
$[\mkern-2.5mu[X]\mkern-2.5mu]$ is returns an $r$-truncated singular-value
decomposition of the matrix $X$, which is a best rank-r approximation with
respect to every unitarily invariant norm _[64]_. If we replace $\widehat{X}$
with $[\mkern-2.5mu[\widehat{X}]\mkern-2.5mu]_{r}$, this error bound still
remains valid.

This means that by applying this Nyström sketching technique the only
approximation error is in the approximation $\widehat{\bar{X}_{t}}$, and all
of the error is constant (i.e. the error does not compound over time) and
comes only from the projection matrix (and the choice of $r$).

##### Memory requirements.

When we implement SpecBM using this matrix sketching procedure we see
significant decrease in time and memory required by SpecBM. Storing the
projection matrix $\Omega$ and $P_{t}$ requires $\mathcal{O}(nr)$ memory.
Storing $\langle C,\bar{X}_{t}\rangle$ and $\mathrm{tr}(\bar{X}_{t})$ requires
$\mathcal{O}(1)$ memory, respectively, and storing $\mathcal{A}\bar{X}_{t}$,
$\nu_{t}$, $y_{t}$, and $\tilde{y}_{t+1}$ requires $\mathcal{O}(m)$ memory.
The primal-dual interior point methods require $\mathcal{O}(k^{2})$ and
$\mathcal{O}(k^{4})$ memory, respectively, and storing $V_{t}$ requires
$\mathcal{O}(nk)$ memory. This means that the entire memory required by SpecBM
is $\mathcal{O}(nr+nk+m+k^{4})$ working memory (not including the memory to
store the problem data) which for many SDPs is much less than explicitly
storing the iterate $\bar{X}_{t}$ which requires $\mathcal{O}(n^{2})$ memory.

##  5 Experiments

We evaluate SpecBM on three different problem types with and without warm-
starting strategies against the scalable semidefinite programming algorithm
CGAL [106, 107] (with sketching). Yurtsever et al. [107] perform an extensive
evaluation against strong SDP solvers, including SeDuMi [81], SDPT3 [86],
Mosek [9], SDPNAL+[103], and the BM factorization heuristic [21]. They show on
several applications that CGAL with sketching is by far the standard against
which to compare for scalable semidefinite programming. Following [107], we
scale the problem data for all problems such that the following condition
holds,

$\|C\|\mathclose{\vphantom{)}}_{\textrm{F}}=\mathrm{tr}(X_{\star})=1.$ (41)

In all three problem types, the constraints exactly determine the trace of all
optimal solutions. Since we know the trace of the optimal solution and we
apply problem scaling (41), we can set $\alpha=2$ for all problem types.

We consider the primal iterate an $\varepsilon$-approximate solution if the
relative primal objective suboptimality and relative infeasibility is less
than $\varepsilon$, i.e.

$\frac{\langle C,X\rangle-\langle C,X_{\star}\rangle}{1+\langle
C,X\rangle}\leq\varepsilon\quad\textrm{and}\quad\frac{\mathrm{dist}(\mathcal{A}X,\mathcal{K})}{1+\|b\|}\leq\varepsilon.$
(42)

Computing the relative infeasibility is simple. As for the relative primal
objective suboptimality, we do not usually have access to the optimal primal
objective value $\langle C,X_{\star}\rangle$, but we can compute an upper
bound. For any $y\in\mathcal{Y}$,

$\displaystyle
f(y)=\alpha\max\\{\lambda_{\textrm{max}}(C-\mathcal{A}\mathclose{\vphantom{)}}^{*}y),0\\}+\langle
b,y\rangle\geq\langle
C-\mathcal{A}\mathclose{\vphantom{)}}^{*}y,X_{\star}\rangle-\langle
b,y\rangle=\langle C,X_{\star}\rangle,$ (43)
$\displaystyle\quad\implies\langle C,X\rangle-\langle
C,X_{\star}\rangle\leq\langle C,X\rangle-f(y).$

We use this upper bound to approximate the relative primal objective
suboptimality in our experiments for SpecBM. We use the upper bound given in
[107] to compute the relative primal objective suboptimality for CGAL.

For both algorithms, we compute eigenvectors and eigenvalues using an
implementation of thick-restart Lanczos [99, 44] with 32 inner iterations and
a maximum of 10 restarts. Both CGAL and SpecBM are implemented using 64-bit
floating point arithmetic. For all experiments, we use a CPU machine with 16
cores and up to 128 GB of RAM. Throughout the presentation of the experimental
results we use $\uparrow$ to indicate that _higher is better_ for that
particular metric and $\downarrow$ to indicate _lower is better_ for that
particular metric.

###  5.1 MaxCut

The MaxCut problem is a fundamental combinatorial optimization problem and its
SDP relaxation is a common test bed for SDP solvers. Given an undirected
graph, the MaxCut is a partitioning of the $n$ vertices into two sets such
that the number of edges between the two subsets is maximized. Formally, the
MaxCut problem can be written as follows

$\operatorname*{maximize}\quad\frac{1}{4}\,x^{\top}\\!L\,x\quad\operatorname*{subject\;to}\quad
x\in\\{\pm 1\\}^{n},$ (44)

where $L$ is the graph Laplacian. This optimization problem is known to be NP-
hard [48], but rounding the solution to the SDP relaxation can yield very
strong approximate solutions [38].

####  5.1.1 SDP Relaxation

For any vector $x\in\\{\pm 1\\}^{n}$, the matrix $X:=xx^{\top}$ is positive
semidefinite and all of its diagonal entries are exactly one. These implicit
constraints allow us to arrive at the MaxCut SDP as follows

$\operatorname*{maximize}\quad\frac{1}{4}\,\mathrm{tr}(LX)\quad\operatorname*{subject\;to}\quad\textrm{diag}(X)=\mathbf{1}\textrm{
and }X\succeq 0,$ (45)

where $\textrm{diag}(\cdot)$ extracts the diagonal of a matrix into a vector
and $\textbf{1}\in\mathbb{R}^{n}$ is the all-ones vector. A matrix solution
$X_{\star}$ to (45) does not readily elicit a graph cut. One way to round
$X_{\star}$ is to compute a best rank-one approximation
$x_{\star}x_{\star}^{\top}$ and use $\textrm{sgn}(x_{\star})$ as an
approximate cut. There are more sophisticated rounding procedures (e.g. [38]),
but this one works quite well in practice.

Table 1: MaxCut data instance statistics. | fe_sphere | hi2010 | fe_tooth | 598a | 144 | az2010 | auto | rgg_n_2_19_s0 | netherlands_osm | 333SP
---|---|---|---|---|---|---|---|---|---|---
$n$ | 16K | 25K | 78K | 111K | 145K | 242K | 449K | 524K | 2.2M | 3.7M
$\textrm{nnz}(L)$ | 115K | 149K | 983K | 1.6M | 2.3M | 1.4M | 7.1M | 7.1M | 7.1M | 25.9M
Figure 1: Convergence time (sec) to moderate relative error tolerance
($\downarrow$). The time (in seconds) for CGAL and SpecBM to achieve an
$\varepsilon$-approximate solution for $\varepsilon=10^{-1}$ with and without
warm-starting on 99% of initial data on ten DIMACS10 instances. The bars
marked with $\dagger$ indicate that an $\varepsilon$-approximate solution was
not achieved in 72 hours. The datasets are sorted left-to-right in ascending
order by $n$, ranging from 16K to 3.7M (more than $10^{13}$ decision variables
for 333SP). We observe that SpecBM achieves a staggering improvement in
convergence over CGAL which fails to reach an accurate solution on 8 out of 10
problem instances. In contrast, SpecBM is able to reach a solution on all of
the problem instances in 28 hours or less without a warm-start initialization.

Notice that (45) contains exactly $n$ equality constraints. Additionally, in
many instances of MaxCut, the graph Laplacian $L$ is sparse and the optimal
solution is low-rank. This means that we can represent the MaxCut instance
with far less than $\mathcal{O}(n^{2})$ memory. These attributes make
sketching the matrix variable $X$ with relatively small $r$ an attractive and
effective approach when considering the MaxCut problem.

####  5.1.2 Rounding

Both CGAL and SpecBM maintain a sketch of the primal variable that, along with
the test matrix, can be used to construct a low-rank approximation of the
primal iterate $\widehat{X}=U\Lambda U^{\top}$ where $U\in\mathbb{R}^{n\times
r}$ has orthonormal columns. Following [107], we evaluate the size of $r$ cuts
given by the column vectors of $\textrm{sgn}(U)$ and choose the largest.

####  5.1.3 Experimental Setup

We evaluate the CGAL and SpecBM with and without warm-starting on ten
instances from the DIMACS10 [11] where $n$ and the number of edges in each
graph are shown in Table 1. We warm-start each method by dropping the last 1%
of vertices from the graph resulting in a graph with 99% of the vertices. We
pad the solution from the 99%-sized problem with zeros and rescale as
necessary to create the warm-start initialization. We use the following
hyperparameters for SpecBM in our experiments: $r=10$, $\rho=0.01$,
$\beta=0.25$, $k_{c}=10$, $k_{p}=1$. In our experiments, we find that
$k_{p}>0$ does improve performance slightly. See the implementation for more
details.

####  5.1.4 Experimental Results

Figure 1 displays the time (in seconds) taken by CGAL and SpecBM, with and
without warm-starting on 99% of the data instance, to acheive an
$\varepsilon$-approximate solution where $\varepsilon=10^{-1}$. In all cases,
we optimize until an $\varepsilon$-approximate solution is found or 72 hours
elapses (including for creating a warm-start initialization). The bars marked
with $\dagger$ indicate that an $\varepsilon$-approximate solution was not
acheived in 72 hours. As can be seen, SpecBM without warm-starting converges
to an $\varepsilon$-approximate solution faster the CGAL with and without
warm-starting. On 8 out of 10 data instances, CGAL fails to achieve an
$\varepsilon$-approximate solution in 72 hours while SpecBM without warm-
starting reaches an $\varepsilon$-approximate solution for all data instances
within 72 hours. On 2 out of 10 data instances, SpecBM fails to utilize the
warm-start initialization to acheive an $\varepsilon$-approximate soluiton in
under 72 hours. In all other cases, warm-starting either helps or at least
does not hurt the convergence time.

(a)

(b)

(c)

Figure 2: Convergence measures on instance 144. We solve instance 144 from
DIMACS10 and plot the primal objective sub-optimality (objective residual,
$\downarrow$), relative infeasibility (infeasibility gap, $\downarrow$), and
weight of the cut ($\uparrow$) produced by the rounding procedure. The warm-
started runs use 99% of the original data to obtain a warm-start
initialization. We observe that SpecBM is able to more reliably leverage a
warm-start initialization.

Figure 2 plots three different convergence measures over time for a large data
instance, 144, for a total time of 12 hours. For this instance, warm-starting
helps achieve a faster convergence rate of both relative primal objective
suboptimality and relative infeasibility. Additionally, the weight of the cut
produced from the primal iterates generated by SpecBM drops in value less when
warm-starting as compared with CGAL. This is an additional indication that
SpecBM is able to utilize the warm-start initialization better than CGAL.
Finally, observe that SpecBM, with and without warm-starting, produces a
marginally better graph cut than CGAL.

###  5.2 Quadratic Assignment Problem

The quadratic assignment problem (QAP) is a very difficult but fundamental
class of combinatorial optimization problems containing the traveling salesman
problem, max-clique, bandwidth problems, and facilities location problems
among others [59]. SDP relaxations have been shown to facilitate finding good
solutions to large QAPs [109]. The QAP can be formally defined as follows

$\operatorname*{minimize}\quad\mathrm{tr}\left(W\Pi
D\Pi^{\top}\right)\quad\operatorname*{subject\;to}\quad\Pi\textrm{ is an
}\texttt{n}\times\texttt{n}\textrm{ permutation matrix},$ (46)

where $W\in\SS^{\texttt{n}}$ is the weight matrix, $D\in\SS^{\texttt{n}}$ is
the distance matrix, and the goal is to optimize for an assignment $\Pi$ which
aligns $W$ and $D$. The number of $\texttt{n}\times\texttt{n}$ permutation
matrices is $\texttt{n}!$, so a brute-force search becomes quickly intractable
as n grows. Generally, QAP instances with $\texttt{n}>30$ are intractable to
solve exactly. In our experiments, using an SDP relaxation and rounding
procedure, we can obtain good solutions to QAPs where n is between 136 and
198.

####  5.2.1 SDP Relaxation

There are many SDP relaxtions for QAPs, but we consider the one presented by
[107] and inspired by [45, 20] which is formulated as follows

$\displaystyle\operatorname*{minimize}\quad\mathrm{tr}\big{(}(D\otimes
W)\,\textsf{Y}\big{)}$ (47)
$\displaystyle\operatorname*{subject\;to}\;\;\,\mathrm{tr}_{1}(\textsf{Y})=I,\;\mathrm{tr}_{2}(\textsf{Y})=I,\;\mathcal{G}(\textsf{Y})\geq
0,$
$\displaystyle\qquad\qquad\quad\;\;\textrm{vec}(\textsf{B})=\textrm{diag}(\textsf{Y}),\;\textsf{B}\textbf{1}=\textbf{1},\;\textbf{1}^{\top}\textsf{B}=\textbf{1}^{\top},\;\textsf{B}\geq
0,$
$\displaystyle\qquad\qquad\quad\;\;X:=\begin{bmatrix}1&\textrm{vec}(\textsf{B})^{\top}\\\
\textrm{vec}(\textsf{B})&\textsf{Y}\end{bmatrix}\succeq
0,\;\mathrm{tr}(\textsf{Y})=\texttt{n}$

where $\otimes$ denotes the Kronecker product, $\mathrm{tr}_{1}(\cdot)$ and
$\mathrm{tr}_{2}(\cdot)$ denote the partial trace over the first and second
systems of Kronecker product respectively, $\mathcal{G}(\textsf{Y})$ extracts
the entries of Y corresponding to the nonzero entries of $D\otimes W$, and
$\textrm{vec}(\cdot)$ stacks the columns of a matrix one on top of the other
to form a vector. In many cases, one of $D$ or $W$ is sparse (i.e.
$\mathcal{O}(n)$ nonzero entries) resulting in $\mathcal{O}(\texttt{n}^{3})$
total constraints for the SDP.

The primal variable $X$ has dimension
$(\texttt{n}^{2}+1)\times(\texttt{n}^{2}+1)$ and as a result the SDP
relaxation has $\mathcal{O}(\texttt{n}^{4})$ decision variables. Given the
aggressive growth in complexity, most SDP based algorithms have difficulty
operating on instances where $\texttt{n}>50$. To reduce the number of decision
variables, we sketch $X$ with $r=\texttt{n}$, resulting in
$\mathcal{O}(\texttt{n}^{3})$ entries in the sketch and the same number of
constraints in most natural instances. We show that this enables CGAL and
SpecBM to scale to instances where $\texttt{n}=198$ (more than 1.5 billion
decision variables and constraints).

####  5.2.2 Rounding

We adopt the rounding method presented in [107] to convert an approximate
solution of (47) into a permutation matrix. Reconstructing an approximation to
the primal variable from the sketch yields $\widehat{X}=U\Lambda U^{\top}\\!$
where $U$ has shape $(\texttt{n}^{2}+1)\times\texttt{n}$. For each column of
$U$, we discard the first entry of the column vector, reshape the remaining
$\texttt{n}^{2}$ entries into an $\texttt{n}\times\texttt{n}$ matrix, and
project the reshaped matrix onto the set of $\texttt{n}\times\texttt{n}$
permutation matrices using the Hungarian method (also known as Munkres’
assignment algorithm) [52, 65, 47]. This procedure yields a feasible point to
the original QAP (46). To get an upper bound for (46), we perform the rounding
procedure on each column of $U$ and choose the permutation matrix $\Pi$ which
minimizes the objective $\mathrm{tr}\\!\left(W\Pi D\Pi^{\top}\right)$.

####  5.2.3 Experimental Setup

We evaluate CGAL and SpecBM on select large instances from QAPLIB [23] and
TSPLIB [72] ranging in size from 136 to 198. Provided with each of these QAP
instances is the known optimum. For many instances, we find that the quality
of the permutation matrix produced by the rounding procedure does not entirely
correlate with the quality of the iterates with respect to the SDP (47). Thus,
we apply the rounding procedure at every iteration of both CGAL and SpecBM.
Our metric for evaluation is relative gap which is computed as follows

$\texttt{relative gap}=\frac{\texttt{upper bound
obtained}-\texttt{optimum}}{\texttt{optimum}}.$

We report the lowest value so far of relative gap as best relative gap. CGAL
[107] is shown to obtain significantly smaller best relative gap than CSDP
[20] and PATH [108] on most instances, and thus, is a strong baseline.

To create a warm-start initialization, we create a slightly smaller QAP by
dropping the final row and column of both $D$ and $W$ (i.e. solve a size
$\texttt{n}-1$ subproblem of the original instance). We use the solution to
slightly smaller problem to set a warm-start initialization for the original
problem, rescaling and padding with zeros where necessary. We stopped
optimizing after one hour. When warm-starting, we optimize the slightly
smaller problem for one hour and then optimize the original problem for one
hour.

Following [107], we apply the following scaling of the problem variables in
(47)

$\|C\|\mathclose{\vphantom{)}}_{\textrm{F}}=\mathrm{tr}(X_{\star})=\|\mathcal{A}\|\mathclose{\vphantom{)}}_{\textrm{op}}=1\quad\textrm{and}\quad\|A_{1}\|\mathclose{\vphantom{)}}_{\textrm{F}}=\cdots=\|A_{m}\|\mathclose{\vphantom{)}}_{\textrm{F}}.$
(48)

We use the following hyperparameters for SpecBM in our experiments:
$\rho=0.005$, $\beta=0.25$, $k_{c}=2$, $k_{p}=0$. In our experiments, we find
that $k_{p}>0$ does not improve performance. See the implementation for more
details.

####  5.2.4 Experimental Results

(a) pr144

(b) kroA150

(c) tai150b

Figure 3: relative gap ($\downarrow$) vs. time. We plot the relative gap
(y-axis) against time in seconds (x-axis) for three instances from QAPLIB and
TSPLIB. We observe that for both algorithms the best rounded solution is found
early in optimization. In addition, we observe that SpecBM is able to more
reliably leverage a warm-start initialization.

(a) pr144

(b) kroA150

(c) tai150b

Figure 4: best relative gap ($\downarrow$) vs. time. We plot the best relative
gap (y-axis) against time in seconds (x-axis) for three instances from QAPLIB
and TSPLIB. We observe that SpecBM is able to more reliably leverage a warm-
start initialization. Figure 5: best relative gap ($\downarrow$). The best
relative gap obtained in one hour of optimization is shown for ten instances
from QAPLIB and TSPLIB. We observe that on most problem instance that SpecBM
produces a better relative gap than CGAL and that a warm-start initialization
helps SpecBM obtain a better relative gap. CGAL is much less reliable in
leveraging a warm-start initialization.

Figure 3 and Figure 4 plot relative gap and best relative gap against time (in
seconds) for three instances, respectively. Figure 3 shows that even for large
instances the best relative gap is obtained within the first few minutes of
the hour of optimization. In these instances, we can also observe that SpecBM
not only obtains a better best relative gap, but also is able to leverage a
warm-start solution to improve performance. Figure 5 shows the best relative
gap obtained for ten data instances over one hour of optimization. It can be
seen that warm-starting does not always yield a better best relative gap, but
SpecBM is better able to leverage a warm-start solution than CGAL. These
results show that warm-starting is a beneficial heuristic for finding a
quality approximate solution to QAPs. The warm-starting technique used in
these experiments is not the only possible warm-starting initialization
strategy. For examples, one could generate several warm-start initializations
for (47) by dropping other rows and columns of $D$ and $W$ besides the last
row and column. Then, after creating n warm-start initializations, we could
take the best upper bound obtained by optimizing (47) from each of those warm-
start initializations.

###  5.3 Interactive Entity Resolution with $\exists$-constraints

Angell et al. [7] introduced a novel SDP relaxation of a combinatorial
optimization problem that arises in automatic knowledge base construction
problem [89, 51, 46]. They consider the problem of interactive entity
resolution with user feedback to correct predictions made by a machine
learning algorithm. Entity resolution (also known as coreference resolution or
record linkage) is the process of identifying which mentions (sometimes
referred to as records) of entities refer to the same entity [16, 6, 1, 2,
102]. Entity resolution is an important problem central to automated knowledge
base construction where the correctness of the resulting knowledge base is
vital to its usefulness. The entity resolution problem is fundamentally a
clustering problem where the prefect clusters contain all mentions of exactly
one entity. Due to the large volumes of raw data frequently involved, machine
learning approaches are often used to perform entity resolution. To correct
inevitable errors in the predictions made by these automated methods, human
feedback can be used to correct errors in entity resolution.

Utilizing human feedback to guide and correct clustering decisions can be
broadly categorized into two groups [12]: active clustering (also know as
semi-supervised clustering) [94, 93, 62, 75, 101] and interactive clustering
[24, 13, 25, 32]333Commonly, _active_ clustering is used to describe
techniques where the algorithm queries the human for feedback about a specific
pair or group of points or clusters (akin to active learning) and
_interactive_ clustering is used to describe techniques where humans observe
the clustering output and provide unelicited feedback, typically in the form
of some type of constraint, to the algorithm.. There are several paradigms of
human feedback used to steer a clustering algorithm [35, 96, 77, 50, 54, 60,
57, 58, 100], the most common of which is pairwise constraints [95] (i.e.
statements about whether two points or mentions must or cannot be clustered
together). Recently, Angell et al. [7] introduced _$\exists$ -constraints _,
a new paradigm of interactive feedback for correcting entity resolution
predictions and presented a novel SDP relaxation as part of a heuristic
algorithm for satisfying $\exists$-constraints in the predicted clustering.
This novel form of feedback allows users to specify constraints stating the
existence of an entity with and without certain features.

####  5.3.1 $\exists$-constraints

To define $\exists$-constraints precisely, we must introduce some terminology
and notation. Let the set of mentions be represented by the vertex set $V$ of
a fully connected graph. We assume that each of the mentions $v\in V$ has an
associated set of discrete features (attributes, relations, other features,
etc.) that is readily available or can be extracted easily using some
automated method. Let $\Phi=\\{\phi_{1},\phi_{2},\ldots\\}$ be the set of
possible features and $\Phi(v)\subseteq\Phi$ be the subset of features
associated with vertex $v\in V$. A clustering
$\mathcal{C}=\\{C_{1},C_{2},\ldots,C_{m}\\}$ of vertices $V$ is a set of
nonempty sets such that $C_{i}\subseteq V$, $C_{i}\cap C_{j}=\emptyset$ for
all $i\neq j$, and $\bigcup_{C\in\mathcal{C}}C=V$. We denote the ground truth
clustering of the vertices (partitioning mentions into entities) as
$\mathcal{C}^{\star}$. Furthermore, for any subset of vertices $S\subseteq V$,
we define the set of features associated with $S$ as $\Phi(S)=\bigcup_{v\in
S}\Phi(v)$. For any cluster $C_{i}$, we use $\Phi(C_{i})$ as the
canonicalization function for a predicted entity, but this framework allows
for more complex (and even learned) canonicalization functions. When providing
feedback, users are able to view the features of predicted clusters and then
provide one or more $\exists$-constraints .

Formally, a $\exists$-constraint $\xi\subseteq\\{+,-\\}\times\Phi$ uniquely
characterizes a constraint asserting there exists a cluster
$C\in\mathcal{C}^{\star}$ such that all of the positive features
$\xi^{+}:=\\{\phi\in\Phi:(+,\phi)\in\xi\\}$ are contained in the features of
$C$ (i.e. $\xi^{+}\subseteq\Phi(C)$) and none of the negative features
$\xi^{-}:=\\{\phi\in\Phi:(-,\phi)\in\xi\\}$ are contained in the features of
$C$ (i.e. $\xi^{-}\cap\Phi(C)=\emptyset$). At any given time, let the set of
$\exists$-constraints be represented by $\Xi$. We say that a subset of nodes
$S\subseteq V$ _satisfies_ a $\exists$-constraint $\xi$ if
$\xi^{+}\subseteq\Phi(S)$ and $\xi^{-}\cap\Phi(S)=\emptyset$. Note that more
than one ground truth cluster can satisfy a $\exists$-constraint . We say
that a subset of vertices $S\subseteq V$ is _incompatible_ with a
$\exists$-constraint $\xi$ if $\Phi(S)\cap\xi^{-}\neq\emptyset$ and denote
this as $S\perp\xi$. Similarly, two $\exists$-constraints $\xi_{a},\xi_{b}$
are incompatible if $\xi_{a}^{+}\cap\xi_{b}^{-}\neq\emptyset$ or
$\xi_{a}^{-}\cap\xi_{b}^{+}\neq\emptyset$, and is also denoted
$\xi_{a}\perp\xi_{b}$. Furthermore, we say a subset of vertices $S\subseteq V$
and a $\exists$-constraint $\xi$ or two $\exists$-constraints
$\xi_{a},\xi_{b}$ are _compatible_ if they are not incompatible and denote
this as $S\not\perp\xi$ and $\xi_{a}\not\perp\xi_{b}$, respectively. Moreover,
we use _compatible_ to describe when a subset of vertices $S\subseteq V$ could
be part of a cluster that _satisfies_ a certain $\exists$-constraint $\xi$.
This definition allows for $\Phi(S)\cap\xi^{+}$ to be empty, but still enables
$S$ to be compatible with $\xi$ just as long as $\Phi(S)\cap\xi^{-}$ is empty.
Observe that if every mention has a unique feature, the $\exists$-constraint
framework fully encompasses must-link and cannot-link constraints.

####  5.3.2 SDP relaxation

Following [7], the general approach to clustering with $\exists$-constraints
is to jointly cluster the vertices and the $\exists$-constraints , so we
include $\exists$-constraints as additional vertices of the graph $V$ (i.e.
$\Xi\subset V$) and create corresponding decision variables for each of the
$\exists$-constraints . Additionally, we create constraints to ensure that
all of the $\exists$-constraints are satisfied and none of the
$\exists$-constraints are incompatible with any of the other members of the
predicted clustering. Given a $\exists$-constraint $\xi$ and positive feature
$\phi\in\xi^{+}$, let $\Gamma(\xi,\phi):=\\{v\in
V\setminus\Xi:\phi\in\Phi(v)\land\Phi(v)\cap\xi^{-}=\emptyset\\}$ be the set
of candidate vertices that could satisfy $\phi\in\xi^{+}$. The resulting
integer linear programming problem is as follows

$\operatorname*{maximize}_{x_{u},x_{v}\in\\{e_{1},e_{2},\ldots,e_{n}\\}}\sum_{(u,v)\in
V\times
V}\\!\\!w_{uv}\,x_{u}^{\top}x_{v}\;\;\operatorname*{subject\;to}\;\;\begin{cases}\sum_{v\in\Gamma(\xi,\phi)}\;x_{v}^{\top}x_{\xi}\geq
1,&\textrm{for all }\phi\in\xi^{+}\\\ x_{v}^{\top}x_{\xi}=0,&\textrm{for all
}v\perp\xi\end{cases}$ (49)

where $w_{uv}\in\mathbb{R}$ is the pairwise similarity between mentions $u$
and $v$ (usually computed using a trained machine learning model trained on
pairs of mentions) and $e_{i}\in\mathbb{R}^{n}$ is the $i^{\textrm{th}}$
standard basis vector. Note that we assume that the $w_{v\xi}=0$ as in [7].
The first constraint ensures that all of the positive features in $\xi^{+}$
are contained in the same cluster as $\xi$ for all $\xi\in\Xi$. The last
constraint ensures that $\xi$ is not in the same cluster as anything
incompatible with it for all $\xi\in\Xi$. This problem is intractable to
optimize exactly, so we relax (49) to an SDP in the standard way as follows

$\operatorname*{maximize}_{X\succeq 0}\;\;\langle
W,X\rangle\;\;\operatorname*{subject\;to}\;\;\begin{cases}X_{vv}=1,&\textrm{for
all }v\in V\\\ X_{uv}\geq 0,&\textrm{for all }u,v\in V\\\
\sum_{v\in\Gamma(\xi,\phi)}\;X_{v\xi}\geq 1,&\textrm{for all
}\phi\in\xi^{+}\\\ X_{v\xi}=0,&\textrm{for all }v\perp\xi\end{cases}$ (50)

The decision variables are changed from standard basis vectors, whose dot
products determine whether or not two elements (vertices or
$\exists$-constraints ) are in the same cluster, to a positive semidefinite
matrix where each entry represents the global affinity two elements have for
one another (corresponding to the row and column of the entry). We constrain
the diagonal of this positive semidefinite matrix to be all ones (representing
each element’s affinity with itself) and we enforce that all entries in the
matrix are non-negative. The constraints ensuring $\exists$-constraint
satisfaction are similar to (49). Let $W$ be the (weighted) adjacency matrix
of the fully connected graph over the $n$ vertices.

####  5.3.3 Rounding

The solution $X_{\star}$ to (50) is a fractional approximation of the problem
we wish to solve. In order to infer a predicted clustering, we need some
procedure to round the fractional solution to an integer solution. Angell et
al. [7] propose first building a hierarchical clustering over $V$ using
$X_{\star}$ as the similarity function. Then, a dynamic programming algorithm
is used to extract a (predicted) flat clustering from the hierarchical
clustering which maximizes the number $\exists$-constraints satisfied and the
(correlation) clustering objective. For a more detailed description of the
rounding procedure, see [7]. The predicted clustering is then observed by the
human users in terms of the features of each predicted cluster enabling the
users to provide any additional $\exists$-constraints to correct the entity
resolution decisions. We simulate this human feedback loop with an oracle
$\exists$-constraint generator.

####  5.3.4 Experimental Setup

We evaluate the performance of the CGAL and SpecBM with and without warm-
starting on the same three author coreference datasets used in [7]. In these
datasets, each mention is an author-paper pair and the goal is to identify
which author-paper pairs were written by the same author. As standard in
author coreference, the datasets are preprocessed into _blocks_ (also known as
_canopies_) based on the author’s first name initial and last name (e.g.
authors “Rajarshi Das” and “Ravi Das” would both be contained in the same
block “r das”, but in a different block than “Jane Smith”, which would be in
“j smith”). Within each block, similarity scores are computed between all
pairs of mentions using a trained pairwise model [7, 83]. We then perform two
additional preprocessing steps to convert the many small dense problems into
one large sparse problem. We aggregate all of these pairwise similarity scores
into a block diagonal weight matrix and fill the remaining entries with
negative one. We then sparsify this highly structured similarity matrix using
a spectral sparsifier [79]. Table 2 details the dataset statistics including
the number of non-zeros in the pairwise similarity matrix after spectral
sparsification, denoted nnz($W$).

Table 2: Author coreference dataset statistics. | # mentions | # blocks | # clusters |  nnz($W$) | # features
---|---|---|---|---|---
PubMed | 315 | 5 | 34 | 3,973 | 14,093
QIAN | 410 | 38 | 77 | 5,158 | 10,366
SCAD-zbMATH | 1,196 | 120 | 166 | 18,608 | 8,203

We simulate $\exists$-constraint generation using the same oracle implemented
in [7]. The oracle has access to the features of the ground truth clusters and
the features of the predicted clusters output by the rounding algorithm and
generates a small $\exists$-constraint which is satisfied by the ground truth
clustering, but is not satisfied by the predicted clustering. The oracle
generates $\exists$-constraints one at a time and added to the optimization
problem until the ground truth clustering is predicted. For both solvers, we
iterate until the relative error tolerance (42) and the following absolute
infeasibility condition are satisfied

$\|\mathcal{A}X-\mathrm{proj}_{\mathcal{K}}(\mathcal{A}X)\|_{\infty}\leq\varepsilon.$
(51)

We use $\varepsilon=10^{-1}$ to determine the stopping condition and find that
this is sufficient for the rounding algorithm to satisfy nearly all of the
$\exists$-constraints generated by the oracle. We also note that the
$\ell^{\infty}$-norm condition is especially important in this application in
order to make sure the relaxed $\exists$-constraints are satisfied. Due to the
small problem size, we simplify the implementation by using variants of CGAL
and SpecBM without sketching.

When warm-starting, we predict the values of the expanded primal variable
corresponding to the newly added $\exists$-constraint by computing
representations of each vertex using a low-rank factorization, performing a
weighted average of the representations of the positive vertices in the newly
added $\exists$-constraint , and expanding the factorization into an
initialization of primal variable for the new problem. We pad the dual
variable with zeros. Since there is no pairwise similarity between
$\exists$-constraints and mentions, we find that scaling up the norms of
$A_{i}$’s corresponding to the constraints $X_{vv}=1$ for all $v$ such that
$\Phi(v)\cap\xi_{\textrm{new}}\neq\emptyset$, where $\xi_{\textrm{new}}$ is
the newly added $\exists$-constraint . We use the following SpecBM
hyperparameters in our experiments: $\rho=0.01$, $\beta=0.25$, $k_{c}=3$,
$k_{p}=0$. See the implementation for more details.

####  5.3.5 Experimental Results

(a) PubMed

(b) QIAN

(c) SCAD-zbMATH

Figure 6: Cumulative SDP solve time ($\downarrow$) vs. number of
$\exists$-constraints . In each plot, the $x$-axis is the number of
$\exists$-constraints generated one after the other over time and the $y$-axis
is the the cumulative solve time (in seconds) for each SDP solver to reach a
relative suboptimality, relative infeasibility, and max absolute infeasibility
of $\varepsilon=10^{-1}$. When warm-starting, both solvers are initialized
using the solution from the previous SDP (with one less $\exists$-constraint
). $\exists$-constraints are generated until the perfect clustering is
predicted. We observe that SpecBM is able to much more reliably leverage a
warm-start initialization. In addition, we observe that the performance gap
between SpecBM and CGAL grows as the problem size grows. Figure 7: Average
warm-start fold change ($\uparrow$). Warm-start fold change represents the
ratio of SDP solve time without warm-starting divided by the SDP solve time
with warm-starting per $\exists$-constraint . The error bars indicate one
standard deviation from the mean. We observe that SpecBM is able to much more
reliably leverage a warm-start initialization. In addition, we observe that
the performance gap between SpecBM and CGAL grows as the problem size grows.

Figure 6 shows the cumulative SDP solve time against the number of
$\exists$-constraints . As shown in the plots, warm-starting enables SpecBM
to converge an $\varepsilon$-approximate solution 2-3x faster than without
warm-starting over time. In contrast, warm-starting actually slows the
convergence of CGAL significantly. We also notice that the larger the problem
size, the bigger the gap in solve times between SpecBM and CGAL without warm-
starting.

Figure 7 shows the average warm-start SDP solve time fold change. Warm-start
fold change greater than one indicates a speedup in solve time while warm-
start fold change less than one indicates a slowdown in solve time. We see
that the warm-start fold change for CGAL is less than or equal to one,
indicating a slowdown consistent with Figure 6. We see that on average warm-
starting affords SpecBM a 20-100x speedup on average per $\exists$-constraint
. Note that this is much faster than the 2-3x we see in Figure 6. This is
due to the fact that sometimes warm-starting does not help SpecBM, which makes
intuitive sense if the warm-start initialization is far away from the solution
set for the new SDP. We believe that with some further experimentation
including additional pairwise learned similarities between mentions and
$\exists$-constraints or different warm-starting strategies we might be able
to mitigate these situations. Regardless, in most cases, warm-starting
provides a significant improvement in convergence time when using SpecBM. We
also note that Figure 7 shows that warm-starting provides more of a benefit
the larger the SDP we are trying to solve.

##  6 Related Work

##### Efficient SDP solving.

Due to their wide-spread applicability and importance, methods for solving
SDPs efficiently have been extensively studied in the literature [84, 67, 68,
3, 18, 34, 106, 104, 29]. Interior point methods [40, 68, 3] are by far the
most widely-used and well-known method for solving SDPs. Interior point
methods enjoy quadratic convergence while having incredibly high per-iteration
complexity, and thus, are impractical for instances having more than a few
hundred variables. ADMM [18] and splitting-based methods [69] have also been
applied to solving SDPs. These methods are generally more applicable than
interior point methods, but our still restricted do to the computationally
expensive full eigenvalue decomposition. First-order methods such as
conditional gradient augmented lagrangian [105], primal-dual subgradient
algorithms [104], the matrix multiplicative weight method [10, 88], and the
mirror-prox algorithm with the quantum entropy mirror map [66] have by far the
lowest per-iteration complexity among all of the methods for solving SDPs.
Along with this lower per-iteration complexity comes a slower convergence rate
for these first order methods.

##### Spectral bundle methods.

Helmberg and Rendl [42] were the first to introduce spectral bundle methods
for equality-constrained SDPs. The biggest difference between the original
spectral bundle method and our approach is the model where $k_{c}=1$ is fixed
and $k_{p}$ is chosen by the user. We find that larger $k_{c}$ provides better
convergence in general and $k_{p}$ is not very helpful (and sometimes hurtful)
to convergence. Several other algorithmic variants of spectral bundle methods
exist including solving inequality-constrained SDPs [41], trust region-based
methods for nonconvex eigenvalue optimization [8], and incorporating second-
order information to speed up empirical convergence [43]. Ding and Grimmer
[28] present a spectral bundle method most similar to ours for solving
equality-constrained SDPs. They prove an impressive set of convergence rates
under various assumptions, including linear convergence if certain strong
assumptions hold. The authors also show that linear convergence is exhibited
in their experiments.

##### Memory-efficient semidefinte programming.

The most famous memory efficient approach to semidefinite programming is the
Burer–Monteiro (BM) factorization heuristic [21]. Several optimization
techniques have been used to optimize the low-rank factorized problem [98, 17,
22, 53, 74, 78]. Ding et al. [30] present an optimal storage approach to
solving SDPs by first approximately solving the dual problem and then using
the dual slack matrix to solve the primal problem efficiently. Yurtsever et
al. [107] implement the conditional gradient augmented Lagrangian method [106]
with the matrix sketching technique in [105] which forms a very strong state-
of-the-art method for scalably solving general SDPs. Ding and Grimmer [28]
implements a spectral bundle method for solving equality constrained SDPs with
the same matrix sketching technique [105].

##  7 Conclusion

In this work, we presented a practical spectral bundle method for solving
large SDPs with both equality and inequality constraints. We proved non-
asymptotic convergence rates under standard assumptions for SpecBM. We showed
empirically that SpecBM is fast and scalable on three practical types of SDP
instances. Additionally, we showed that SpecBM can more reliably leverage a
warm-start initialization to accelerate convergence. Lastly, we make our
standalone implementation in pure JAX available for wide-spread use.

Given these results, we believe there is promise in the following future
directions. First, the SpecBM algorithm presented in this work is limited to
equality and inequality constraints, but we have provided a framework via the
dual slack variable to incorporate other types of constraints (e.g. norm
constraints). Second, the success of warm-starting shown empirically in this
work indicates that SpecBM might be useful as a part of a solution to solving
mixed-integer SDPs when combined with a branch-and-cut approach. Lastly, it
might be useful to analyze theoretically when warm-starting should or should
not be beneficial to speeding up convergence. For example, can we run a small
test to check and see whether or not warm-starting is better than starting
from a cold-start variable initialization.

## Acknowledgements

This work is funded in part by the Center for Data Science and the Center for
Intelligent Information Retrieval at the University of Massachusetts Amherst,
and in part by the National Science Foundation under grants IIS-1763618 and
IIS-1922090, and in part by the Chan Zuckerberg Initiative under the project
Scientific Knowledge Base Construction, and in part by IBM Research AI through
the AI Horizons Network. The work reported here was performed in part by the
Center for Data Science and the Center for Intelligent Information Retrieval,
and in part using high performance computing equipment obtained under a grant
from the Collaborative R&D Fund managed by the Massachusetts Technology
Collaborative. Rico Angell is partially supported by the National Science
Foundation Graduate Research Fellowship under grant 1938059. Any opinions,
findings and conclusions or recommendations expressed in this material are
those of the authors and do not necessarily reflect those of the sponsor(s).

## References

  * Agarwal et al. [2021] Dhruv Agarwal, Rico Angell, Nicholas Monath, and Andrew McCallum.  Entity linking and discovery via arborescence-based supervised clustering.  _arXiv preprint arXiv:2109.01242_ , 2021.
  * Agarwal et al. [2022] Dhruv Agarwal, Rico Angell, Nicholas Monath, and Andrew McCallum.  Entity linking via explicit mention-mention coreference modeling.  In _Proceedings of the 2022 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies_ , 2022.
  * Alizadeh [1995] Farid Alizadeh.  Interior point methods in semidefinite programming with applications to combinatorial optimization.  _SIAM journal on Optimization_ , 1995.
  * Alizadeh et al. [1997] Farid Alizadeh, Jean-Pierre A Haeberly, and Michael L Overton.  Complementarity and nondegeneracy in semidefinite programming.  _Mathematical programming_ , 1997.
  * Alizadeh et al. [1998] Farid Alizadeh, Jean-Pierre A Haeberly, and Michael L Overton.  Primal-dual interior-point methods for semidefinite programming: convergence rates, stability and numerical results.  _SIAM Journal on Optimization_ , 1998.
  * Angell et al. [2020] Rico Angell, Nicholas Monath, Sunil Mohan, Nishant Yadav, and Andrew McCallum.  Clustering-based inference for biomedical entity linking.  _arXiv preprint arXiv:2010.11253_ , 2020.
  * Angell et al. [2022] Rico Angell, Nicholas Monath, Nishant Yadav, and Andrew Mccallum.  Interactive correlation clustering with existential cluster constraints.  In _Proceedings of the 39th International Conference on Machine Learning_ , 2022.
  * Apkarian et al. [2008] Pierre Apkarian, Dominikus Noll, and Olivier Prot.  A trust region spectral bundle method for nonconvex eigenvalue optimization.  _SIAM Journal on Optimization_ , 2008.
  * ApS [2019] Mosek ApS.  Mosek optimization toolbox for matlab.  _User’s Guide and Reference Manual, Version_ , 2019.
  * Arora et al. [2005] Sanjeev Arora, Elad Hazan, and Satyen Kale.  Fast algorithms for approximate semidefinite programming using the multiplicative weights update method.  In _46th Annual IEEE Symposium on Foundations of Computer Science (FOCS’05)_ , 2005.
  * [11] David A. Bader, Henning Meyerhenke, Peter Sanders, and Dorothea Wagner.  Graph partitioning and graph clustering.  _10th DIMACS Implementation Challenge Workshop. February 13-14, 2012._ URL https://www.cise.ufl.edu/research/sparse/matrices/DIMACS10/index.html.
  * Bae et al. [2020] Juhee Bae, Tove Helldin, Maria Riveiro, Sławomir Nowaczyk, Mohamed-Rafik Bouguelia, and Göran Falkman.  Interactive clustering: A comprehensive review.  _ACM Computing Surveys (CSUR)_ , 2020.
  * Basu et al. [2010] Sumit Basu, Danyel Fisher, Steven Drucker, and Hao Lu.  Assisting users with clustering tasks by combining metric learning and classification.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 24, 2010.
  * Beck [2015] Amir Beck.  On the convergence of alternating minimization for convex programming with applications to iteratively reweighted least squares and decomposition schemes.  _SIAM Journal on Optimization_ , 2015.
  * Bertsimas [1995] Dimitris Bertsimas.  The achievable region method in the optimal control of queueing systems; formulations, bounds and policies.  _Queueing systems_ , 1995.
  * Binette and Steorts [2022] Olivier Binette and Rebecca C Steorts.  (Almost) all of entity resolution.  _Science Advances_ , 2022.
  * Boumal et al. [2014] Nicolas Boumal, Bamdev Mishra, P-A Absil, and Rodolphe Sepulchre.  Manopt, a matlab toolbox for optimization on manifolds.  _The Journal of Machine Learning Research_ , 2014.
  * Boyd et al. [2011] Stephen Boyd, Neal Parikh, Eric Chu, Borja Peleato, Jonathan Eckstein, et al.  Distributed optimization and statistical learning via the alternating direction method of multipliers.  _Foundations and Trends® in Machine learning_ , 2011\.
  * Bradbury et al. [2018] James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne, et al.  Jax: composable transformations of python+ numpy programs.  2018\.
  * Bravo Ferreira et al. [2018] José FS Bravo Ferreira, Yuehaw Khoo, and Amit Singer.  Semidefinite programming approach for the quadratic assignment problem with a sparse graph.  _Computational Optimization and Applications_ , 2018.
  * Burer and Monteiro [2003] Samuel Burer and Renato DC Monteiro.  A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization.  _Mathematical programming_ , 2003.
  * Burer and Monteiro [2005] Samuel Burer and Renato DC Monteiro.  Local minima and convergence in low-rank semidefinite programming.  _Mathematical programming_ , 2005.
  * Burkard et al. [1997] Rainer E Burkard, Stefan E Karisch, and Franz Rendl.  Qaplib–a quadratic assignment problem library.  _Journal of Global optimization_ , 1997.
  * Chuang and Hsu [2014] Jason Chuang and Daniel J Hsu.  Human-centered interactive clustering for data analysis.  In _Conference on Neural Information Processing Systems (NIPS). Workshop on Human-Propelled Machine Learning_ , 2014.
  * Coden et al. [2017] Anni Coden, Marina Danilevsky, Daniel Gruhl, Linda Kato, and Meena Nagarajan.  A method to accelerate human in the loop clustering.  In _Proceedings of the 2017 SIAM International Conference on Data Mining_ , pages 237–245. SIAM, 2017.
  * Dathathri et al. [2020] Sumanth Dathathri, Krishnamurthy Dvijotham, Alexey Kurakin, Aditi Raghunathan, Jonathan Uesato, Rudy R Bunel, Shreya Shankar, Jacob Steinhardt, Ian Goodfellow, Percy S Liang, et al.  Enabling certification of verification-agnostic networks via memory-efficient semidefinite programming.  _Advances in Neural Information Processing Systems_ , 2020.
  * Díaz and Grimmer [2023] Mateo Díaz and Benjamin Grimmer.  Optimal convergence rates for the proximal bundle method.  _SIAM Journal on Optimization_ , 2023.
  * Ding and Grimmer [2023] Lijun Ding and Benjamin Grimmer.  Revisiting spectral bundle methods: Primal-dual (sub) linear convergence rates.  _SIAM Journal on Optimization_ , 2023.
  * Ding and Udell [2021] Lijun Ding and Madeleine Udell.  On the simplicity and conditioning of low rank semidefinite programs.  _SIAM Journal on Optimization_ , 2021.
  * Ding et al. [2021] 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_ , 2021.
  * Drusvyatskiy and Wolkowicz [2017] Dmitriy Drusvyatskiy and Henry Wolkowicz.  The many faces of degeneracy in conic optimization.  _Foundations and Trends® in Optimization_ , 2017.
  * Dubey et al. [2010] Avinava Dubey, Indrajit Bhattacharya, and Shantanu Godbole.  A cluster-level semi-supervision model for interactive clustering.  In _Joint European Conference on Machine Learning and Knowledge Discovery in Databases_. Springer, 2010.
  * Feltenmark and Kiwiel [2000] Stefan Feltenmark and Krzysztof C Kiwiel.  Dual applications of proximal bundle methods, including lagrangian relaxation of nonconvex problems.  _SIAM Journal on Optimization_ , 2000.
  * Friedlander and Macedo [2016] Michael P Friedlander and Ives Macedo.  Low-rank spectral optimization via gauge duality.  _SIAM Journal on Scientific Computing_ , 2016.
  * Frome et al. [2007] Andrea Frome, Yoram Singer, Fei Sha, and Jitendra Malik.  Learning globally-consistent local distance functions for shape-based image retrieval and classification.  In _2007 IEEE 11th International Conference on Computer Vision_ , 2007\.
  * Frostig et al. [2018] Roy Frostig, Matthew James Johnson, and Chris Leary.  Compiling machine learning programs via high-level tracing.  _Systems for Machine Learning_ , 2018.
  * Gittens [2013] Alex Gittens.  _Topics in randomized numerical linear algebra_.  California Institute of Technology, 2013.
  * Goemans and Williamson [1995] Michel X Goemans and David P Williamson.  Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming.  _Journal of the ACM (JACM)_ , 1995.
  * Halko et al. [2011] Nathan Halko, Per-Gunnar Martinsson, and Joel A Tropp.  Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions.  _SIAM review_ , 2011.
  * Helmberg [1994] Christoph Helmberg.  An interior point method for semidefinite programming and max-cut bounds.  1994\.
  * Helmberg and Kiwiel [2002] Christoph Helmberg and Krzysztof C Kiwiel.  A spectral bundle method with bounds.  _Mathematical Programming_ , 2002.
  * Helmberg and Rendl [2000] Christoph Helmberg and Franz Rendl.  A spectral bundle method for semidefinite programming.  _SIAM Journal on Optimization_ , 2000.
  * Helmberg et al. [2014] Christoph Helmberg, Michael L Overton, and Franz Rendl.  The spectral bundle method with second-order information.  _Optimization Methods and Software_ , 2014.
  * Hernández et al. [2007] Vicente Hernández, Jose E Román, Andrés Tomás, and Vicente Vidal.  Krylov-schur methods in slepc.  _Universitat Politecnica de Valencia, Tech. Rep. STR-7_ , 2007.
  * Huang et al. [2014] Qixing Huang, Yuxin Chen, and Leonidas Guibas.  Scalable semidefinite relaxation for maximum a posterior estimation.  In _International Conference on Machine Learning_. PMLR, 2014.
  * Iv et al. [2022] Robert Iv, Alexandre Passos, Sameer Singh, and Ming-Wei Chang.  Fruit: Faithfully reflecting updated information in text.  In _Proceedings of the 2022 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies_ , 2022.
  * Jonker and Volgenant [1988] Roy Jonker and Ton Volgenant.  A shortest augmenting path algorithm for dense and sparse linear assignment problems.  In _DGOR/NSOR: Papers of the 16th Annual Meeting of DGOR in Cooperation with NSOR/Vorträge der 16. Jahrestagung der DGOR zusammen mit der NSOR_. Springer, 1988.
  * Karp [1972] Richard M. Karp.  _Reducibility among Combinatorial Problems_ , pages 85–103.  Springer US, 1972.
  * Kiwiel [2000] Krzysztof C Kiwiel.  Efficiency of proximal bundle methods.  _Journal of Optimization Theory and Applications_ , 2000.
  * Klein et al. [2002] Dan Klein, Sepandar D Kamvar, and Christopher D Manning.  From instance-level constraints to space-level constraints: Making the most of prior knowledge in data clustering.  Technical report, 2002.
  * Krishnamurthy [2015] Jayant Krishnamurthy.  Learning to understand natural language with less human effort.  Technical report, CARNEGIE-MELLON UNIV PITTSBURGH PA SCHOOL OF COMPUTER SCIENCE, 2015.
  * Kuhn [1955] Harold W Kuhn.  The hungarian method for the assignment problem.  _Naval research logistics quarterly_ , 1955.
  * Kulis et al. [2007] Brian Kulis, Arun C Surendran, and John C Platt.  Fast low-rank semidefinite programming for embedding and clustering.  In _Artificial Intelligence and Statistics_. PMLR, 2007.
  * Kulis et al. [2009] Brian Kulis, Sugato Basu, Inderjit Dhillon, and Raymond Mooney.  Semi-supervised graph clustering: a kernel approach.  _Machine learning_ , 2009.
  * Lemarechal et al. [1981] Claude Lemarechal, Jean-Jacques Strodiot, and André Bihain.  On a bundle algorithm for nonsmooth optimization.  In _Nonlinear programming 4_. Elsevier, 1981.
  * Li et al. [2017] Huamin Li, George C Linderman, Arthur Szlam, Kelly P Stanton, Yuval Kluger, and Mark Tygert.  Algorithm 971: An implementation of a randomized algorithm for principal component analysis.  _ACM Transactions on Mathematical Software (TOMS)_ , 2017.
  * Li et al. [2008] Zhenguo Li, Jianzhuang Liu, and Xiaoou Tang.  Pairwise constraint propagation by semidefinite programming for semi-supervised classification.  In _Proceedings of the 25th international conference on Machine learning_ , 2008.
  * Li et al. [2009] Zhenguo Li, Jianzhuang Liu, and Xiaoou Tang.  Constrained clustering via spectral regularization.  In _2009 IEEE Conference on Computer Vision and Pattern Recognition_ , 2009.
  * Loiola et al. [2007] Eliane Maria Loiola, Nair Maria Maia De Abreu, Paulo Oswaldo Boaventura-Netto, Peter Hahn, and Tania Querido.  A survey for the quadratic assignment problem.  _European journal of operational research_ , 2007.
  * Lu and Carreira-Perpinan [2008] Zhengdong Lu and Miguel A Carreira-Perpinan.  Constrained spectral clustering through affinity propagation.  In _2008 IEEE Conference on Computer Vision and Pattern Recognition_ , pages 1–8. IEEE, 2008.
  * Luo and Tseng [1993] Zhi-Quan Luo and Paul Tseng.  Error bounds and convergence analysis of feasible descent methods: a general approach.  _Annals of Operations Research_ , 1993.
  * Mazumdar and Saha [2017] Arya Mazumdar and Barna Saha.  Clustering with noisy queries, 2017.
  * Mifflin [1977] Robert Mifflin.  An algorithm for constrained optimization with semismooth functions.  _Mathematics of Operations Research_ , 1977.
  * Mirsky [1960] Leon Mirsky.  Symmetric gauge functions and unitarily invariant norms.  _The quarterly journal of mathematics_ , 1960.
  * Munkres [1957] James Munkres.  Algorithms for the assignment and transportation problems.  _Journal of the society for industrial and applied mathematics_ , 1957\.
  * Nemirovski [2004] Arkadi Nemirovski.  Prox-method with rate of convergence o (1/t) for variational inequalities with lipschitz continuous monotone operators and smooth convex-concave saddle point problems.  _SIAM Journal on Optimization_ , 2004.
  * Nesterov [1989] Ju E Nesterov.  Self-concordant functions and polynomial-time methods in convex programming.  _Report, Central Economic and Mathematic Institute, USSR Acad. Sci_ , 1989.
  * Nesterov and Nemirovskii [1994] Yurii Nesterov and Arkadii Nemirovskii.  _Interior-point polynomial algorithms in convex programming_.  SIAM, 1994.
  * O’Donoghue et al. [2021] Brendan O’Donoghue, Eric Chu, Neal Parikh, and Stephen Boyd.  SCS: Splitting Conic Solver, version 3.1.1.  https://github.com/cvxgrp/scs, November 2021.
  * Overton [1992] Michael L Overton.  Large-scale optimization of eigenvalues.  _SIAM Journal on Optimization_ , 1992.
  * Raghunathan et al. [2018] Aditi Raghunathan, Jacob Steinhardt, and Percy S Liang.  Semidefinite relaxations for certifying robustness to adversarial examples.  _Advances in neural information processing systems_ , 2018.
  * Reinelt [1995] Gerhard Reinelt.  Tsplib95.  _Interdisziplinäres Zentrum für Wissenschaftliches Rechnen (IWR), Heidelberg_ , 1995.  URL http://comopt.ifi.uni-heidelberg.de/software/TSPLIB95/.
  * Rosen et al. [2019] David M Rosen, Luca Carlone, Afonso S Bandeira, and John J Leonard.  Se-sync: A certifiably correct algorithm for synchronization over the special euclidean group.  _The International Journal of Robotics Research_ , 2019.
  * Sahin et al. [2019] Mehmet Fatih Sahin, Ahmet Alacaoglu, Fabian Latorre, Volkan Cevher, et al.  An inexact augmented lagrangian framework for nonconvex optimization with nonlinear constraints.  _Advances in Neural Information Processing Systems_ , 2019.
  * Sato and Iwayama [2009] Yusuke Sato and Makoto Iwayama.  Interactive constrained clustering for patent document set.  In _Proceedings of the 2nd international workshop on Patent information retrieval_ , pages 17–20, 2009.
  * Schacke [2004] Kathrin Schacke.  On the kronecker product.  _Master’s thesis, University of Waterloo_ , 2004.
  * Shental et al. [2004] Noam Shental, Aharon Bar-Hillel, Tomer Hertz, and Daphna Weinshall.  Computing gaussian mixture models with em using equivalence constraints.  _Advances in neural information processing systems_ , 2004.
  * Souto et al. [2022] Mario Souto, Joaquim D Garcia, and Álvaro Veiga.  Exploiting low-rank structure in semidefinite programming by approximate operator splitting.  _Optimization_ , 2022.
  * Spielman and Srivastava [2008] Daniel A Spielman and Nikhil Srivastava.  Graph sparsification by effective resistances.  In _Proceedings of the fortieth annual ACM symposium on Theory of computing_ , 2008.
  * Sremac et al. [2021] Stefan Sremac, Hugo J Woerdeman, and Henry Wolkowicz.  Error bounds and singularity degree in semidefinite programming.  _SIAM Journal on Optimization_ , 2021.
  * [81] Jos F Sturm.  Using sedumi 1.02, a matlab toolbox for optimization over symmetric cones.  _Optimization methods and software_.
  * Sturm [2000] Jos F Sturm.  Error bounds for linear matrix inequalities.  _SIAM Journal on Optimization_ , 2000.
  * Subramanian et al. [2021] Shivashankar Subramanian, Daniel King, Doug Downey, and Sergey Feldman.  S2AND: A Benchmark and Evaluation System for Author Name Disambiguation.  In _Proceedings of the ACM/IEEE Joint Conference on Digital Libraries in 2021_ , 2021.
  * Todd [2001] Michael J Todd.  Semidefinite optimization.  _Acta Numerica_ , 2001.
  * Todd et al. [1998] Michael J Todd, Kim-Chuan Toh, and Reha H Tütüncü.  On the nesterov–todd direction in semidefinite programming.  _SIAM Journal on Optimization_ , 1998.
  * Toh et al. [1999] Kim-Chuan Toh, Michael J Todd, and Reha H Tütüncü.  Sdpt3—a matlab software package for semidefinite programming, version 1.3.  _Optimization methods and software_ , 1999.
  * Tropp et al. [2017] Joel A Tropp, Alp Yurtsever, Madeleine Udell, and Volkan Cevher.  Fixed-rank approximation of a positive-semidefinite matrix from streaming data.  _Advances in Neural Information Processing Systems_ , 2017.
  * Tsuda et al. [2005] Koji Tsuda, Gunnar Rätsch, and Manfred K Warmuth.  Matrix exponentiated gradient updates for on-line learning and bregman projection.  _Journal of Machine Learning Research_ , 2005.
  * Uhlen et al. [2010] Mathias Uhlen, Per Oksvold, Linn Fagerberg, Emma Lundberg, Kalle Jonasson, Mattias Forsberg, Martin Zwahlen, Caroline Kampf, Kenneth Wester, Sophia Hober, et al.  Towards a knowledge-based human protein atlas.  _Nature biotechnology_ , 2010.
  * Vandenberghe and Boyd [1999] Lieven Vandenberghe and Stephen Boyd.  Applications of semidefinite programming.  _Applied Numerical Mathematics_ , 1999.
  * Vandenberghe et al. [1998a] Lieven Vandenberghe, Stephen Boyd, and Abbas El Gamal.  Optimizing dominant time constant in RC circuits.  _IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems_ , 1998a.
  * Vandenberghe et al. [1998b] Lieven Vandenberghe, Stephen Boyd, and Shao-Po Wu.  Determinant maximization with linear matrix inequality constraints.  _SIAM journal on matrix analysis and applications_ , 1998b.
  * Vikram and Dasgupta [2016] Sharad Vikram and Sanjoy Dasgupta.  Interactive bayesian hierarchical clustering.  In _International Conference on Machine Learning_ , 2016.
  * Viswanathan et al. [2023] Vijay Viswanathan, Kiril Gashteovski, Carolin Lawrence, Tongshuang Wu, and Graham Neubig.  Large language models enable few-shot clustering.  _arXiv preprint arXiv:2307.00524_ , 2023.
  * Wagstaff and Cardie [2000] Kiri Wagstaff and Claire Cardie.  Clustering with instance-level constraints.  _AAAI/IAAI_ , 2000.
  * Wagstaff et al. [2001] Kiri Wagstaff, Claire Cardie, Seth Rogers, Stefan Schrödl, et al.  Constrained k-means clustering with background knowledge.  In _Icml_ , 2001.
  * Waldspurger and Waters [2020] Irene Waldspurger and Alden Waters.  Rank optimality for the burer–monteiro factorization.  _SIAM journal on Optimization_ , 2020.
  * Wang et al. [2017] Po-Wei Wang, Wei-Cheng Chang, and J Zico Kolter.  The mixing method: low-rank coordinate descent for semidefinite programming with diagonal constraints.  _arXiv preprint arXiv:1706.00476_ , 2017.
  * Wu et al. [1999] Kesheng Wu, Andrew Canning, HD Simon, and L-W Wang.  Thick-restart lanczos method for electronic structure calculations.  _Journal of Computational Physics_ , 1999.
  * Xing et al. [2002] Eric Xing, Michael Jordan, Stuart J Russell, and Andrew Ng.  Distance metric learning with application to clustering with side-information.  _Advances in neural information processing systems_ , 2002.
  * Xiong et al. [2016] Caiming Xiong, David M Johnson, and Jason J Corso.  Active clustering with model-based uncertainty reduction.  _IEEE transactions on pattern analysis and machine intelligence_ , 2016.
  * Yadav et al. [2021] Nishant Yadav, Nicholas Monath, Rico Angell, and Andrew McCallum.  Event and entity coreference using trees to encode uncertainty in joint decisions.  In _Proceedings of the Fourth Workshop on Computational Models of Reference, Anaphora and Coreference_ , 2021.
  * Yang et al. [2015] Liuqin Yang, Defeng Sun, and Kim-Chuan Toh.  Sdpnal+: a majorized semismooth newton-cg augmented lagrangian method for semidefinite programming with nonnegative constraints.  _Mathematical Programming Computation_ , 2015.
  * Yurtsever et al. [2015] Alp Yurtsever, Quoc Tran Dinh, and Volkan Cevher.  A universal primal-dual convex optimization framework.  _Advances in Neural Information Processing Systems_ , 2015.
  * Yurtsever et al. [2017] Alp Yurtsever, Madeleine Udell, Joel Tropp, and Volkan Cevher.  Sketchy decisions: Convex low-rank matrix optimization with optimal storage.  In _Artificial intelligence and statistics_. PMLR, 2017.
  * Yurtsever et al. [2019] Alp Yurtsever, Olivier Fercoq, and Volkan Cevher.  A conditional-gradient-based augmented lagrangian framework.  In _International Conference on Machine Learning_ , 2019.
  * Yurtsever et al. [2021] Alp Yurtsever, Joel A Tropp, Olivier Fercoq, Madeleine Udell, and Volkan Cevher.  Scalable semidefinite programming.  _SIAM Journal on Mathematics of Data Science_ , 2021.
  * Zaslavskiy et al. [2008] Mikhail Zaslavskiy, Francis Bach, and Jean-Philippe Vert.  A path following algorithm for the graph matching problem.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 2008.
  * Zhao et al. [1998] Qing Zhao, Stefan E Karisch, Franz Rendl, and Henry Wolkowicz.  Semidefinite programming relaxations for the quadratic assignment problem.  _Journal of Combinatorial Optimization_ , 1998.

##  Appendix A Solving Iteration Subproblems

In this section, we detail the primal-dual path-following interior point
methods used to solve the small semidefinite program to compute the value
$\hat{f}_{t}(\tilde{y}_{t+1})$ and the small quadratic semidefinite program
(24). Before we can derive the algorithms, it is necessary to define the
$\mathrm{svec}$ operator and symmetric Kronecker product [5, 76, 85].

For any matrix $A\in\SS^{n}$, the vector
$\mathrm{svec}(A)\in\mathbb{R}\mathclose{\vphantom{\big{)}}}^{n+1\choose 2}$
is defined as

$\mathrm{svec}(A)=\left[a_{11},\sqrt{2}a_{21},\ldots,\sqrt{2}a_{n1},a_{22},\sqrt{2}a_{32},\ldots,\sqrt{2}a_{n2},\ldots,a_{nn}\right]^{\top}.$

The $\mathrm{svec}$-operator is a structure preserving map between $\SS^{n}$
and $\mathbb{R}\mathclose{\vphantom{\big{)}}}^{n+1\choose 2}$ where the
constant $\sqrt{2}$ multiplied by some of the entries ensures that

$\langle
A,B\rangle=\mathrm{tr}(AB)=\mathrm{svec}(A)^{\top}\mathrm{svec}(B),\quad\forall\,A,B\in\SS^{n}.$

For any $M\in\mathbb{R}\mathclose{\vphantom{)}}^{n\times n}$, let
$\textrm{vec}(M)$ be the map from $\mathbb{R}\mathclose{\vphantom{)}}^{n\times
n}$ to $\mathbb{R}\mathclose{\vphantom{)}}^{n^{2}}$ defined by stacking the
columns of $M$ into a single $n^{2}$-dimensional vector. It is also useful to
define the matrix $U\in\mathbb{R}\mathclose{\vphantom{\big{)}}}^{{n+1\choose
2}\times n^{2}}$ which maps $\textrm{vec}(A)\mapsto\mathrm{svec}(A)$ for any
$A\in\SS^{n}$. Let $u_{ij,kl}$ be the entry in the row which defines element
$a_{ij}$ in $\mathrm{svec}(A)$ and the column that is multiplied with the
element $a_{kl}$ in $\textrm{vec}(A)$. Then

$u_{ij,kl}=\begin{cases}1&\quad i=j=k=l\\\ \frac{1}{\sqrt{2}}&\quad i=k\neq
j=l\textrm{, or }i=l\neq j=k\\\ 0&\quad\textrm{o.w.}\end{cases}$

As an example, the case when $n=2$:

$U=\begin{bmatrix}1&0&0&0\\\ 0&\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}}&0\\\
0&0&0&1\\\ \end{bmatrix}.$

Note that $U$ is a unique matrix with orthonormal rows and has the following
property

$U^{\top}U\textrm{vec}(A)=U^{\top}\mathrm{svec}(A)=\textrm{vec}(A),\quad\forall\,A\in\SS^{n}.$

The symmetric Kronecker product $\otimes_{s}$ can be defined for any two
square matrices $G,H\in\mathbb{R}^{n\times n}$ by its action on a vector
$\mathrm{svec}(A)$ for $A\in\SS^{n}$ as follows

$(G\otimes_{s}H)\,\mathrm{svec}(A)=\frac{1}{2}\mathrm{svec}(HAG^{\top}+GAH^{\top}).$

Alternatively, but equivalently [76], the symmetric Kronecker product can be
defined more explicitly using the matrix $U$ defined above as follows

$G\otimes_{s}H=\frac{1}{2}U(G\otimes H+H\otimes G)U^{\top},$

where $\otimes$ is the standard Kronecker product. We use this latter
definition in our implementation.

###  A.1 Computing $\hat{f}_{t}(\tilde{y}_{t+1})$

The value $\hat{f}_{t}(\tilde{y}_{t+1})$ is the optimum value of the following
optimization problem (remember that $\nu$ can be dropped since
$\tilde{y}_{t+1}$ is always feasible)

$\displaystyle\operatorname*{maximize}\quad\langle
C-\mathcal{A}\mathclose{\vphantom{)}}^{*}\tilde{y}_{t+1},\eta\bar{X}_{t}+V_{t}SV_{t}^{\top}\rangle+\langle
b,\tilde{y}_{t+1}\rangle$
$\displaystyle\operatorname*{subject\;to}\;\;\;\eta\geq 0$
$\displaystyle\qquad\qquad\quad\;\;S\succeq 0$
$\displaystyle\qquad\qquad\quad\;\;\eta+\mathrm{tr}(S)\leq\alpha$

This amounts to computing $\langle
C-\mathcal{A}\mathclose{\vphantom{)}}^{*}\tilde{y}_{t+1},\eta_{\star}\bar{X}_{t}+V_{t}S_{\star}V_{t}^{\top}\rangle+\langle
b,\tilde{y}_{t+1}\rangle$ where $(\eta_{\star},S_{\star})$ is a solution to
the following (small) semidefinite program

$\displaystyle\operatorname*{minimize}\quad
g_{1}^{\top}\mathrm{svec}(S)+\eta\,g_{2}$ (52)
$\displaystyle\operatorname*{subject\;to}\;\;\;\eta\geq 0$
$\displaystyle\qquad\qquad\quad\;\;S\succeq 0$
$\displaystyle\qquad\qquad\quad\;\;1-v_{I}^{\top}\mathrm{svec}(S)-\eta\geq 0$

where

$\displaystyle g_{1}$
$\displaystyle=\alpha\,\mathrm{svec}(V_{t}^{\top}\\!(\mathcal{A}\mathclose{\vphantom{)}}^{*}\tilde{y}_{t+1}-C)V_{t}),$
$\displaystyle g_{2}$
$\displaystyle=\frac{\alpha}{\mathrm{tr}\\!\left(\bar{X}\right)}\langle\bar{X}_{t},\mathcal{A}\mathclose{\vphantom{)}}^{*}\tilde{y}_{t+1}-C)\rangle,$
$\displaystyle v_{I}$ $\displaystyle=\mathrm{svec}(I).\qquad(I\textrm{ is the
}k\times k\textrm{ identity matrix})$

We will use a primal-dual interior point method to solve (52). We follow the
well known technique for deriving primal-dual interior point methods [40]. We
start by defining the Lagrangian of the dual barrier problem of (52),

$\displaystyle L_{\mu}(S,\eta,T,\zeta,\omega)$
$\displaystyle=g_{1}^{\top}\mathrm{svec}(S)+\eta
g_{2}-\mathrm{svec}(S)^{\top}\mathrm{svec}(T)-\eta\zeta$ (53)
$\displaystyle\qquad-\omega(1-v_{I}^{\top}\mathrm{svec}(S)-\eta)+\mu(\log\det(T)+\log\zeta+\log\omega),$

where we introduce a dual slack matrix $T\succeq 0$ as complementary to $S$, a
dual slack scalar $\zeta\geq 0$ as complementary to $\eta$, a Lagrange
multiplier $\omega\geq 0$ for the trace constraint inequality, and a barrier
parameter $\mu>0$. Notice that we have moved from needing to optimize a
constrained optimization problem to an unconstrained optimization problem. The
saddle point solution of (53) is given by the solution of the KKT-conditions,
which reduces to just the first-order optimality conditions since the problem
is unconstrained. The first-order optimality conditions of (53) are the
following system of equations

$\displaystyle\nabla_{S}L_{\mu}$ $\displaystyle=g_{1}-\mathrm{svec}(T)+\omega
v_{I}=0$ (54) $\displaystyle\nabla_{\eta}L_{\mu}$
$\displaystyle=g_{2}-\zeta+\omega=0$ (55) $\displaystyle\nabla_{T}L_{\mu}$
$\displaystyle=S-\mu T^{-1}=0$ (56) $\displaystyle\nabla_{\zeta}L_{\mu}$
$\displaystyle=\eta-\mu\zeta^{-1}=0$ (57)
$\displaystyle\nabla_{\omega}L_{\mu}$
$\displaystyle=1-v_{I}^{\top}\mathrm{svec}(S)-\eta-\mu\omega^{-1}=0$ (58)

By the strict concavity of $\log\det T$, $\log\zeta$, and $\log\omega$, there
exists a unique solution
$(S_{\mu},\eta_{\mu},T_{\mu},\zeta_{\mu},\omega_{\mu})$ to this system of
equations for any value of the barrier parameter $\mu>0$. The sequence of
these solutions as $\mu\to 0$ forms the _central trajectory_ (also known as
the central _path_). For a point $(S,\eta,T,\zeta,\omega)$ on the central
trajectory, we can use any combination of (56), (57), and/or (58) to solve for
$\mu$,

$\mu=\frac{\langle
S,T\rangle}{k}=\eta\zeta=\omega(1-v_{I}^{\top}\mathrm{svec}(S)-\eta)=\frac{\langle
S,T\rangle+\eta\zeta+\omega(1-v_{I}^{\top}\mathrm{svec}(S)-\eta)}{k+2}.$ (59)

The fundamental idea of primal-dual interior point methods is to use Newton’s
method to follow the central path to a solution of (52). Before we can apply
Newton’s method, we must linearize the non-linear equations (56), (57), and
(58) into an equivalent linear formulation. There are several ways one could
linearize these equations and the choice of linearization significantly
impacts the algorithm’s behavior (see [5] or [40] for more on the choice of
linearization). We choose one standard method linearization, detailed in the
following system of equations

$\mathrm{F}_{\mu}(\theta)=\mathrm{F}_{\mu}(S,\eta,T,\zeta,\omega):=\begin{pmatrix}\vspace{0.1em}g_{1}-\mathrm{svec}(T)+\omega
v_{I}\\\ \vspace{0.1em}g_{2}-\zeta+\omega\\\ \vspace{0.1em}ST-\mu I\\\
\vspace{0.1em}\eta\zeta-\mu\\\
\vspace{0.1em}\omega(1-v_{I}^{\top}\mathrm{svec}(S)-\eta)-\mu\\\
\end{pmatrix}=0.$ (60)

The solution $\theta_{\star}$ to this system of equations
$\mathrm{F}_{\mu}(\theta)=0$ satisfies the first-order optimality conditions
(54)-(58) and is the optimal solution to the barrier problem. We utilize
Newton’s method to take steps in the update direction $\Delta\theta=(\Delta
S,\Delta\eta,\Delta T,\Delta\zeta,\Delta\omega)$ towards $\theta_{\star}$. The
update direction $\Delta\theta$ determined by Newton’s method must satisfy the
following equation

$\mathrm{F}_{\mu}(\theta)+\nabla\mathrm{F}_{\mu}(\Delta\theta)=0.$

Hence, the update direction $\Delta\theta$ is the solution to the following
system of equations (after the same standard linearization has been applied to
the following system as in (60))

$\displaystyle-\mathrm{svec}(\Delta T)+\Delta\omega\,v_{I}$
$\displaystyle=\mathrm{svec}(T)-g_{1}-\omega\,v_{I}$ (61)
$\displaystyle-\Delta\zeta+\Delta\omega$ $\displaystyle=\zeta-g_{2}-\omega$
(62)
$\displaystyle\omega^{-1}(1-v_{I}^{\top}\mathrm{svec}(S)-\eta)\Delta\omega-
v_{I}^{\top}\mathrm{svec}(\Delta S)-\Delta\eta$
$\displaystyle=\mu\omega^{-1}+v_{I}^{\top}\mathrm{svec}(S)+\eta-1$ (63)
$\displaystyle(T\otimes_{s}S^{-1})\,\mathrm{svec}(\Delta
S)+\mathrm{svec}(\Delta T)$
$\displaystyle=\mu\,\mathrm{svec}(S^{-1})-\mathrm{svec}(T)$ (64)
$\displaystyle\zeta\eta^{-1}\Delta\eta+\Delta\zeta$
$\displaystyle=\mu\eta^{-1}-\zeta$ (65)

We can solve this system efficiently by analytically eliminating all variables
except $\mathrm{svec}(\Delta S)$. We can compute $\mathrm{svec}(\Delta S)$ by
solving the following linear matrix equation using an off-the-shelf linear
system solver

$\displaystyle\left(T\otimes_{s}S^{-1}+\frac{\zeta\eta^{-1}}{\kappa_{1}\zeta\eta^{-1}+1}v_{I}v_{I}^{\top}\right)\mathrm{svec}(\Delta
S)=v_{I}g_{2}-v_{I}\mu\eta^{-1}-g_{1}+\mu\,\mathrm{svec}(S^{-1})$ (66)
$\displaystyle\qquad\qquad+v_{I}\zeta\eta^{-1}(-\kappa_{1}\zeta\eta^{-1}-1)^{-1}\left(-\kappa_{1}(\mu\eta^{-1}-g_{2}-\omega)+\mu\omega^{-1}+v_{I}^{\top}\mathrm{svec}(S)+\eta-1\right)$

where $\kappa_{1}:=\omega^{-1}(1-v_{I}^{\top}\mathrm{svec}(S)-\eta)$. Then,
computing the rest of the update directions $\Delta\eta$,
$\mathrm{svec}(\Delta T)$, $\Delta\zeta$, and $\Delta\omega$ amounts to back-
substituting the solution to (66) for $\mathrm{svec}(\Delta S)$ through the
analytical variable elimination equations.

Given how to compute the update directions, the primal-dual interior point
method proceeds as follows. Initialize the variables
$\theta=(S,\eta,T,\zeta,\omega)$ to an arbitrary strictly feasible point (i.e.
$S\succ 0$, $\eta>0$, $T\succ 0$, $\zeta>0$, and $\omega>0$). Starting from
this primal-dual pair we compute an estimate of the barrier parameter as
follows

$\mu\leftarrow\frac{\langle
S,T\rangle+\eta\zeta+\omega(1-v_{I}^{\top}\mathrm{svec}(S)-\eta)}{2(k+2)},$

where, as done in [40], we use (59) and divide by two. Then, we compute the
update direction $\Delta\theta$ as described above and perform a backtracking
line search to find a step size $\delta\in(0,1]$ such that
$\theta+\delta\,\Delta\theta$ is again strictly feasible. Lastly, following
[42], we compute a non-increasing estimate of the barrier parameter

$\mu\leftarrow\min\left\\{\mu_{\textrm{prev}},\gamma\frac{\langle
S,T\rangle+\eta\zeta+\omega(1-v_{I}^{\top}\mathrm{svec}(S)-\eta)}{2(k+2)}\right\\}\quad\textrm{where}\quad\gamma=\begin{cases}1&\textrm{if
}\delta\leq\frac{1}{5}\\\ \frac{5}{10}-\frac{4}{10}\delta^{2}&\textrm{if
}\delta>\frac{1}{5}\end{cases}.$

We iterate over these steps until the barrier parameter $\mu$ is small enough
(e.g. $\mu<10^{-7}$).

###  A.2 Solving (24)

The subproblem (24) can be rewritten as the following (small) quadratic
semidefinite program

$\displaystyle\operatorname*{maximize}\quad\langle
C-\mathcal{A}\mathclose{\vphantom{)}}^{*}y_{t},\eta\bar{X}_{t}+V_{t}SV_{t}^{\top}\rangle+\langle
b-\tilde{\nu},y_{t}\rangle-\frac{1}{2\rho}\left\|b-\tilde{\nu}-\mathcal{A}(\eta\bar{X}_{t}+V_{t}SV_{t}^{\top})\right\|_{2}^{2}$
$\displaystyle\operatorname*{subject\;to}\;\;\;\eta\geq 0$
$\displaystyle\qquad\qquad\quad\;\;S\succeq 0$
$\displaystyle\qquad\qquad\quad\;\;\eta+\mathrm{tr}(S)\leq\alpha$

For this subproblem, unlike the subproblem solved in subsection A.1, we are
solving for $\eta$ and $S$ to compute the candidate iterate $\tilde{y}_{t+1}$,
update the primal variable, and update the model. This (small) quadratic
semidefinite program is equivalent to the following optimization problem

$\displaystyle\operatorname*{minimize}\quad\frac{1}{2}\,\mathrm{svec}(S)^{\top}Q_{11}\,\mathrm{svec}(S)+\eta\,q_{12}^{\top}\,\mathrm{svec}(S)+\frac{1}{2}\eta^{2}q_{22}+h_{1}^{\top}\mathrm{svec}(S)+\eta\,h_{2}$
(67) $\displaystyle\operatorname*{subject\;to}\;\;\;\eta\geq 0$
$\displaystyle\qquad\qquad\quad\;\;S\succeq 0$
$\displaystyle\qquad\qquad\quad\;\;1-v_{I}^{\top}\mathrm{svec}(S)-\eta\geq 0$

where

$\displaystyle Q_{11}$
$\displaystyle=\frac{\alpha^{2}}{\rho}\sum_{i=1}^{m}\mathrm{svec}(V_{t}^{\top}\\!A_{i}V_{t})\,\mathrm{svec}(V_{t}^{\top}\\!A_{i}V_{t})^{\top}$
$\displaystyle q_{12}$
$\displaystyle=\frac{\alpha^{2}}{\rho\,\mathrm{tr}\\!\left(\bar{X}_{t}\right)}\,\mathrm{svec}\left(V_{t}^{\top}\\!\mathcal{A}^{*}\mathcal{A}\,\bar{X}_{t}V_{t}\right)$
$\displaystyle q_{22}$
$\displaystyle=\frac{\alpha^{2}}{\rho\,\mathrm{tr}\\!\left(\bar{X}_{t}\right)^{2}}\left\langle\mathcal{A}\bar{X}_{t},\mathcal{A}\bar{X}_{t}\right\rangle$
$\displaystyle h_{1}$
$\displaystyle=\alpha\,\,\mathrm{svec}\left(V_{t}^{\top}\\!\left(\mathcal{A}\mathclose{\vphantom{)}}^{*}y_{t}-C-\frac{1}{\rho}\mathcal{A}\mathclose{\vphantom{)}}^{*}(b-\tilde{\nu})\right)V_{t}\right)$
$\displaystyle h_{2}$
$\displaystyle=\frac{\alpha}{\mathrm{tr}\\!\left(\bar{X}_{t}\right)}\left\langle\bar{X},\,\mathcal{A}\mathclose{\vphantom{)}}^{*}y_{t}-C-\frac{1}{\rho}\mathcal{A}\mathclose{\vphantom{)}}^{*}(b-\tilde{\nu})\right\rangle$
$\displaystyle v_{I}$ $\displaystyle=\mathrm{svec}(I)$

We follow the same derivation procedure as in subsection A.1, so we proceed by
including the details which differ from the previous section. The Lagrangian
of the dual barrier problem of (67) is as follows

$\displaystyle L_{\mu}(S,\eta,T,\zeta,\omega)$
$\displaystyle=\frac{1}{2}\,\mathrm{svec}(S)^{\top}Q_{11}\,\mathrm{svec}(S)+\eta\,q_{12}^{\top}\,\mathrm{svec}(S)$
(68)
$\displaystyle\qquad+\frac{1}{2}\eta^{2}q_{22}+h_{1}^{\top}\mathrm{svec}(S)+\eta\,h_{2}-\mathrm{svec}(S)^{\top}\mathrm{svec}(T)-\eta\zeta$
$\displaystyle\qquad-\omega(1-v_{I}^{\top}\mathrm{svec}(S)-\eta)+\mu(\log\det(T)+\log\zeta+\log\omega).$

The first-order optimality conditions of (68) yields the following system of
equations (after the same standard linearization)

$\mathrm{F}_{\mu}(S,\eta,T,\zeta,\omega):=\begin{pmatrix}\vspace{0.1em}Q_{11}\,\mathrm{svec}(S)+\eta\,q_{12}+h_{1}-\mathrm{svec}(T)+\omega
v_{I}\\\
\vspace{0.1em}q_{12}^{\top}\,\mathrm{svec}(S)+\eta\,q_{22}+h_{2}-\zeta+\omega\\\
\vspace{0.1em}ST-\mu I\\\ \vspace{0.1em}\eta\zeta-\mu\\\
\vspace{0.1em}\omega(1-v_{I}^{\top}\mathrm{svec}(S)-\eta)-\mu\\\
\end{pmatrix}=:\begin{pmatrix}\vspace{0.1em}\mathsf{F}_{1}\\\
\vspace{0.1em}\mathsf{F}_{2}\\\ \vspace{0.1em}\mathsf{F}_{3}\\\
\vspace{0.1em}\mathsf{F}_{4}\\\ \vspace{0.1em}\mathsf{F}_{5}\\\
\end{pmatrix}=0.$ (69)

The Newton’s method step direction $(\Delta S,\Delta\eta,\Delta
T,\Delta\zeta,\Delta\omega)$ is determined via the following linearized system

$\displaystyle Q_{11}\mathrm{svec}(\Delta
S)+\Delta\eta\,q_{12}-\mathrm{svec}(\Delta T)+\Delta\omega\,v_{I}$
$\displaystyle=-\mathsf{F}_{1}$ (70) $\displaystyle
q_{12}^{\top}\,\mathrm{svec}(\Delta S)+\Delta\,\eta
q_{22}-\Delta\zeta+\Delta\omega$ $\displaystyle=-\mathsf{F}_{2}$ (71)
$\displaystyle\omega^{-1}(1-v_{I}^{\top}\mathrm{svec}(S)-\eta)\Delta\omega-
v_{I}^{\top}\mathrm{svec}(\Delta S)-\Delta\eta$
$\displaystyle=\mu\omega^{-1}+v_{I}^{\top}\mathrm{svec}(S)+\eta-1$ (72)
$\displaystyle(T\otimes_{s}S^{-1})\,\mathrm{svec}(\Delta
S)+\mathrm{svec}(\Delta T)$
$\displaystyle=\mu\,\mathrm{svec}(S^{-1})-\mathrm{svec}(T)$ (73)
$\displaystyle\zeta\eta^{-1}\Delta\eta+\Delta\zeta$
$\displaystyle=\mu\eta^{-1}-\zeta$ (74)

We can solve this system efficiently by analytically eliminating all variables
except $\mathrm{svec}(\Delta S)$. We can compute $\mathrm{svec}(\Delta S)$ by
solving the following linear matrix equation using an off-the-shelf linear
system solver

$\displaystyle\left(Q_{11}+T\otimes_{s}S^{-1}-(\kappa_{1}\kappa_{2}+1)^{-1}\big{(}q_{12}(\kappa_{1}q_{12}+v_{I})^{\top}\\!+v_{I}(q_{12}-\kappa_{2}v_{I})^{\top}\big{)}\right)\mathrm{svec}(\Delta
S)$ (75)
$\displaystyle\qquad=q_{12}\left((\kappa_{1}\kappa_{2}+1)^{-1}\left(\mu\omega^{-1}+v_{I}^{\top}\mathrm{svec}(S)+\eta-1+\kappa_{1}(\textsf{F}_{2}-\mu\eta^{-1}+\zeta)\right)\right)$
$\displaystyle\qquad\qquad+v_{I}\left((\kappa_{1}\kappa_{2}+1)^{-1}\left(\textsf{F}_{2}-\mu\eta^{-1}+\zeta-\kappa_{2}(\mu\omega^{-1}+v_{I}^{\top}\mathrm{svec}(S)+\eta-1)\right)\right)$
$\displaystyle\qquad\qquad-\textsf{F}_{1}+\mu\,\mathrm{svec}(S^{-1})-\mathrm{svec}(T),$

where $\kappa_{1}:=\omega^{-1}(1-v_{I}^{\top}\mathrm{svec}(S)-\eta)$ and
$\kappa_{2}:=\zeta\eta^{-1}+q_{22}$. Then, computing the rest of the update
directions $\Delta\eta$, $\mathrm{svec}(\Delta T)$, $\Delta\zeta$, and
$\Delta\omega$ amounts to back-substituting the solution to (75) for
$\mathrm{svec}(\Delta S)$ through the analytical variable elimination
equations.

We make a single change to the initialization as compared with the interior
point method presented in subsection A.1. Since this interior point method can
be executed multiple times in a row with only the value of $\tilde{\nu}$
changing as the first step in the alternating maximization algorithm, we warm-
start initialize $\theta$ with the previous execution’s $\theta_{\star}$. We
observe non-negligible convergence improvements over arbitrary initialization
as this interior point method procedure is called in sequence.

##  Appendix B Proofs of Theoretical Results

###  B.1 Proof of Fact 3.2

The function $f$ is continuous and proper over the domain $\mathcal{Y}$. It is
a well-known fact that the level sets of continuous proper functions are
compact. Hence, the union of level sets $\\{y\in\mathcal{Y}:f(y_{\star})\leq
f(y)\leq f(y_{0})\\}$ is compact, and therefore, $\sup_{\,t\geq
0}\\{\mathrm{dist}(y_{t},\mathcal{Y}_{*})\\}<\infty$. The subgradients take
the form $g_{t+1}=b+\alpha\mathcal{A}(vv^{\\!\top})\in\partial
f(\tilde{y}_{t+1})$ where $v\in\mathbb{R}^{n}$ is a unit-normed vector, and
thus, $\|g_{t+1}\|<\infty$ for all $t$.

###  B.2 Proof of Lemma 3.3

Let
$\mathcal{F}:=\\{y\in\mathbb{R}^{m}:\mathcal{A}\mathclose{\vphantom{)}}^{*}y-C\in\SS^{n}_{+}\\}$.
Observe that $\mathcal{Y}_{\star}=\mathcal{Y}\cap\mathcal{F}$. Since
$\mathcal{Y}_{\star}$ is compact, the Hölderian error bound [82, 31, 80]
ensures that

$\mathrm{dist}(y,\mathcal{Y}_{\star})\leq\mathcal{O}\left(\mathrm{dist}^{2^{-d}}(y,\mathcal{Y})+\mathrm{dist}^{2^{-d}}(y,\mathcal{F})\right),$

where $d\in\\{0,1,\ldots,m\\}$ is the singularity degree of the SDP. This
implies that for $y\in Y$,

$\displaystyle\mathrm{dist}^{2^{d}}(y,\mathcal{Y}_{\star})$
$\displaystyle\leq\mathcal{O}\left(\mathrm{dist}(y,\mathcal{F})\right)$
$\displaystyle\leq\mathcal{O}\left(\alpha\max\\{\lambda_{\textrm{max}}(C-\mathcal{A}\mathclose{\vphantom{)}}^{*}y),0\\}\right)$
$\displaystyle\leq\mathcal{O}\left(f(y)-f(y_{\star})\right).$

###  B.3 Proof of Lemma 3.4

For any descent step $t$, it is easy to verify from 22 and (25) that

$\mathrm{proj}_{\mathcal{K}}(\mathcal{A}X_{t+1})-\mathcal{A}X_{t+1}\leq\rho(y_{t}-y_{t+1}),$

and therefore,

$\|\mathcal{A}X_{t+1}-\mathrm{proj}_{\mathcal{K}}(\mathcal{A}X_{t+1})\|^{2}\leq\rho^{2}\|y_{t}-y_{t+1}\|^{2}.$

To complete the proof we utilize the fact that
$\hat{f}_{t}(y_{t+1})=\min_{y\in\mathcal{Y}}\hat{f}_{t}(y)+\frac{\rho}{2}\|y-y_{t}\|^{2}$,
the fact that $\hat{f}_{t}(y)\leq f(y)$ for all $y\in\mathcal{Y}$, and the
definition of a descent step, which gives

$\frac{\rho}{2}\|y_{t+1}-y_{t}\|^{2}\leq\hat{f}_{t}(y_{t})-\hat{f}_{t}(y_{t+1})\leq
f(y_{t})-\hat{f}_{t}(y_{t+1})\leq\frac{f(y_{t})-f(y_{t+1})}{\beta}\leq\frac{f(y_{t})-f(y_{\star})}{\beta}.$

Assume Slater’s condition holds and $y_{\star}$ is unique. Then, using Lemma
3.3 gives

$\displaystyle\|\mathcal{A}X_{t+1}-\mathrm{proj}_{\mathcal{K}}(\mathcal{A}X_{t+1})\|^{2}$
$\displaystyle\leq\rho^{2}\|y_{t}-y_{t+1}\|^{2}$
$\displaystyle\leq\rho^{2}\left(\|y_{t}-y_{\star}\|^{2}+\|y_{t+1}-y_{\star}\|^{2}\right)$
$\displaystyle\leq\mathcal{O}\\!\left((f(y_{t})-f(y_{\star}))^{2}\right)+\mathcal{O}\\!\left((f(y_{t+1})-f(y_{\star}))^{2}\right)$
$\displaystyle\leq\mathcal{O}\\!\left((f(y_{t})-f(y_{\star}))^{2}\right).$

###  B.4 Proof of Lemma 3.5

By definition of strong duality, we know that for any
$X_{\star}\in\mathcal{X}_{\star}$ that $\langle C,X_{\star}\rangle=\langle
b,y_{\star}\rangle$, and equivalently, $\langle
X_{\star},\mathcal{A}\mathclose{\vphantom{)}}^{*}y_{\star}-C\rangle=0$. Then,
the dual objective gap is bounded as follows

$\displaystyle\langle b,y_{t}-y_{\star}\rangle$
$\displaystyle=\langle\mathcal{A}X_{\star},y_{t}-y_{\star}\rangle$
$\displaystyle=\langle
X_{\star},\mathcal{A}\mathclose{\vphantom{)}}^{*}(y_{t}-y_{\star})\rangle$
$\displaystyle=\langle
X_{\star},(\mathcal{A}\mathclose{\vphantom{)}}^{*}y_{t}-C)-(\mathcal{A}\mathclose{\vphantom{)}}^{*}y_{\star}-C)\rangle$
$\displaystyle=\langle
X_{\star},\mathcal{A}\mathclose{\vphantom{)}}^{*}y_{t}-C\rangle$
$\displaystyle\geq-\|X_{\star}\|\mathclose{\vphantom{)}}_{*}\max\\{\lambda_{\max}(C-\mathcal{A}\mathclose{\vphantom{)}}^{*}y_{t}),0\\}.$

We now utilize this bound to obtain the desired result

$f(y_{t})-f(y_{\star})=\langle
b,y_{t}-y_{\star}\rangle+\alpha\max\\{\lambda_{\max}(C-\mathcal{A}\mathclose{\vphantom{)}}^{*}y_{t}),0\\}\geq\|X_{\star}\|\mathclose{\vphantom{)}}_{*}\max\\{\lambda_{\max}(C-\mathcal{A}\mathclose{\vphantom{)}}^{*}y_{t}),0\\},$

where the last inequality is implied by the assumption that $\alpha\geq
2\mathbb{N}(\mathcal{X}_{\star})\geq
2\|X_{\star}\|\mathclose{\vphantom{)}}_{*}$.

###  B.5 Proof of Lemma 3.6

We start by rewriting the primal-dual gap as follows

$\displaystyle\langle C,X_{t+1}\rangle-\langle b,y_{t+1}\rangle$
$\displaystyle=\langle
C,X_{t+1}\rangle-\langle\mathcal{A}X_{t+1},y_{t+1}\rangle+\langle\mathcal{A}X_{t+1}-b,y_{t+1}\rangle$
$\displaystyle=\langle
X_{t+1},C-\mathcal{A}\mathclose{\vphantom{)}}^{*}y_{t+1}\rangle+\langle\mathcal{A}X_{t+1}-b,y_{t+1}\rangle.$

We will now bound the absolute values of the two resulting terms to bound the
desired quantity. Using the Cauchy-Schwarz inequality and Lemma 3.4 it can be
seen

$\displaystyle|\langle\mathcal{A}X_{t+1}-b,y_{t+1}\rangle|$
$\displaystyle\leq\|\mathcal{A}X_{t+1}-b\|\|y_{t+1}\|$
$\displaystyle\leq\mathcal{O}\\!\left(\|y_{t+1}\|\sqrt{f(y_{t})-f(y_{\star})}\right)$
$\displaystyle\leq\mathcal{O}\\!\left(\sqrt{f(y_{t})-f(y_{\star})}\right),$

where the last inequality comes from Fact 3.2. Assume Slater’s condition holds
and $y_{\star}$ is unique. Then, using Lemma 3.3 gives

$|\langle\mathcal{A}X_{t+1}-b,y_{t+1}\rangle|\leq\mathcal{O}\left(f(y_{t})-f(y_{\star})\right).$

Since $X_{t+1}\succeq 0$ and $\mathrm{tr}(X_{t+1})\leq\alpha$ by construction,
we can use Lemma 3.5 to yield

$\displaystyle|\langle
X_{t+1},C-\mathcal{A}\mathclose{\vphantom{)}}^{*}y_{t+1}\rangle|$
$\displaystyle\leq\mathcal{O}(\max\\{C-\mathcal{A}\mathclose{\vphantom{)}}^{*}y_{t+1},0\\})$
$\displaystyle\leq\mathcal{O}(f(y_{t+1})-f(y_{\star}))$
$\displaystyle\leq\mathcal{O}(f(y_{t})-f(y_{\star})).$

The immediately yields the desired result.

###  B.6 Proof of Theorem 3.1

The overall proof strategy is to use Theorem 2.1 to show that the penalized
dual gap, $f(y_{t})-f(y_{*})$, converges at a worst case rate of
$\mathcal{O}(1/\varepsilon^{3})$, which improves to
$\mathcal{O}(1/\varepsilon)$ if Slater’s condition or strict complementarity
hold. Then, we can use Lemmas 3.4, 3.5, and 3.6 to obtain the stated
convergence of primal feasibility, dual feasibility, and primal-dual
optimality.

To apply the result of Theorem 2.1, we showed that the norms of the iterates
and subgradients of $f$ are bounded (Fact 3.2), and so all we need to do is
show that the model satisfies the conditions (4), (5), and (6).

Let $v$ be a maximum eigenvector
$C-\mathcal{A}\mathclose{\vphantom{)}}^{*}\tilde{y}_{t+1}$ if
$\lambda_{\max}(C-\mathcal{A}\mathclose{\vphantom{)}}^{*}\tilde{y}_{t+1})>0$
and $v=0$, otherwise. Then denote
$g_{t+1}=b+\alpha\mathcal{A}(vv^{\top})\in\partial f(\tilde{y}_{t+1})$ as the
subgradient of $f$ corresponding $v$ at the candidate iterate
$\tilde{y}_{t+1}$. Let
$s_{t+1}=\rho(y_{t}-\tilde{y}_{t+1})\in\partial\hat{f}_{t}(\tilde{y}_{t+1})$
be the aggregate subgradient. Recall that at iteration $t+1$ the model
approximates $\mathcal{X}$ by the low-dimensional spectral set

$\widehat{\mathcal{X}}_{t+1}:=\left\\{\eta\bar{X}_{t+1}+V_{t+1}SV_{t+1}^{\top}:\eta\,\mathrm{tr}(\bar{X}_{t+1})+\mathrm{tr}(S)\leq\alpha,\eta\geq
0,S\in\SS^{k}_{+}\right\\}.$

and therefore the penalized dual objective and model respectively take the
following forms for any $y\in\mathcal{Y}$,

$f(y)=\sup_{X\in\mathcal{X}}\langle
C-\mathcal{A}\mathclose{\vphantom{)}}^{*}y,X\rangle+\langle
b,y\rangle\quad\textrm{and}\quad\hat{f}_{t+1}(y)=\sup_{X\in\widehat{\mathcal{X}}_{t+1}}\langle
C-\mathcal{A}\mathclose{\vphantom{)}}^{*}y,X\rangle+\langle b,y\rangle.$ (76)

##### Verifying (4).

The condition (4) follows immediately from (76), since
$\widehat{\mathcal{X}}_{t+1}\subseteq\mathcal{X}$.

##### Verifying (5).

Since $V_{t+1}$ spans $v$, there exists a vector $s\in\mathbb{R}^{k}$ such