# Momentum-inspired Low-Rank Coordinate Descent
for Diagonally Constrained SDPs

David S. Hippocampus
Department of Computer Science
Cranberry-Lemon University
Pittsburgh, PA 15213
<EMAIL_ADDRESS>
Use footnote for providing further information about author (webpage,
alternative address)— _not_ for acknowledging funding agencies.

###### Abstract

We present a novel, practical and provable approach for solving diagonally
constrained semi-definite programming (SDP) problems at scale _using
accelerated non-convex programming_. Our algorithm non-trivially combines
acceleration motions from convex optimization with coordinate power iteration
and matrix factorization techniques. The algorithm is extremely simple to
implement, and adds only a single extra hyperparameter – momentum – for which
we provide tuning bounds that hold in theory and practice. From a theoretical
perspective, we prove that our method admits local linear convergence in the
neighborhood of the optimum, and always converges to a first-order critical
point. We showcase the experimental merits of our method on three major
application domains: MaxCut, MaxSAT, and MIMO signal detection. In all cases,
our methodology provides significant speedups over non-convex and convex SDP
solvers – $5\times$ faster than state-of-the-art non-convex solvers and $9$ to
$10^{3}\times$ faster than convex SDP solvers – while maintaining comparable
solution quality.

##  1 Introduction

Background. This work focuses on efficient ways to solve semidefinite
programming instances (SDPs) with diagonal constraints:

$\displaystyle\min_{X\succeq 0}$ $\displaystyle\left\langle C,X\right\rangle$
s.t. $\displaystyle X_{ii}=1,~{}i=1,\dots,n.$ (1)

Here, $X\in\mathbb{S}^{n}$ is the symmetric optimization variable of dimension
$n\times n$, and $C\in\mathbb{S}^{n}$ is a problem-dependent cost matrix. The
above formulation usually appears in practice as the convex relaxation of
quadratic form optimization over discrete variables.

$\displaystyle\underset{x\in\mathcal{D}^{n}}{\min}$ $\displaystyle\langle
x,Cx\rangle,$ (2)

where $\mathcal{D}$ is a discrete set on the unit imaginary circle defined as
follows.

$\displaystyle\mathcal{D}=\left\\{e^{\frac{\imath 2\pi
m}{M}}~{}:~{}m=0,1,\dots,M-1\right\\}.$

Here, $\imath:=\sqrt{-1}$ denotes the imaginary unit. Different values for $M$
define different realizations of the discrete set $\mathcal{D}$. For example,
$M=2$ reduces $\mathcal{D}$ to the binary set $\\{\pm 1\\}$, while
$\mathcal{D}$ becomes $\\{\pm 1,\pm\imath\\}$ for $M=4$.111Many other
realizations of the set $\mathcal{D}$ exist. For example, in the case of
Max-$3$-Cut goemans2004approximation; so2007approximating, we can define
$\mathcal{D}=\\{1,e^{\frac{\imath 2\pi}{3}},e^{\frac{\imath 4\pi}{3}}\\}$. The
formulation in (2) appears in many applications, including MaxCut ($M=2$) or
Max-$k$-Cut karp1972reducibility; barahona1988application;
goemans1995improved; deza1994applications; hartmann1996cluster;
shi2000normalized; mei2017solving; frieze1995improved;
goemans2004approximation; so2007approximating, MaxSAT ($M=2$)
goemans1995improved, maximal margin classification ($M=2$)
gieseke2013polynomial, semi-supervised learning ($M=2$) wang2013semi,
correlation clustering on general and random dot-product graphs ($M=2$)
veldt2017correlation, community detection ($M=2$) hajek2016achieving;
abbe2018community, and (quantized) phase synchronization in communication
networks boumal2016nonconvex; zhong2018near. Other problems that can be cast
as (quantized) quadratic form optimization over the unit complex torus include
phase recovery waldspurger2015phase, angular synchronization
singer2011angular, and optimization problems in communication systems
heath1998simple; love2003grassmannian; motedayen2006optimal;
kyrillidis2014fixed; kyrillidis2011rank.

Typically, solving (2) is computationally expensive due to the _presence of
discrete structures_. Discrete algorithmic solutions have been developed to
solve (2) with near-optimal performance and rigorous theoretical guarantees by
exploiting the problem’s structure marti2009advanced; optimization2014inc;
mosek2015mosek; krislock2017biqcrunch. However, the go-to techniques for
solving (2) are oftentimes continuous. For example, methods exist that relax
combinatorial entities in (2) and find a problem solution with non-linear,
continuous optimization techniques. Such techniques can be roughly categorized
into convex grotschel2012geometric; gartner2012approximation;
nesterov1989self; nesterov1994interior; tran2014inexact; tran2016single and
non-convex (_but still continuous_) approaches burer2003nonlinear;
bhojanapalli2018smoothed; bhojanapalli2016dropping; kyrillidis2018provable;
wang2017mixing.

More recently, numerous algorithms have been proposed for solving large-scale
SDP instances. wang2017mixing present a low-rank coordinate descent approach
to solving MaxCut and MaxSat SDP instances that is 10-100$\times$ faster than
state-of-the-art solvers. Similarly, yurtsever2019scalable solve MaxCut
instances on a laptop with almost $8$ million vertices and $10^{13}$ matrix
entries, while goto2019combinatorial propose a simulated bifurcation algorithm
that can solve all-to-all, $100,000$-node MaxCut problems with continuous
weights. Despite such developments, however, SDPs are still hard problems to
solve in practice. _This work aims to contribute along the path of efficient
algorithms for solving large-scale SDPs._

This Paper. We present an algorithmic prototype for solving diagonally
constrained SDP problems at scale. Based on the property of SDP that admits
rank-$k$ solutions alizadeh1997complementarity; barvinok1995problems;
pataki1998rank as detailed in Lemma 2.1, we focus on solving a non-convex, but
equivalent, formulation of (1):

$\displaystyle\min_{V\in\mathbb{R}^{k\times n}}f(V)\triangleq$
$\displaystyle\left\langle C,V^{\top}V\right\rangle$ s.t.
$\displaystyle\|v_{i}\|_{2}=1,~{}\forall i.$ (3)

Here, $V\in\mathbb{R}^{k\times n}$ is the optimization variable, and the set
of linear constraints on $X$ in (1) translate to $V\in(S^{k-1})^{n}$ (i.e.,
each column vector $v_{i}$ satisfies $\|v_{i}\|_{2}=1$). Our algorithm, named
Mixing Method++, adds momentum to the coordinate power iteration method for
solving SDPs by incorporating ideas from accelerated convex optimization and
matrix factorization. As its predecessor wang2017mixing, the algorithm is
simple to implement, requiring only one additional hyperparameter (momentum).
From an empirical perspective, we test the algorithm on MaxCut, MaxSAT, and
MIMO signal detection applications, where we demonstrate significant speedups
in comparison to other state-of-the-art solvers (e.g., $9\times$ to
$10^{2}\times$ speedup on MaxSAT instances) with comparable or improved
solution quality. Additionally, we prove that our accelerated method admits
local linear convergence in the neighborhood of the optimum and always
converges to a first order stationary point. To the best of our knowledge,
this is the first theoretical result that incorporates momentum into a
coordinate-descent approach for solving SDPs.

##  2 Background

Notation. We use $\mathbb{S}^{n}$ to denote the cone of positive semi-definite
(PSD) matrices. We denote the $k$-dimensional sphere by $S^{k-1}$, where $k-1$
indicates the dimension of the manifold. The product of $n$, $k$-dimensional
spheres is denoted by $(S^{k-1})^{n}$. For any $v\in\mathbb{R}^{n}$ and
$M\in\mathbb{R}^{n\times m}$, we define $\|v\|_{2}$ and $\|M\|_{F}$ as the
$\ell_{2}$ and Frobenius norms. For a matrix $V\in\mathbb{R}^{k\times n}$, we
denote its $i$-th column interchangeably as $V_{:,i}$ or $v_{i}$. For
$a\in\mathbb{R}^{n}$, we define $D_{a}$ as the $n\times n$ diagonal matrix
with entries $(D_{a})_{ii}=a_{i}$. We denote the vector of all $1$’s as
$\mathds{1}$ and the Hadamard product as $\odot$. Singular values are denoted
as $\sigma(\cdot)$, while $\sigma_{\texttt{nnz}}(\cdot)$ represents minimum
non-zero singular values. The matrix $C$ corresponds to the cost matrix of the
optimization problem (1).

Low-rank Property of SDPs. Focusing on (1), the number of linear constraints
on $X$ is far less than the number of variables within $X$. This constitutes
such SDPs as _weakly constrained_ , which yields the following result
barvinok1995problems; pataki1998rank.

######  Lemma 2.1.

The SDP problem in (1) has a solution with rank
$k=\left\lceil\sqrt{2n}\,\right\rceil$.

One can enforce low-rank solutions to weakly constrained SDPs by defining
$X=V^{\top}V$ burer2003nonlinear. Here, $V$ belongs to $\mathbb{R}^{k\times
n}$, where $k(k+1)/2>n$ boumal2016nonconvex. This low-rank parameterization,
shown in (3), yields a non-convex formulation of the SDP problem in (1) that
shares the same global solution. This formulation can be solved efficiently
because the conic constraint in (1) is automatically satisfied.222Handling a
PSD constraint requires an eigenvalue decomposition of $n\times n$ matrices,
leading to a $O(n^{3})$ overhead per iteration.

Mixing Method. Consider again the low-rank SDP parameterization given in (3).
To connect this formulation with its discrete form in (2), assume $k=1$, which
yields the following expression:

$\displaystyle\min_{v\in\mathbb{R}^{1\times n}}$ $\displaystyle vCv^{\top}$
s.t. $\displaystyle|v_{i}|=1,~{}\forall i,$ (4)

Notice that $v$ is normalized entry-wise during each optimization step. To
solve (4), one can consider using related algorithms such as power iteration
mises1929praktische, as shown below.

$\displaystyle u$ $\displaystyle\leftarrow Cv$ $\displaystyle v^{+}$
$\displaystyle\leftarrow\texttt{normalize}(u),$

Here, $\texttt{normalize}(u)$ projects $u$ such that $|v^{+}_{i}|=1$ for each
entry of $v^{+}$. This is the crux of the _Mixing Method_ wang2017mixing,
which applies a coordinate power iteration routine to sequentially update each
column of $V$ and solve (3); see Algorithm 1. Variants of Algorithm 1 that
utilize random coordinate selection instead sequential updates for the columns
of $V$ have also been explored erdogdu2018convergence.

In (3), diagonal elements of $C$ do not contribute to the optimization due to
the constraint $\|v_{i}\|_{2}=1$, and thus can be considered as zero. When
considering only the $i$-th column of $V$, The objective $\left\langle
C,V^{\top}V\right\rangle$ takes the form
$2v_{i}^{\top}\sum_{j=1}^{n}c_{ij}v_{j}$. The full gradient of the objective
with respect to $v_{i}$ can be computed in closed form at each iteration,
yielding the coordinate power iteration routine outlined in Algorithm 1.
Interestingly, this procedure also comes with theoretical guarantees
wang2017mixing, which we include below for completeness. 333Most of global
guarantees in wang2017mixing hold for a variant of the Mixing Method that
resembles coordinate gradient descent with step size, rather than the
coordinate power iteration method.

######  Theorem 2.1.

The Mixing Method converges linearly to the global optimum of (3), under the
assumption that the initial point is close enough to the optimal solution.

Algorithm 1 Mixing Method wang2017mixing

1: Input: $C\in\mathbb{S}^{n}$, $V\in(S^{k-1})^{n}$.

2: while not yet converged do

3: for $i=1$ to $n$ do

4: $v_{i}\leftarrow\texttt{normalize}\left(-\sum_{j\neq i}c_{ij}v_{j}\right)$

5: end for

6: end while

Algorithm 2 Mixing Method++

1: Input: $C\in\mathbb{S}^{n}$, $V\in(S^{k-1})^{n}$, $0\leq\beta<1$.

2: while not yet converged do

3: for $i=1$ to $n$ do

4: $u_{i}\leftarrow\texttt{normalize}\left(-\sum_{j\neq i}c_{ij}v_{j}\right)$

5:
$v_{i}\leftarrow\texttt{normalize}\left(u_{i}+\beta\left(u_{i}-v_{i}\right)\right)$

6: end for

7: end while

##  3 Mixing Method++

Towards a Momentum-inspired Update Rule. In this work, we introduce Mixing
Method++, which significantly improves upon the performance of Mixing Method
with the addition of a momentum step. Our goal is to employ acceleration
techniques polyak1987introduction; nesterov2013introductory to improve the
overall performance of Mixing Method. For gradient descent, a classical
acceleration technique is the Heavy-Ball method polyak1987introduction, which
iterates as follows:

$\displaystyle w^{t+1}=w^{t}-\eta\nabla f(w^{t})+\beta(w^{t}-w^{t-1}).$

Here, $w$ is the optimization variable, $f(\cdot)$ is a differentiable loss
function, $\eta$ is the step size, and $\beta$ is the momentum parameter.
Intuitively, the heavy-ball method exploits the history of previous updates by
moving towards the direction $w^{t}-w^{t-1}$, weighted by the momentum
parameter $\beta$. If we naively incorporate a similar momentum term into the
power iteration, we arrive at the following recursion:

$\displaystyle
w^{t+1}=\texttt{normalize}\left(Cw^{t}+\beta(Cw^{t}-w^{t})\right),$

where there is no notion of a step size $\eta$. By adapting this accelerated
power iteration scheme to the Mixing Method update rule, we arrive at the
following formulation:

$\displaystyle\hat{v}_{i}=\texttt{normalize}(g_{i}+\beta(g_{i}-v_{i})),$

where we denote $g_{i}=-\sum_{j\neq i}c_{ij}v_{j}$. Because $g_{i}$ is not
normalized before the addition of the momentum term, it is possible that
$g_{i}$ and $v_{i}$ are of significantly different magnitude. As such, we
consider the normalization of $g_{i}$ as an intermediate step and adopt a
double-projection approach per iteration; see Algorithm 2. Such double
projection ensures that $u_{i}$ and $v_{i}$ are of comparable magnitude when
the momentum term is added, resulting in significantly-improved momentum
inertia.

Properties of Mixing Method++. The proposed algorithm is comprised of a two-
step update procedure. The last normalization projects
$u_{i}+\beta\left(u_{i}-v_{i}\right)$ onto the unit sphere, which ensures the
updated matrix is feasible.In contrast, the first normalization ensures that
$-\sum_{j=1}^{n}c_{ij}v_{j}$ and $v_{i}$ are comparable in magnitude, as the
former may have a norm greater than $1$. The overall computational cost of the
added momentum step is negligible, making the step-wise computational
complexity of Mixing Method++ roughly equal to that of the Mixing Method.

We show that Mixing Method++ converges linearly when $V$ lies in a
neighborhood of the global optimum. Despite the similarity of our convergence
rates to previous work wang2017mixing, achieving acceleration in theory is not
always feasible, even in convex, non-stochastic cases; see
assran2020convergence; devolder2014first; ghadimi2015global;
polyak1987introduction; lessard2016analysis; loizou2017momentum. However,
acceleration can oftentimes be observed in practice nonetheless. Our numerical
experiments empirically reveal the presence of acceleration, and Mixing
Method++ with $\beta=0.8$ yields a robust algorithm with significant speed-ups
in comparison to Mixing Method. This difference is particularly noticeable in
large instances of relaxations for (2), as illustrated in Section 5.

##  4 Convergence properties of Mixing Method++

We begin by presenting some notation for brevity and readability. $f^{\star}$
is used to denote the optimal value of the objective function $f(V)=\langle
C,V^{T}V\rangle$. The matrices $V$ and $\hat{V}$ refer to the current and the
resulting iterates from the inner iteration of Algorithm 2. For each
$i\in[n]$, we define the family of maps $Z_{i}:\mathbb{R}^{k\times
n}\to\mathbb{R}^{k\times n}$ as follows.

$V\longmapsto\begin{pmatrix}|&&|&|&|&&|\\\
\hat{v_{1}}&\ldots&\hat{v}_{i-1}&v_{i}&v_{i+1}&\ldots&v_{n}\\\
|&&|&|&|&&|\end{pmatrix}.$ (5)

In words, $Z_{i}(V)$ corresponds to the partially updated $V$ matrix (i.e.,
within the inner loop of Algorithm 2) before updating $v_{i}$ to
$\hat{v}_{i}$. Thus, only the vector $v_{i}$ changes from $Z_{i}(V)$ to
$Z_{i+1}(V)$. Notice that the last iteration, $Z_{n+1}(V)$, corresponds to the
end of the inner cycle such that $Z_{n+1}(V)=\hat{V}$.

We now present our main theoretical results. Due to space limitations, we
defer full proofs to the Appendix and provide only a sketch of Theorem 4.1. We
adopt Assumption 3.1 from wang2017mixing.

######  Assumption 4.1.

For $i\in[n]$, assume $\|\sum_{j=1}^{n}c_{ij}v_{j}\|_{2}$ do not degenerate in
the procedure. That is, all norms are always greater than or equal to a
constant value $\delta>0$.

Put simply, Assumption 4.1 states that the update rule will never yield
$u_{i}$ with zero norm. We find that this assumption consistently holds in
practice for both Mixing Method and Mixing Method++ applied over several
objectives.444We provide empirical analysis of Assumption 4.1 in Appendix C.
Note that if Assumption 4.1 does not hold, we can exploit the simple trick No
further assumptions are required for the second normalization step of the
algorithm provided that $0\leq\beta<1$, as shown by the following observation.

######  Observation 4.1.

$\|u_{i}-\beta v_{i}\|_{2}\not=0$, whenever $0\leq\beta<1$.

###### Proof.

Based on upper and lower triangle inequality bounds and the fact that $u_{i}$
and $v_{i}$ are both in $S^{k-1}$, we get the inequality
$1\leq\left\lVert\left(1+\beta\right)u_{i}-\beta v_{i}\right\rVert_{2}\leq
1+2\beta$. The result follows immediately. ∎

Our first technical lemma lower bounds the decrease in the objective over
inner iterations of Algorithm 2, thus ensuring that the algorithm makes
progress as $V$ becomes $\hat{V}$.

######  Lemma 4.1.

Let $\hat{V},V$ be the next and previous iteration of Algorithm (2),
respectively. $y_{i}=\|g_{i}\|_{2}$, where $g_{i}$ is a linear combination of
the $i$th row of $C$ and the columns in $V$ (i.e., $g_{i}=-\sum_{j\neq
i}c_{ij}v_{j}$); see Appendix A. $w_{i}=\|u_{i}-\beta v_{i}\|_{2}$. Then,

$\begin{split}f(V)-f(\hat{V})\geq\frac{1-\beta}{1+\beta}\sum_{i=1}^{n}y_{i}\,\|v_{i}-\hat{v}_{i}\|_{2}^{2}.\end{split}$
(6)

Furthermore, $f$ is non-increasing whenever $0\leq\beta<1$.

Next, we lower bound the magnitude of the update from $V$ to $\hat{V}$ based
on the objective residual, showing that the magnitude of our update will be
large until we approach an optimum.

######  Lemma 4.2.

Under the assumption 4.1, there exists positive constants $\gamma,\zeta$ such
that

$\|V-\hat{V}\|_{F}^{2}\geq\big{(}\zeta-\gamma\,\|\Delta\|_{F}\big{)}\,(f(V)-f^{*}),$

where $\Delta:=(D_{w}D_{y}+\beta D_{y})-(D_{w^{*}}D_{y^{*}}+\beta D_{y^{*}})$.

Finally, we upper bound the value of $\|\Delta\|^{2}_{2}$ with respect to the
objective residual, which later enables the definition of a neighborhood for
local linear convergence within Theorem 4.1.

######  Lemma 4.3.

Under Assumption 4.1, there exists a constant $\tau>0$ such that

$\|(w+\beta\mathds{1})\odot y-(w^{*}+\beta\mathds{1})\odot
y^{*}\|_{2}^{2}\leq\tau(f(V)-f^{*}).$

By drawing upon the results in our technical lemmas, we then prove that Mixing
Method++ achieves linear convergence within a local neighborhood of the
optimum.

######  Theorem 4.1 (Local linear convergence).

Let $f^{\star}$ represent the optimal value of the objective function, and let
$\delta$ be the non-degenerative lower-bound of Assumption 4.1. Define the
neighborhood,

$\mathcal{N}_{\star}:=\left\\{V\in(S^{k-1})^{n}\,:\,\tau(f(V)-f^{\star})\leq\left(\tfrac{\zeta-\kappa}{\gamma}\right)^{2}\right\\},$

for constant values $\zeta,\tau,\gamma$ and $\kappa$ presented in Lemmas 4.1,
4.2, and 4.3. If $V\in\mathcal{N}_{\star}$, we have

$f(\hat{V})-f^{\star}\leq\left(1-\rho\right)\,\left(f(V)-f^{\star}\right),$

where $\hat{V}$ represents the new state of $V$ after performing the inner
iterations of _Mixing Method++_ and $\rho\in[0,1]$ is a constant. This result
shows that Algorithm (2) decreases the distance to the optimal value
$f^{\star}$ linearly when $V\in\mathcal{N}_{\star}$.

Proof sketch. We now overview the proof for Theorem 4.1; see Appendix A for
the full proof.

  * •

Lemma 4.1 shows that the objective function is non-increasing provided that
$\beta$ is bounded by $1$. Therefore,
$f(V)-f(\hat{V})\geq\delta\frac{1-\beta}{1+\beta}\|V-\hat{V}\|_{F}$. This
property provides the theoretical bound that restricts the search space for
momentum $\beta<1$.

  * •

Using the definition of $S^{*}=\left(1+\beta\right)\left(C+D_{y^{*}}\right)$,
the fact that $S^{*}\succeq 0$, and lower bounds on the trace of matrix
products, we invoke Lemma 4.2 to derive the following lower bound
$\|V-\hat{V}\|_{F}\geq(\zeta-\gamma\,\|\Delta\|_{F})\,(f(V)-f^{*})$.

  * •

Let $V$ be in the neighborhood $\mathcal{N}_{*}$ for some positive constant
$\kappa$. Given that the objective is non-increasing, all following iterations
of $V$ derived from Mixing Method++ must lie in $\mathcal{N}_{*}$. More
specifically, the $k$-folding ($k$-th composition) of $Z_{n+1}$ on $V$ lies in
$\mathcal{N}_{*}$ for all $k\geq 0$ if $V\in\mathcal{N}_{\star}$, implying
that $Z_{n+1}(V)\in\mathcal{N}_{\star}$ when $V\in\mathcal{N}_{\star}$.

  * •

Using Lemma 4.3 and the fact that $V$ is in $\mathcal{N}_{*}$, we show that
$\zeta-\gamma\,\|\Delta\|_{F}\geq\kappa$.

  * •

Using inequalities in the arguments above, we derive
$f(V)-f(\hat{V})\geq\delta\kappa\frac{1-\beta}{1+\beta}(f(V)-f^{*})$, which by
algebraic manipulation implies the following, where we substitute
$\rho=\delta\kappa\frac{1-\beta}{1+\beta}$ to yield the final result.

$\left(1-\delta\kappa\frac{1-\beta}{1+\beta}\right)(f(V)-f^{*})\geq
f(\hat{V})-f^{*}.$

##  5 Experimental Results

We test Mixing Method++ on three, well-known applications of the objective
(1): MaxCut, MaxSAT, and MIMO signal detection. We use the same formulation in
wang2017mixing to solve MaxCut and MaxSAT, while for MIMO we follow the
experimental setup of liu2017discrete. We implement Mixing Method++ in C with
$\beta=0.8$. 555We extensively tested our theoretical bound on $\beta$ in
experiments. After testing all values in the set $[0,1)$ with an interval of
$0.01$, we observe that $\beta=0.8$ indeed yields the best performance on over
90% of instances. Each experiment is run on a single core in a homogeneous
Linux cluster with 2.63 GHz CPU and 16 GB of RAM. For MaxCut, we compare
Mixing Method++ to the CGAL CGAL, Sedumi sturm1999using, SDPNAL+
yang2015sdpnal, SDPLR burer2003nonlinear, MoSeK aps2019mosek, and Mixing
Method wang2017mixing solvers, while for MaxSAT we compare against the Mixing
Method and Loandra loandra. For MIMO signal detection, we compare the
performance of Mixing Method++ to that of the Mixing Method.

We often report performance in terms of the objective residual, which is
defined as the difference between objective values given by each solver and
the lowest objective value achieved by any solver. More formally, let
$\texttt{Obj}_{ij}$ be the value of the objective function given by optimizer
$i$ on SDP instance $j$. Then, the objective residual of optimizer $i$ on
instance $j$ is defined as
$|\texttt{Obj}_{ij}-\mathop{\min}\limits_{i^{\prime}\in
O}\texttt{Obj}_{i^{\prime}j}|$, where $O$ is the set of all solvers used on
problem instance $j$.

###  5.1 MaxCut

Solver |  | # Solved instances |  | median(log(objective residual)) |  |  Acceleration ($\times$)
---|---|---|---|---|---|---
Sedumi |  | 60 |  | -4.30 |  | 360.20
MoSeK |  | 74 |  | -3.24 |  | 348.37
SDPNAL+ |  | 52 |  | -4.33 |  | 1316.35
CGAL |  | 92 |  | 0.93 |  | 49.19
SDPLR |  | 95 |  | -2.89 |  | 9.14
Mixing Method |  | 106 |  | -2.80 |  | 5.26
Mixing Method++ |  | 111 |  | -3.87 |  | 1
Table 1: Results on 203 MaxCut instances. Mixing Method++ solved most
instances efficiently with decent precision. E.g., Mixing Method++ is 5.26 and
360.20 times faster than Mixing Method and Sedumi, respectively. The median of
log(objective residual) of Mixing Method++ is -3.87, meaning that the median
precision of Mixing Method++ on all cases is $10^{-3.87}$.

The SDP relaxation of the MaxCut problem is given by

$\underset{\|v_{i}\|=1,}{\mathrm{max}}\,\frac{1}{2}\sum_{ij}c_{ij}\left(\tfrac{1-v_{i}^{\top}v_{j}}{2}\right),$

where $C$ is the adjacency matrix. Thus, we aim to minimize the following
objective function

$f_{\text{MaxCut}}=\left\langle C,\tfrac{\mathds{1}_{n\times
n}-V^{\top}V}{4}\right\rangle~{}~{}\text{s.t.}~{}~{}\|v_{i}\|_{2}=1~{}~{}\forall
i,$

where $V\in\mathbb{R}^{k\times n}$ and $\mathds{1}_{n\times n}$ is the
$n\times n$ matrix of all $1$s. We consider the following benchmarks within
our MaxCut experiments.

  * •

Gset: 67 binary-valued matrices, produced by an autonomous, random graph
generator. The dimension $n$ varies from $800$ to $10^{4}$.

  * •

Dimacs10: 152 symmetric matrices with $n$ varying from $39$ to $50,912,018$ ,
chosen for the $10$-th Dimacs Implementation Challenge. We consider the 136
instances with dimension $n\leq 8\cdot 10^{6}$.

Figure 1: G40: $n=2000$, number of edges $=11766$. Mixing Method++ converges
within $1$s while Mixing Method uses about $4$s. MoSeK and Sedumi need
$10^{2}$ and $10^{3}$ seconds respectively to converge.

Figure 2: uk: $n=4824$, number of edges$=6837$. Consider the CPU time when
each solver reaches objective residual $10^{-2}$. Mixing Method++ uses only
20s while Mixing Method needs 200s. MoSeK and Sedumi consume more than
$10^{3}$ and $10^{4}$ seconds respectively.

Table 1 lists the number of solved instances, median logarithm of objective
residual, and the pairwise acceleration ratio for different SDP solvers tested
on MaxCut.666When an instance is solved by both Mixing Method++ and solver $s$
within the 24 hour time limit, the pairwise acceleration ratio is defined as
the CPU time taken by solver $s$ divided by that of Mixing Method++. In
comparison to the other solvers, Mixing Method++ solves the most MaxCut
instances (i.e., 111 of 203 total instances) within the imposed 24 hour time
limit and provides a $5\times$ to $1316\times$ speedup. Moreover, the
efficiency of Mixing Method++ does not harm its precision, which is evident in
its median logarithm objective residual of $-3.87$. Such performance is
preceded only by Sedumi $(-4.30)$ and SDPNAL+ $(-4.33)$, which are each more
than $300\times$ slower than Mixing Method++.

We also study two specific cases of MaxCut: G40 from Gset and uk from
Dimacs10. In Figures 1 and 2, we plot objective residual with respect to CPU
time of Mixing Method++ and other SDP solvers on the G40 and uk instances,
respectively.777We do not include CGAL and SDPLR in the plots, because the
intermediate solution given by those methods is infeasible, which makes
intermediate objective values difficult to quantify. In both MaxCut instances,
Mixing Method++ converges significantly faster than the other solvers that
were considered. Supplemental results on the MaxCut objective are provided in
Appendix D.1.

###  5.2 MaxSAT

Methods |  | Avg. Approx. Ratio
---|---|---
Loandra |  | 0.945
Mixing Method |  | 0.975
Mixing Method++ |  | 0.977
Table 2: Results on MaxSAT. Mixing Method++ achieves the best approximation
ratio with a time limit of 300s. 850 problems are tested in total. Problems
are derived from crafted (331), random (454), and industrial domains (65).

Let $s_{ij}\in\\{\pm 1,0\\}$ be the sign of variable $i$ in clause $j$. Then,
solving the problem below provides an upper bound to the exact MaxSAT
solution.

$\underset{\|v_{i}\|=1,}{\mathrm{max}}\sum_{j=1}^{m}\left(1-\tfrac{\|Vs_{j}\|^{2}-(|s_{j}|-1)^{2}}{4|s_{j}|}\right),\vspace{-0.05cm}$

where $s_{j}=[s_{1j},...,s_{nj}]^{\top}$. Therefore, the MaxSAT SDP relaxation
can then be solved by minimizing the following, related objective.

$f_{\text{MaxSAT}}=\left\langle
C,{V^{\top}V}\right\rangle~{}~{}\text{s.t.}~{}~{}\|v_{i}\|_{2}=1~{}~{}\forall
i,$

where $V\in\mathbb{R}^{k\times n}$ and
$C=\sum_{j=1}^{m}\frac{s_{j}s_{j}^{\top}}{4|s_{j}|}$.

We consider 850 problem instances from the 2016-2019 MaxSAT competition within
our experiments. These benchmarks are categorized into three tracks: _random,
crafted_ and _industrial_. We evaluate Mixing Method++ as a partial MaxSAT
solver888A partial solver is only required to provide the assignment it finds
with the least number of violated clauses, while a complete solver must also
provide proof that such an assignment is optimal. and compare its performance
to that of the Mixing Method and Loandra, which was the best partial solver in
the 2019 MaxSAT competition. For cases where the optimal solution is not
known, we use the best solution given by all solvers as ground truth.

The average approximation ratio of each solver across all MaxSAT instances is
shown in Table 2. Mixing Method++ achieves the best average approximation
ratio. We also consider two specific MaxSat instances: s3v90c800-6 from the
random track and wb.4m8s4.dimacs.filtered from the industrial track. In
Figures 3 and 4, we plot both the objective residual and cost (i.e., the
number of unsatisfied clauses) with respect to CPU time for s3v90c800-6 and
wb.4m8s4.dimacs.filtered instances, respectively. In both cases, Mixing
Method++ converges faster than vanilla Mixing Method. Furthermore, Mixing
Method++ returns the best solution of all solvers on the s3v90c800-6 instance.
Although Loandra provides the optimal solution to the wb.4m8s4.dimacs.filtered
quickly, our method still produces a better solution than the Mixing Method in
less time. Supplemental results on the MaxSAT objective are provided in
Appendix D.2.

Figure 3: s3v90c800-6.cnf ($n=90$, $m=800$). Mixing Method++ not only
converges faster than Mixing Method, but provides better solution (30)
compared to Mixing Method (31) and Loandra (33). The optimal solution has cost
26.

Figure 4: wb_4m8s4_dimacs_filtered.cnf ($n=463,080$, $m=1,759,150$). Loandra
reaches optimal (230) in 1 sec, while Mixing Method and Mixing Method++ return
a solution with cost $28,306$ and $25,210$ respectively in 1 hour.

###  5.3 MIMO

The MIMO signal detection setting is defined as follows:

$\operatorname*{arg\,min}\limits_{x\in\mathcal{S}}\|y-Hx\|_{2}^{2},$

where $y\in\mathbb{R}^{m}$ is the signal vector, $H\in\mathbb{R}^{m\times n}$
is the channel matrix, and $x\in\mathcal{S}$ is the transmitted signal. It is
assumed that $y=Hx+v$, where $v\sim\mathcal{N}(0,\sigma_{v}^{2})$. For our
experiments, we consider the constellation $\\{\pm 1\\}$, which yields
$\mathcal{S}=\\{\pm 1\\}^{n}$. The MIMO SDP relaxation can then be solved by
minimizing the following objective.

$f_{\text{MIMO}}=\langle
C,V^{T}V\rangle~{}~{}\text{s.t.}~{}~{}\|v_{i}\|_{2}=1~{}~{}\forall i,$ (7)

where $V\in\mathbb{R}^{k\times(n+1)}$ and $C=\begin{bmatrix}H^{T}H&-H^{T}y\\\
-y^{T}H&y^{T}y\end{bmatrix}$.

We compare Mixing Method++ to the Mixing Method on several simulated instances
of the MIMO SDP objective. In all cases, $H$ is sampled i.i.d. from a standard
normal distribution, and each coordinate of $x$ is sampled independently from
the constellation $\\{\pm 1\\}$. Our noise vector $v$ is sampled from a
centered normal distribution with variance $\sigma_{v}^{2}$. Experiments are
performed for numerous different settings of $\sigma_{v}^{2}$, determined by
the signal-to-noise ratio (SNR) (i.e.,
$\sigma_{v}^{2}=\frac{mn}{\text{SNR}}$). We perform tests with
$\text{SNR}\in\\{8,16\\}$ and problem sizes
$(m,n)\in\\{(16,16),(32,32),(64,32)\\}$. We plot the objective residual
achieved by both Mixing Method and Mixing Method++ on different MIMO instances
with respect to CPU time in Fig. 5. As can be seen, Mixing Method++ matches or
exceeds the performance of the Mixing Method in all experimental settings.

Figure 5: Convergence comparison of Mixing Method++ with $\beta=0.8$ and
Mixing Method on MIMO signal detection problems with different SNRs. The size
of the channel matrix is listed above each subplot.

###  5.4 The Effect of Sparsity

Interestingly, we empirically observe that Mixing Method++ achieves different
levels of acceleration based on the sparsity of the objective matrix $C$. To
further understand the impact of sparsity on the performance of Mixing
Method++, we compare Mixing Method++ to the Mixing Method for objectives with
different sparsity levels. 999Here, we define sparsity as the number of zero-
valued elements divided by the total number of elements. In particular, we
construct a MaxCut objective for a graph with 500 nodes and consider adjacency
matrices with three different sparsity levels: $0.8$, $0.5$, and $0.2$. In
Figure 6, we plot the objective residual of both Mixing Method and Mixing
Method++ on these problem instances. As can be observed, the acceleration
provided by Mixing Method++ is more significant when $C$ is sparse.
Furthermore, decreasing the value of $\beta$ when $C$ is dense seems to
improve the performance of Mixing Method++ when tested on instances of the
MaxCut objective.

Figure 6: Convergence comparison of Mixing Method++ (blue line) with
$\beta=0.8$ and Mixing Method (orange line) on 500-node MaxCut problem
instances with different sparsity levels.

###  5.5 Adaptive Momentum

Within our experiments, we observe that Mixing Method tends to decrease the
objective residual faster than Mixing Method++ in the early stage iterations.
Inspired by this observation and previous works that propose adaptive momentum
schedules o2015adaptive; chen2019demon; wang2020scheduled, we explore whether
dynamically adjusting the value of $\beta$ between iterations can further
improve the performance of Mixing Method++. We empirically observe that the
following momentum schedule yields the best performance:

$\displaystyle\beta_{t}=\beta\left(1-\exp(-\alpha t/T)\right),$

where $\beta=0.8$ is the default momentum value, $t$ is the current iteration,
$T$ is the total number of iterations, and $\alpha$ is a constant such that
$\exp(-\alpha)\approx 0$. Intuitively, this schedule begins with a low
momentum and quickly increases the momentum value to $\beta=0.8$, thus forming
a warm-up schedule for momentum. As illustrated in Figure 7, such a momentum
schedule improves the convergence speed of the algorithm.

Figure 7: Comparison of Mixing Method (orange line), Mixing Method++ with
fixed momentum $\beta=0.8$ (blue line), Mixing Method++ with an increasing
Scheduled Momentum from $0$ to $0.8$ (green line) in terms of the objective
cost per iteration.

##  6 Conclusion and Future Directions

We present a novel approach, Mixing Method++, to solve diagonally constrained
SDPs. Mixing Method++ inherits the simplicity—and non-trivially preserves
theoretical guarantees from—its predecessor. In practice, it yields not only
more efficient convergence on nearly all tested instances, but also
improvements in solution quality. Mixing Method++ adds one addition
hyperparameter $\beta$ for which we provide a theoretical upper bound (i.e.,
$\beta<1$). Using $\beta=0.8$ experimentally leads to a robust algorithm,
which outperforms Mixing Method and numerous other state-of-the-art SDP
solvers. Theoretically validating the experimental acceleration provided by
Mixing Method++ is still an open problem, which could potentially be handled
with the help of Lyapunov Analysis wilson2016lyapunov; jin2018accelerated.

##  Appendix A Local linear convergence

We now give detailed proof of Theorem 4.1, which was sketched in Section 4.

Let $g_{i}$ be the linear combination in step 4 of Algorithm 1, before
projecting onto the unit sphere:

$\begin{split}g_{i}&=\sum_{j<i}c_{ij}\hat{v}_{j}+\sum_{j>i}c_{ij}v_{j}.\end{split}$

We also define $y_{i}=\|g_{i}\|_{2}$ and $f(V)\triangleq\langle
C,V^{T}V\rangle$ as the objective value for some viable solution $V$.
Therefore, step 4 of Algorithm 1 can be expressed as: $y_{i}u_{i}=-g_{i}$.
Thus, the changes presented in the Algorithm 2 can be written as the
combination of 4 and the momentum, in step 5:

$\begin{split}u_{i}=-\frac{g_{i}}{y_{i}},\hskip 10.0pt\text{ and }\hskip
10.0pt\hat{v}_{i}=\frac{\left(1+\beta\right)u_{i}-\beta
v_{i}}{w_{i}},\end{split}$ (8)

where $w_{i}=\left\lVert\left(1+\beta\right)u_{i}-\beta
v_{i}\right\rVert_{2}$. This leads to the following one line relation,

$w_{i}y_{i}\hat{v}_{i}=-\left(1+\beta\right)g_{i}-\beta y_{i}v_{i},~{}\forall
i\in[n].$ (9)

Herein, the matrices $V$ and $\hat{V}$ refer to the current and the end state
of the inner iteration in algorithm 2. For each $i\in[n]$, we define the
family of maps $Z_{i}:(S^{k-1})^{n}\to(S^{k-1})^{n}$ as follows

$Z_{i}(V)\,e_{j}=\begin{cases}\hat{v}_{j},\hskip 10.0pt\text{if $j<i$}\\\
v_{j},\hskip 10.0pt\text{if $j\geq i$},\end{cases}$ (10)

where $\\{e_{j}\\}$ is the canonical basis in $\mathbb{R}^{n}$. It is trivial
to see that $Z_{i}(V)$ corresponds to the matrix in the inner cycle before
updating $v_{i}$ to $\hat{v}_{i}$. Thus, only the vector $v_{i}$ has been
changed from $Z_{i}(V)$ to $Z_{i+1}(V)$. Also notice that $Z_{n+1}(V)$
corresponds to the end of the inner cycle, i.e., $Z_{n+1}(V)=\hat{V}$.

### Matrix Representation

The transition from the state $v_{i}$ to $\hat{v}_{i}$ is expressed in the
relation (9), then for any $i\in[n]$ in each inner cycle we have,

$\displaystyle\begin{split}\left(1+\beta\right)\sum_{j<i}c_{ij}\hat{v}_{j}+w_{i}y_{i}\hat{v}_{i}=-\left(1+\beta\right)\sum_{j>i}c_{ij}v_{j}-\beta
y_{i}v_{i}.\end{split}$ (11)

Let $L$ be the strictly lower triangular matrix of $C$ so (11) leads us to the
representation

$\displaystyle\left(\left(1+\beta\right)L+D_{w}D_{y}\right)\,\hat{V}^{\top}=-\left(\left(1+\beta\right)L^{\top}+\beta
D_{y}\right)\,V^{\top}.$ (12)

### Proof of Lemma 4.1

###### Proof.

Let $\beta$ be in $[0,1]$. When fixing all the others variables, $f$ with
respect to $v_{i}$ is given by

$\langle
C,V^{\top}V\rangle=2v_{i}^{\top}\left(\sum_{j}c_{ij}v_{j}\right)+\text{constant}.$

Since only $v_{i}$ is changed, during the $i$th coordinate-wise updating, the
only part of the objective that will change is $2v_{i}^{\top}g_{i}$. That is,

$f(Z_{i}(V))-f(Z_{i+1}(V))=2g_{i}^{\top}\left(v_{i}-\hat{v}_{i}\right).$

The updating rule of algorithm 2 is given by the recursion (9), so using this
into the above equation yields

$\begin{split}f(Z_{i}(V))-f(Z_{i+1}(V))&=\frac{2y_{i}\,(w_{i}-\beta)}{1+\beta}(1-v_{i}^{\top}\hat{v}_{i}),\\\
&=\frac{y_{i}\,(w_{i}-\beta)}{1+\beta}\|v_{i}-\hat{v}_{i}\|_{2}^{2}.\end{split}$

From observation 4.1, we have $w_{i}-\beta\geq 1-\beta$. Thus,

$\begin{split}f(Z_{i}(V))-f(Z_{i+1}(V))\geq
y_{i}\,\frac{1-\beta}{1+\beta}\|v_{i}-\hat{v}_{i}\|_{2}^{2}.\end{split}$

Adding over all the $i$, leads us to the desired conclusion

$\begin{split}f(V)-f(\hat{V})&=\sum_{i=1}^{n}\,f(Z_{i}(V))-f(Z_{i+1}(V))\\\
&\geq\frac{1-\beta}{1+\beta}\sum_{i=1}^{n}y_{i}\,\|v_{i}-\hat{v}_{i}\|_{2}^{2}.\end{split}$

∎

The following observation will be useful later on.

######  Observation A.1.

For all $i\in[n]$, $w_{i}^{*}=1$

###### Proof.

Since $v^{*}_{i}$ is a fixed point, using this fact in Eq. (9) yields
$(w^{*}_{i}+\beta)v_{i}^{*}=u_{i}^{*}$. Taking the norm,
$w_{i}^{*}+\beta=1+\beta$, since $u^{*}_{i}$ and $v^{*}_{i}$ are in the unit
sphere. ∎

Let $S$ be defined as $\left(1+\beta\right)C+D_{w}D_{y}+\beta D_{y}$ and let
$S^{*}$ be $\left(1+\beta\right)C+D_{w^{*}}D_{y^{*}}+\beta D_{y^{*}}$; namely,
the corresponding value of $S$ when $V$ is optimal in the optimization (3). We
have the following observation.

######  Observation A.2.

$S^{*}$ is PSD.

###### Proof.

By observation A.1 we have $w^{*}_{i}=1$, so

$S^{*}=\left(1+\beta\right)C+D_{w^{*}}D_{y^{*}}+\beta
D_{y^{*}}=\left(1+\beta\right)\left(C+D_{y}^{*}\right),$

where according to lemma 3.12, and since $V$ is optimal, it follows that
$\frac{S^{*}}{1+\beta}\succeq 0$, hence $S^{*}\succeq 0$. ∎

### Proof of Lemma 4.2

###### Proof.

By (12), we have

$\displaystyle\hat{V}=\,-V\left(\left(1+\beta\right)L+\beta
D_{y}\right)\left(\left(1+\beta\right)L^{\top}+D_{w}D_{y}\right)^{-1}.$

This leads to

$\displaystyle\begin{split}V-\hat{V}=&V(\left(1+\beta\right)C+D_{w}D_{y}+\beta
D_{y})\\\ &(\left(1+\beta\right)L^{\top}+D_{w}D_{y})^{-1}\\\
=&VS\underbrace{(\left(1+\beta\right)L^{\top}+D_{w}D_{y})^{-1}}_{R},\end{split}$

Using the definition $\Delta\triangleq(D_{w}D_{y}+\beta
D_{y})-(D_{w^{*}}D_{y^{*}}+\beta D_{y^{*}})$, and the fact that
$S=S^{*}+\Delta$ in the previous equation, we get

$\begin{split}\|V-\hat{V}\|_{F}^{2}&\geq\|VS^{*}R\|_{F}^{2}+2\text{tr}\big{(}V^{\top}VS^{*}RR^{\top}\Delta\big{)}.\end{split}$
(13)

We claim that

$\|(w+\beta\mathds{1})\odot y-(w^{*}+\beta\mathds{1})\odot
y^{*}\|_{2}I_{n}+\Delta\succeq 0.$

To see this, first notice $D_{w}D_{y}+\beta D_{y}=D_{(w+\beta\mathds{1})\odot
y}$. Thus,

$\Delta=D_{(w+\beta\mathds{1})\odot y}-D_{(w^{*}+\beta\mathds{1})\odot
y^{*}}.$

Besides,

$\|(w+\beta\mathds{1})\odot y-(w^{*}+\beta\mathds{1})\odot
y^{*}\|_{2}=\|\Delta\|_{F}.$

So the claim follows since $\beta\geq 0$ and termwise

$\begin{split}\|\Delta\|_{F}+\Delta_{i}\geq|\Delta_{i}|+\Delta_{i}\geq
0.\end{split}$

Now, for the first part in (13), we get

$\begin{split}\|VS^{*}R\|_{F}^{2}&=\text{tr}(\underbrace{{S^{*}}^{\top}}_{\succeq
0}\underbrace{V^{\top}V}_{\succeq 0}\underbrace{S^{*}}_{\succeq
0}\underbrace{RR^{\top}}_{\succeq 0})\\\
&\geq\sigma_{\min}^{2}(R)\sigma_{\texttt{nnz}}(S^{*})\text{tr}(V^{T}VS^{*}).\end{split}$

For the second in (13), we get

$\begin{split}\text{tr}\big{(}V^{\top}VS^{*}RR^{\top}\Delta\big{)}&\geq-\|\Delta\|_{F}\operatorname*{\text{Tr}}\big{(}V^{\top}VS^{*}RR^{\top}\big{)}\\\
&\geq-\|\Delta\|_{F}\,\sigma_{\max}^{2}(R)\text{tr}\big{(}V^{\top}VS^{*}\big{)},\end{split}$

this last inequality is given by our previous claim, the fact that
$S^{*}\succeq 0$ (observation A.2) and

$\begin{split}\operatorname*{\text{Tr}}(V^{\top}VS^{*}RR^{\top})&=\operatorname*{\text{Tr}}(V^{\top}VS^{*}\underbrace{UD^{2}U^{\top}}_{RR^{\top}})\\\
&=\operatorname*{\text{Tr}}(V^{\top}VS^{*}D^{2})\\\
&\leq\sigma^{2}_{\max}(R)\,\operatorname*{\text{Tr}}(V^{\top}VS^{*}).\end{split}$

Using the lower bounds in (13) gives us

$\begin{split}\|V&-\hat{V}\|_{F}^{2}\geq\big{(}\sigma_{\min}^{2}(R)\sigma_{\texttt{nnz}}(S^{*})-2\|\Delta\|_{F}\sigma_{\max}^{2}(R)\big{)}\text{tr}\big{(}V^{\top}VS^{*}\big{)}.\end{split}$

To conclude, the term

$\begin{split}\operatorname*{\text{Tr}}(V^{\top}VS^{*})&=\operatorname*{\text{Tr}}(V^{\top}VC)-\operatorname*{\text{Tr}}(-V^{\top}VD_{y^{*}})\\\
&=f(V)-\operatorname*{\text{Tr}}(-V^{\top}VD_{y^{*}})\\\
&\stackrel{{\scriptstyle i)}}{{=}}f(V)+\sum\limits_{i}\underbrace{\langle
v_{i},v_{i}\rangle}_{=1}\cdot y^{*}_{i}\\\
&=f(V)+\langle\mathds{1},y^{*}\rangle\\\ &\stackrel{{\scriptstyle
ii)}}{{\geq}}f(V)-f^{*}.\end{split}$

where $i)$ follows from the fact that $D_{y^{*}}$ is diagonal. To understand
$ii)$, consider the dual to the diagonally constrained SDP problem (1) (i.e.,
see Lemma 3.12 of [wang2017mixing] for more details).

$\displaystyle\max\limits_{y}-\langle\mathds{1},y\rangle,~{}\text{such that
}C+D_{y}\succeq 0$ (14)

We know that $C+D_{y^{*}}=\frac{S^{*}}{1+\beta}\succeq 0$ from Observation
A.2. Therefore, $y^{*}$ is a feasible solution to the dual problem (14),
revealing that $f^{*}\geq-\langle\mathds{1},y^{*}\rangle$. This fact allows us
to derive the inequality given by $ii)$. Then, if we define
$\zeta=\sigma^{2}_{\min}(R)\sigma_{nnz}(S^{*})$ and
$\gamma=2\sigma^{2}_{\max}(R)$ with $\sigma_{\min}(R)=1/(w\odot y)_{\max}$ and
$\sigma_{\max}(R)=1/(w\odot y)_{\min}$, we arrive at the final result.

$\displaystyle\|V-\hat{V}\|^{2}_{F}\geq\left(\zeta-\gamma\|\Delta\|_{F}\right)(f(V)-f^{\star})$

It should be noted that both $\zeta$ and $\gamma$ are strictly positive
constants because $R$ is defined as an invertible matrix and, therefore, has
all non-zero eigenvalues. ∎

### Proof of Lemma 4.3

###### Proof.

Recall that
$\Delta_{i}=(w_{i}+\beta)y_{i}-(w_{i}^{\star}+\beta)y_{i}^{\star}$. We can
unroll the value of $\Delta_{i}\hat{v}_{i}$ into the following expression.

$\displaystyle\big{(}(w_{i}y_{i}-w_{i}^{*}y_{i}^{*})+\beta(y_{i}-y_{i}^{*})\big{)}\,\hat{v_{i}}$
(15)

Then, we derive the following equality from (9).

$\begin{split}w_{i}y_{i}\hat{v_{i}}&=-\left(1+\beta\right)g_{i}-\beta
y_{i}v_{i}\\\ &=-\left(1+\beta\right)Z_{i}c_{i}-\beta y_{i}v_{i}\end{split}$
(16)

where we denote $Z_{i}=Z_{i}(V)$ for brevity and $c_{i}$ as the $i$th column
of $C$. We denote
$s_{i}^{*}=\left(1+\beta\right)c_{i}+w_{i}^{*}y_{i}^{*}e_{i}+\beta
y_{i}^{*}e_{i}$. By adding (15) to both sides of (16), we arrive at the
following, expanded version of the update rule in (9).

$\begin{split}w_{i}y_{i}\hat{v_{i}}-w_{i}^{*}y_{i}^{*}\hat{v_{i}}+\beta
y_{i}\hat{v_{i}}-\beta
y_{i}^{*}\hat{v_{i}}&=-\left(1+\beta\right)Z_{i}c_{i}-\beta
y_{i}v_{i}-w_{i}^{*}y_{i}^{*}\hat{v_{i}}+\beta y_{i}\hat{v_{i}}-\beta
y_{i}^{*}\hat{v_{i}}\\\ &\stackrel{{\scriptstyle
i)}}{{=}}-Z_{i}s_{i}^{*}+\big{(}w_{i}^{*}y_{i}^{*}-\beta(y_{i}-y_{i}^{*})\big{)}(v_{i}-\hat{v_{i}}).\end{split}$
(17)

where $i)$ follows from the definition of $s_{i}^{\star}$ and the fact that
$Z_{i}e_{i}=v_{i}$. By combining Assumption 4.1 and Lemma 4.1, it is known
that

$f(V)-f^{*}\geq\frac{1-\beta}{1+\beta}\,\delta\,\|v_{i}-\hat{v_{i}}\|_{2}^{2}.$
(18)

Besides,

$\begin{split}f(V)-f^{*}&\geq f(Z_{i})-f^{*}\\\
&=\operatorname*{\text{Tr}}(Z_{i}^{\top}Z_{i}S^{*})\\\
&\geq\sigma^{-1}_{\max}(S^{*})\|Z_{i}S^{*}\|_{F}^{2}\\\
&\geq\sigma^{-1}_{\max}(S^{*})\|Z_{i}s_{i}^{*}\|_{2}^{2}\end{split}$ (19)

Therefore, based on (18) and (19), it is known that
$\|Z_{i}s_{i}^{\star}\|_{2}^{2}$ and $\|v_{i}-\hat{v}_{i}\|^{2}_{2}$ are both
upper bounded by a positive factor of $f(V)-f^{\star}$. We take the inner
product of (15) with itself to yield the following.

$\displaystyle\Delta_{i}^{2}$ $\displaystyle\stackrel{{\scriptstyle
i)}}{{=}}\big{(}(w_{i}y_{i}-w_{i}^{*}y_{i}^{*})+\beta(y_{i}-y_{i}^{*})\big{)}^{2}$
$\displaystyle\stackrel{{\scriptstyle
ii)}}{{\leq}}2\|Z_{i}s_{i}^{*}\|^{2}+2(y^{*}_{i}-\beta
y_{i})^{2}\|v_{i}-\hat{v_{i}}\|^{2}$
$\displaystyle=2\|Z_{i}s_{i}^{*}\|^{2}+2({y^{*}_{i}}^{2}-2\beta
y_{i}^{*}y_{i}+{y_{i}}^{2})\,\|v_{i}-\hat{v_{i}}\|^{2}$

where $i)$ follows from the fact that $\hat{v}_{i}$ is a unit vector and $ii)$
follows from combining observation A.1 with (17). From here, one can notice
that

$\|(w+\beta\mathds{1})\odot y-(w^{*}+\beta\mathds{1})\odot
y^{*}\|_{2}^{2}=\|\Delta\|_{F}^{2}=\sum\limits_{i=1}^{n}\Delta_{i}^{2}$

. We combine this expression with the upper bounds for
$\|Z_{i}s_{i}^{\star}\|^{2}_{2}$ and $\|v_{i}-\hat{v}_{i}\|^{2}_{2}$ derived
in (18) and (19) to arrive at the final result, where $\tau$ is defined as
some positive constant.

$\displaystyle\|(w+\beta\mathds{1})\odot y-(w^{*}+\beta\mathds{1})\odot
y^{*}\|_{2}^{2}\leq\tau(f(V)-f^{*})$

∎

To finish we show that linear convergence is an immediate consequence of the
previous lemmas.

### Proof of Theorem 4.1

###### Proof.

Based on Lemma 4.3, a neighborhood can be selected around the optimum such
that the value of $f(V)-f^{\star}$ is sufficiently small. In particular, we
define this neighborhood through the selection of a positive constant $\kappa$
such that the following inequality holds, where $\tau$ is defined in Lemma
4.3.

$\displaystyle\left(\frac{\zeta-\kappa}{\gamma}\right)^{2}\geq\tau\left(f(V)-f^{\star}\right)$
(20)

Within this expression, $\zeta$ and $\gamma$ are both defined in Lemma 4.2.
Because $\zeta$ and $\gamma$ are both strictly positive (i.e., see Lemma 4.2),
$\left(\frac{\zeta-\kappa}{\gamma}\right)$ is known to be strictly positive so
long as $\kappa\in[0,\zeta)$. Therefore, there always exists a value of
$\kappa$ such that (20) will be true within a sufficiently small neighborhood
around the optimum. From here, we combine (20) with the inequality from Lemma
4.3 to yield the following.

$\left(\frac{\zeta-\kappa}{\gamma}\right)^{2}\geq\tau(f(V)-f^{*})\geq\|\Delta\|_{F}^{2},$

Then, because $\left(\frac{\zeta-\kappa}{\gamma}\right)$ is strictly positive,
we can manipulate the above expression to yield the following.

$\zeta-\gamma\|\Delta\|_{F}\geq\kappa.$ (21)

The above inequality holds for all successive iterations of Algorithm 2.

Then, based on the combination of Lemma 4.1 and assumption 4.1,
$f(V)-f(\hat{V})\geq\delta\frac{1-\beta}{1+\beta}\|V-\hat{V}\|_{F}^{2}$, and
thus

$\begin{split}f(V)-f(\hat{V})&\geq\delta\frac{1-\beta}{1+\beta}\|V-\hat{V}\|_{F}^{2}\\\
&\stackrel{{\scriptstyle
i)}}{{\geq}}\delta\frac{1-\beta}{1+\beta}\kappa(f(V)-f^{*})\end{split}$

where $i)$ is derived by combining (21) and Lemma 4.2. This expression in turn
implies the following.

$\left(1-\delta\kappa\frac{1-\beta}{1+\beta}\right)(f(V)-f^{*})\geq
f(\hat{V})-f^{*}.$

thus giving us the desired linear convergence. By substituting
$\rho=\delta\kappa\frac{1-\beta}{1+\beta}$, we arrive at the final result. ∎

##  Appendix B Convergence to first-order critical point

Mixing method++ not only has local linear convergence, but it always converges
to a first-order critical point, per following Theorem:

######  Theorem B.1.

Let $V^{\ell}=Z_{n+1}^{\ell}(V)$, see definition (10), for all
$\ell\in\mathbb{N}$; i.e., the collection of points generated by Mixing
Method++ after finishing the inner loop. Then, under the assumption 4.1,
$\\{V^{\ell}\\}_{\ell}$ converges to a limit point $\bar{V}$ and the
riemannian gradient at $\bar{V}$ is zero (Mixing Method++ converges to a
first-order critical point).

###### Proof.

For $\beta$ in $[0,1)$ we know by lemma 4.1 that the objective function is
decreasing. Besides, $(S^{k-1})^{n}$ is compact, indeed limit point compact,
so there exists a limit point $\bar{V}$ such that
$f(\bar{V})=\lim_{k}f(V^{k})$ by continuity of $f$. It is clear that $\bar{V}$
is a fixed point in relation (12), so

$\bar{V}C=-\frac{\bar{V}D_{(\bar{w}+\beta\mathds{1})\odot\bar{y}}}{1+\beta}.$
(22)

Since the constraint set is a product of spheres, its corresponding tangent
space (henceforth denoted by the letter $\mathrm{T}$) at the point $\bar{V}$
is given by the product of tangent spaces of $S^{k-1}$; namely

$\mathrm{T}_{\bar{V}}\left(S^{k-1}\right)^{n}=\prod_{j=1}^{n}\mathrm{T}_{\bar{v_{j}}}\,\left(S^{k-1}\right),$
(23)

where $\bar{v}_{j}$ is the $j$th column of $\bar{V}$.

The above result gives us a way to characterize the tangent space. To this
end, it is actually enough to work with the tangent space of $S^{k-1}$.

By definition, for any $v_{j}$, the vector $\dot{v}_{j}\in\mathbb{R}^{k}$ is
in the tangent space of $S^{k-1}$ if and only if there exists a curve
$\varphi:I\to S^{k-1}$ such that $\varphi(0)=v_{j}$ and
$\dot{\varphi}(0)=\dot{v}_{j}$, where $\dot{\varphi}$ is the derivative of
$\varphi$. Then $\varphi(t)\in S^{k-1}$ if satisfies
$\langle\varphi(t),\varphi(t)\rangle=1$ for all $t\in I$. Differentiating on
both sides leads us to
$0=\langle\dot{\varphi}(t),\varphi(t)\rangle+\langle\varphi(t),\dot{\varphi}(t)\rangle$.
Evaluating at $t=0$, we get

$2\langle\dot{v}_{j},v_{j}\rangle=0,$

thus
$\\{\dot{v}_{j}\in\mathbb{R}^{k}:\langle\dot{v}_{j},v_{j}\rangle=0\\}\subseteq\mathrm{T}_{v_{j}}S^{k-1}$,
since both subspaces are of the same dimension,
$\\{\dot{v}_{j}\in\mathbb{R}^{k}:\langle\dot{v}_{j},v_{j}\rangle=0\\}=\mathrm{T}_{v_{j}}S^{k-1}$.

Now by (23) and the previous analysis, we finally arrive to the
characterization of the tangent space; namely

$\mathrm{T}_{\bar{V}}\left(S^{k-1}\right)^{n}=\left\\{\dot{\bar{V}}\in\mathbb{R}^{k\times
n}:\langle\dot{\bar{v}}_{j},\bar{v}_{j}\rangle=0\hskip 5.0pt\forall
j\in[n]\,\right\\},$

and $\dot{\bar{v}}_{j}$ ($v_{j}$) corresponds to the $j$th column of
$\dot{\bar{V}}$ ($\bar{V}$), respectively.

Let $\mathrm{P}^{\perp}:\mathbb{R}^{k\times
n}\to\mathrm{T}_{\bar{V}}\left(S^{k-1}\right)^{n}$ be the projector operator
from the euclidean space to the tangent space at $\bar{V}$, defined for any
$W\in\mathbb{R}^{k\times n}$ as follows,

$W\xrightarrow{\mathrm{P}^{\perp}}\begin{pmatrix}w_{1}-\langle\bar{v}_{1},w_{1}\rangle\bar{v}_{1}&\ldots&w_{n}-\langle\bar{v}_{n},w_{n}\rangle\bar{v}_{n}\end{pmatrix}.$

Let the riemannian gradient of $f$ denoted by $\mathrm{grad}f(\cdot)$, and the
gradient of $f$ defined in the entire euclidean domain as $\nabla f(\cdot)$,
so using tools from matrix manifold the riemannian gradient is given by

$\mathrm{grad}f(\bar{V})=\mathrm{P}^{\perp}(\nabla f(\bar{V})).$

It is easy to see that $\nabla f(\bar{V})=2\bar{V}C$, and by relation (22),
this yields us to the equivalent version $\nabla
f(\bar{V})=-\frac{2}{1+\beta}\bar{V}D_{(\bar{w}+\beta\mathds{1})\odot\bar{y}}$.

To conclude, let us fix a coordinate, say $j\in[n]$, so $\left(\nabla
f(\bar{V})\right)_{j}=-\frac{2}{1+\beta}(w_{j}y_{j}+\beta y_{j})\bar{v}_{j}$
and by definition of the projection map we get

$\begin{split}\left(\mathrm{P}^{\perp}(\nabla
f(\bar{V})\right)_{j}=\frac{2}{1+\beta}(w_{j}y_{j}+\beta
y_{j})(-\bar{v}_{j}+\bar{v}_{j})=0.\end{split}$

So for all the coordinates $\left(\mathrm{P}^{\perp}(\nabla
f(\bar{V})\right)_{j}=0$, and thus $\mathrm{P}^{\perp}(\nabla
f(\bar{V})=\mathrm{grad}f(\bar{V})=0$. ∎

##  Appendix C Empirical Analysis of Assumption 4.1

Figure 8: The minimum norm of the columns of $V$ throughout iterations of
Mixing Method and Mixing Method++ on the MaxCut and MIMO objectives.

The theoretical analysis of both Mixing Method and Mixing Method++ relies on
Assumption 4.1, which is a non-degeneracy assumption on the norms of the
columns of $V$. This assumption states that norms of the columns of $V$ being
updated should never be zero (i.e., $\|\sum_{j\neq
i}c_{ij}v_{j}\|_{2}>0~{}~{}\forall i$). To determine whether this assumption
is valid in practice, we run both Mixing Method and Mixing Method++ solvers on
random instances of the MaxCut and MIMO objectives. For every iteration, we
track the smallest value of $\|\sum_{j\neq i}c_{ij}v_{j}\|_{2}$ that is
observed during each algorithm’s inner loop. We then plot these minimum norms
throughout all iterations. The results of these experiments are shown in
Figure 8. As can be seen, Assumption 4.1 holds in practice on MaxCut and MIMO
objectives for both Mixing Method and Mixing Method++.

##  Appendix D Supplemental Results

###  D.1 MaxCut

Figure 9: CPU time vs Size on MaxCut. We run all SDP solvers with time limit
24h ($8.64\times 10^{4}$s). Mixing Method++ solves most cases and is the
fastest solver.

Figure 10: Objective residual vs size on MaxCut. Mixing Method++ shows
competitive performance. As an approximate method, CGAL has large objective
residual.

Here, we demonstrate the performance properties of Mixing Method++ on
differently-sized instances of the MaxCut objective. In Figures 9 and 10, we
report the CPU time and objective residual of solving MaxCut instances as a
function of the number of nodes in the graph $n$. As can be seen, Mixing
Method++ decreases required CPU time to solve MaxCut problems of all sizes in
comparison to numerous other SDP solvers. In fact, Mixing Method++ is able to
solve nearly all MaxCut instances that were tested within the imposed 24 hour
time limit. Furthermore, the objective residual obtained by Mixing Method++
remains competitive throughout for all problem sizes and often surpasses the
performance of the other solvers, thus demonstrating that Mixing Method++
yields improved solution quality in addition to providing acceleration.

###  D.2 MaxSAT

Figure 11: MaxSAT: cost of solution (time limit $=300$s). Though Loandra, the
champion of MAXSAT 2019 is the best solver for large industrial instances,
Mixing Method++ improves the quality of solution compared to Mixing Method.

Here, we consider MaxSAT problems within the industrial category of the
2016-2019 MaxSAT competition. In particular, we report the cost of solution
(i.e., the number of unsatisfied clauses) achieved by each solver for
different sizes of problem instances. As can be seen in Figure 11, Loandra
(i.e., the best solver for large industrial problem instances) achieves the
best cost of solution across problem sizes. However, Mixing Method++
consistently improves upon the performance of the Mixing Method in terms of
the number of unsatisfied clauses.