apxtheorem[theorem]Theorem apxlemma[theorem]Lemma apxfact[theorem]Fact

#  A Faster Interior-Point Method for Sum-of-Squares Optimization

Shunhua Jiang
Columbia University
<EMAIL_ADDRESS>Supported by NSF CAREER award CCF-1844887. Bento Natura
London School of Economics
<EMAIL_ADDRESS>This project has received funding from the European
Research Council (ERC) under the European Union’s Horizon 2020 research and
innovation programme (grant agreement no. 757481–ScaleOpt). Omri Weinstein
The Hebrew University and Columbia University
<EMAIL_ADDRESS>Supported by NSF CAREER award CCF-1844887 and ISF grant
#3011005535.

###### Abstract

We present a faster interior-point method for optimizing sum-of-squares (SOS)
polynomials, which are a central tool in polynomial optimization and capture
convex programming in the Lasserre hierarchy. Let $p=\sum_{i}q^{2}_{i}$ be an
$n$-variate SOS polynomial of degree $2d$. Denoting by $L:=\binom{n+d}{d}$ and
$U:=\binom{n+2d}{2d}$ the dimensions of the vector spaces in which $q_{i}$’s
and $p$ live respectively, our algorithm runs in time $\tilde{O}(LU^{1.87})$.
This is polynomially faster than state-of-art SOS and semidefinite programming
solvers, which achieve runtime
$\tilde{O}(L^{0.5}\min\\{U^{2.37},L^{4.24}\\})$.

The centerpiece of our algorithm is a dynamic data structure for maintaining
the inverse of the Hessian of the SOS barrier function under the _polynomial
interpolant basis_ , which efficiently extends to multivariate SOS
optimization, and requires maintaining spectral approximations to low-rank
perturbations of _elementwise (Hadamard) products_. This is the main challenge
and departure from recent IPM breakthroughs using inverse-maintenance, where
low-rank updates to the slack matrix readily imply the same for the Hessian
matrix.

##  1 Initialization

There exist standard techniques to transform a convex program to a form that
has an easily obtainable strictly feasible point, see e.g. [ytm94]. We follow
the initialization procedure presented by [cls19] and [jiang2020faster] and
adapt to SOS optimization. Similar initialization lemma exists for WSOS
optimization.

Let the matrix $P\in\mathbb{R}^{U\times L}$ and the operator
$\Lambda:\mathbb{R}^{U}\to\mathbb{R}^{L\times L}$ that
$\Lambda(s)=P^{\top}\operatorname{diag}(s)P$ be defined as in the interpolant
basis paragraph of Section 6.

######  Lemma 1.1 (Initialization).

Given an instance of $\eqref{SOS_program_primal_dual}$ that fulfills Slater’s
condition, and let $R$ be an upper bound on the $\ell_{1}$-norm of the primal
feasible solutions, i.e. all primal feasible $x$ of (SOS) fulfill
$\|x\|_{1}\leq R$, and let $\delta\in(0,1)$. We define

$\overline{A}=\begin{bmatrix}A&0&\frac{1}{R}b-A\overline{g}^{0}\\\
1_{U}^{\top}&1&0\end{bmatrix}\in\mathbb{R}^{(m+1)\times(U+2)},\;\overline{b}=\begin{bmatrix}\frac{1}{R}b\\\
1+\langle
1_{U},\overline{g}^{0}\rangle\end{bmatrix}\in\mathbb{R}^{m+1},\text{and
}\overline{c}=\begin{bmatrix}\frac{\delta}{\|c\|_{\infty}}c\\\ 0\\\
1\end{bmatrix}\in\mathbb{R}^{U+2},$

and let

$\overline{x}^{0}=\begin{bmatrix}\overline{g}^{0}\\\ 1\\\
1\end{bmatrix}\in\mathbb{R}^{U+2},\;\overline{y}^{0}=\begin{bmatrix}0_{m}\\\
-1\end{bmatrix}\in\mathbb{R}^{m+1},\text{and
}\overline{s}^{0}=\begin{bmatrix}1_{U}+\frac{\delta}{\|c\|_{\infty}}c\\\ 1\\\
1\end{bmatrix}\in\mathbb{R}^{U+2},$

where
$\overline{g}^{0}=g_{\Sigma^{*}}(\overline{s}^{0}_{[:U]})\in\mathbb{R}^{U}$
for the gradient function
$g_{\Sigma^{*}}(s):=\operatorname{diag}(P(P^{\top}\operatorname{diag}(s)P)^{-1}P^{\top})$
that maps from $\mathbb{R}^{U}$ to $\mathbb{R}^{U}$. This defines the
auxiliary primal-dual system

$\begin{aligned} \min\;&\langle\overline{c},\overline{x}\rangle\quad\\\
\overline{A}\overline{x}&=\overline{b}\\\
\overline{x}&\in\Sigma_{n,2d}\times\mathbb{R}_{\geq 0}^{2}\,,\\\
\end{aligned}\quad\quad\quad\begin{aligned}
\max\;&\langle\overline{y},\overline{b}\rangle\\\
\overline{A}^{\top}\overline{y}+\overline{s}&=\overline{c}\\\
\overline{s}&\in\Sigma_{n,2d}^{*}\times\mathbb{R}_{\geq 0}^{2}\,.\\\
\end{aligned}$ (Aux-SOS)

Then $(\overline{x}^{0},\overline{y}^{0},\overline{s}^{0})$ are feasible to
the auxiliary system (Aux-SOS).

Further, under the canonical barrier (we use $\overline{a}_{i}$ to denote the
$i$-th column of $\overline{A}$):

$\overline{F}_{\eta}(\overline{y})=-\eta\langle\overline{y},\overline{b}\rangle-\log\det\Big{(}\Lambda\big{(}(\overline{c}-\overline{A}^{\top}\overline{y})_{[:U]}\big{)}\Big{)}-\log(\overline{c}_{U+1}-\langle\overline{a}_{U+1},\overline{y}\rangle)-\log(\overline{c}_{U+2}-\langle\overline{a}_{U+2},\overline{y}\rangle),$
(1)

we have that
$\|\overline{g}_{\eta^{0}}(\overline{y}^{0})\|_{\overline{H}(\overline{y}^{0})^{-1}}=0$
for $\eta^{0}=1$.

Further, for any solution $(\overline{x},\overline{y},\overline{s})$ to (Aux-
SOS) with duality gap $\leq\delta^{2}$, its restriction
$\widehat{x}:=\overline{x}_{[:U]}$ fulfills

$\displaystyle\langle c,\widehat{x}\rangle$
$\displaystyle\leq\min_{Ax=b,x\in\Sigma_{n,2d}}\langle c,x\rangle+\delta\cdot
R\|c\|_{\infty},$ $\displaystyle\|A\widehat{x}-b\|_{1}$ $\displaystyle\leq
8\delta L\cdot(LR\|A\|_{\infty}+\|b\|_{1}),$ $\displaystyle\widehat{x}$
$\displaystyle\in\Sigma_{n,2d}.$

###### Proof.

Let $\operatorname{proj}$ denote the orthogonal projection matrix that maps a
matrix onto its image space. Then note that

$\displaystyle\|\overline{g}^{0}\|_{1}$
$\displaystyle=\Big{\|}\operatorname{diag}\Big{(}P\Big{(}P^{\top}\operatorname{diag}(1_{U}+\frac{\delta}{\|c\|_{\infty}}c)P\Big{)}^{-1}P^{\top}\Big{)}\Big{\|}_{1}$
(2)
$\displaystyle=\operatorname{tr}\Big{(}\operatorname{diag}\big{(}1+\frac{\delta}{\|c\|_{\infty}}c\big{)}^{-1}\cdot\operatorname{proj}\big{(}\operatorname{diag}(1_{U}+\frac{\delta}{\|c\|_{\infty}}c)^{1/2}\cdot
P\big{)}\Big{)}$ $\displaystyle=(1\pm\delta)L,$

where the last inequality follows as
$\operatorname{tr}(\operatorname{proj}(P))=\operatorname{rk}(P)=L$.

Straightforward calculations show that
$\overline{A}\overline{x}^{0}=\overline{b}$ and
$\overline{A}^{\top}\overline{y}^{0}+\overline{s}^{0}=\overline{c}$. Further,
note that $\overline{s}^{0}>0$ as $\delta<1$ and therefore
$\overline{s}^{0}\in\mathbb{R}_{+}^{U+2}\subset\Sigma_{n,2d}^{*}\times\mathbb{R}_{+}^{2}$.
The containment $\mathbb{R}_{+}^{U}\subset\mathbb{R}_{\geq
0}^{U}\subset\Sigma_{n,2d}^{*}$ is clear by the characterisation
$\Sigma_{n,2d}^{*}=\\{s:\Lambda(s)=P^{\top}\operatorname{diag}(s)P\succeq
0\\}$.

Further
$\overline{x}^{0}_{[:U]}=g_{\Sigma^{*}}(s^{0}_{[:U])})=\operatorname{diag}(P(P^{\top}\operatorname{diag}(\overline{s}^{0}_{[:U]})P)^{-1}P^{\top})\in\Sigma_{n,2d}$
is e.g. shown in Proposition 3.3.3 in
[renegarMathematicalViewInteriorPoint2001]. We make it explicit by showing
that $\langle\overline{x}^{0}_{[:U]},\hat{s}\rangle\geq 0$ for any
$\hat{s}\in\Sigma_{n,2d}^{*}$. Indeed,

$\displaystyle\langle\overline{x}^{0}_{[:U]},\hat{s}\rangle$
$\displaystyle=\langle\operatorname{diag}(P(P^{\top}\operatorname{diag}(\overline{s}^{0}_{[:U]})P)^{-1}P^{\top}),\hat{s}\rangle$
$\displaystyle=\operatorname{tr}\big{(}(P(P^{\top}\operatorname{diag}(\overline{s}^{0}_{[:U]})P)^{-1}P^{\top})\cdot\operatorname{diag}(\hat{s})\big{)}$
$\displaystyle=\operatorname{tr}{\big{(}(P^{\top}\operatorname{diag}(\overline{s}^{0}_{[:U]})P)^{-1}(P^{\top}\operatorname{diag}(\hat{s})P)\big{)}}$
$\displaystyle\geq 0,$

where the last inequality follows as the argument in $\operatorname{tr}$ is a
product of positive semi-definite matrices.

Define $\overline{s}=\overline{c}-\overline{A}^{\top}\overline{y}$. The
corresponding barrier (as defined in (1)) is

$\displaystyle\overline{F}_{\eta}(\overline{y})$
$\displaystyle=-\eta\langle\overline{y},\overline{b}\rangle-\log\det\Big{(}\Lambda\big{(}(\overline{c}-\overline{A}^{\top}\overline{y})_{[:U]}\big{)}\Big{)}-\log(\overline{c}_{U+1}-\langle\overline{a}_{U+1},\overline{y}\rangle)-\log(\overline{c}_{U+2}-\langle\overline{a}_{U+2},\overline{y}\rangle)$
$\displaystyle=-\eta\langle\overline{y},\overline{b}\rangle-\log\det\Big{(}\Lambda\big{(}\overline{s}_{[:U]}\big{)}\Big{)}-\log(\overline{c}_{U+1}-y_{m+1})-\log(\overline{c}_{U+2}-\langle\frac{1}{R}b-A\overline{g}^{0},\overline{y}_{[:m]}\rangle),$

and gradient of the system is

$\displaystyle\overline{g}_{\eta}(\overline{y})$
$\displaystyle=\begin{bmatrix}-\frac{\eta}{R}b+Ag_{\Sigma^{*}}(\overline{s}_{[:U]})+\frac{1}{\overline{s}_{U+2}}(\frac{1}{R}b-A\overline{g}^{0})\\\
-\eta(1+\langle 1_{U},\overline{g}^{0}\rangle)+\langle
1_{U},g_{\Sigma^{*}}(\overline{s}_{[:U]})\rangle+\frac{1}{\overline{s}_{U+1}}\end{bmatrix}.$

Simple algebraic calculations now show that for $\eta^{0}=1$ we have
$\overline{g}_{\eta^{0}}(\overline{y}^{0})=0_{m+1}$, and so in particular
$\|\overline{g}_{\eta^{0}}(\overline{y}^{0})\|_{\overline{H}(\overline{y}^{0})^{-1}}\leq\epsilon_{N}$
as necessary in Section 9.1. It remains to show that near-optimal solutions to
the modified problem correspond to near-optimal solutions to the original
problem.

Let $\mathrm{OPT}$ resp. $\overline{\mathrm{OPT}}$ be the objective values of
the original resp. modified program:

$\mathrm{OPT}:=\min_{Ax=b,x\in\Sigma_{n,2d}}\langle
c,x\rangle,\qquad\overline{\mathrm{OPT}}:=\min_{\overline{A}x=\overline{b},\overline{x}\in\Sigma_{n,2d}\times\mathbb{R}^{2}}\langle\overline{c},\overline{x}\rangle.$

Given any optimal $x^{*}\in\mathbb{R}^{U}$ of the original primal of (SOS),
consider the following $\overline{x}^{*}\in\mathbb{R}^{U+2}$ fulfilling
$\overline{A}\overline{x}^{*}=\overline{b},\overline{x}^{*}\in\Sigma_{n,2d}\times\mathbb{R}_{\geq
0}^{2}$:

$\overline{x}^{*}=\begin{bmatrix}\frac{1}{R}x^{*}\\\ \langle
1_{U},\overline{g}^{0}\rangle+1-\frac{1}{R}\langle 1_{U},x^{*}\rangle\\\
0\end{bmatrix}.$

Note that indeed $\overline{x}^{*}_{U+1}\in\mathbb{R}_{\geq 0}$ as by (2) we
have $x^{*}_{U+1}=\langle 1_{U},\overline{g}^{0}\rangle+1-\frac{1}{R}\langle
1_{U},x^{*}\rangle\geq(1-\delta)L+1-1\geq 0$ for $\delta<1$ and by choice of
$R$. Therefore

$\overline{\mathrm{OPT}}\leq\langle\overline{c},\overline{x}^{*}\rangle=\begin{bmatrix}\frac{\delta}{\|c\|_{\infty}}c^{\top}&0&1\end{bmatrix}\cdot\begin{bmatrix}\frac{1}{R}x^{*}\\\
\langle 1_{U},\overline{g}^{0}\rangle+1-\frac{1}{R}\langle
1_{U},x^{*}\rangle\\\ 0\end{bmatrix}=\frac{\delta\cdot\langle
c,x^{*}\rangle}{R\|c\|_{\infty}}=\frac{\delta\cdot\mathrm{OPT}}{R\|c\|_{\infty}}.$
(3)

Now given any feasible primal-dual solution
$(\overline{x},\overline{y},\overline{s})$ to (Aux-SOS) with duality gap
$\leq\delta^{2}$ we have

$\frac{\delta}{\|c\|_{\infty}}\langle
c,\overline{x}_{[:U]}\rangle\leq\frac{\delta}{\|c\|_{\infty}}\langle
c,\overline{x}_{[:U]}\rangle+\overline{x}_{U+2}=\langle\overline{c},\overline{x}\rangle\leq\overline{\mathrm{OPT}}+\delta^{2}\leq\frac{\delta}{R\|c\|_{\infty}}\mathrm{OPT}+\delta^{2},$
(4)

where the first inequality uses $\overline{x}_{U+2}\in\mathbb{R}_{\geq 0}$ and
the last inequality follows from (3). So, for
$\widehat{x}:=R\cdot\overline{x}_{[:U]}$, using (4) we have that

$\langle c,\widehat{x}\rangle=R\langle
c,\overline{x}_{[:U]}\rangle=\frac{R\|c\|_{\infty}}{\delta}\frac{\delta}{\|c\|_{\infty}}\langle
c,\overline{x}_{[:U]}\rangle\leq\frac{R\|c\|_{\infty}}{\delta}\Bigl{(}\frac{\delta}{R\|c\|_{\infty}}\mathrm{OPT}+\delta^{2}\Bigr{)}=\mathrm{OPT}+\delta
R\|c\|_{\infty}.$

Recall that $\mathbb{R}_{\geq 0}^{U}\subset\Sigma_{n,2d}^{*}$ and therefore
$\Sigma_{n,2d}\subset\mathbb{R}_{\geq 0}^{U}$. In particular
$\overline{x}_{[:U]}\geq 0$ and so using the equality $\langle
1_{U},\overline{x}_{[:U]}\rangle+\overline{x}_{U+1}=1+\langle
1_{U},\overline{g}^{0}\rangle$ (which follows from
$\overline{A}\overline{x}=\overline{b}$ since
$(\overline{x},\overline{y},\overline{s})$ is feasible) and
$\overline{x}_{U+1}\geq 0$ we can further bound

$\|\overline{x}_{[:U]}\|_{1}=\langle 1_{U},\overline{x}_{[:U]}\rangle\leq
1+\langle 1_{U},\overline{g}^{0}\rangle\leq 1+(1+\delta)L\leq 2L\,,$

where the penultimate inequality used (2). So we can further bound

$\frac{1}{\|c\|_{\infty}}\langle
c,\overline{x}_{[:U]}\rangle\geq-\frac{1}{\|c\|_{\infty}}\|c\|_{\infty}\|\overline{x}_{[:U}\|_{1}\geq-2L.$
(5)

Therefore using (4) and (5) we get that

$\overline{x}_{U+2}\leq\frac{\delta}{R\|c\|_{\infty}}\mathrm{OPT}+\delta^{2}-\frac{\delta}{\|c\|_{\infty}}\langle
c,\overline{x}_{[:U]}\rangle\leq\frac{\delta}{R\|c\|_{\infty}}\mathrm{OPT}+\delta^{2}+2\delta
L\leq\delta+\delta^{2}+2\delta L\leq 4\delta L.$ (6)

where the penultimate inequality used $\mathrm{OPT}\leq R\|c\|_{\infty}$ and
the last inequality uses $\delta\leq 1$.

Feasibility of $\overline{x}$ to $\overline{A}\overline{x}=\overline{b}$ shows
that
$A\overline{x}_{[:U]}+(\frac{1}{R}b-A\overline{g}^{0})\overline{x}_{U+2}=\frac{1}{R}b$
and we therefore have

$\displaystyle\|A\widehat{x}-b\|_{1}$
$\displaystyle=\|RA\overline{x}_{[:U]}-b\|_{1}=\|(b-RA\overline{g}^{0})\overline{x}_{U+2}\|_{1}\leq
4\delta L\cdot(R\|A\|_{\infty}\|\overline{g}^{0}\|_{1}+\|b\|_{1})$
$\displaystyle\leq 8\delta L\cdot(LR\|A\|_{\infty}+\|b\|_{1}).$

where the penultimate inequality used (6), and the last inequality follows
from (2). ∎

It is worth noting that the system
$\overline{A}\overline{x}=\overline{b},\overline{x}\in\Sigma_{n,2d}\times\mathbb{R}_{\geq
0}^{2}$ is not an instance of $\eqref{SOS_program_primal_dual}$ anymore. This
is in contrast for similar initialization techniques for Linear Programming
[cls19] or semi-definite Programming [jiang2020faster]. Nonetheless it is not
hard to see that optimization over product cones as in this case
$\Sigma_{n,2d}\times\mathbb{R}_{\geq 0}^{2}$ can be done efficiently if the
optimization over the factors can be done efficiently, see e.g. the
application in [lsz19]. Alternatively note that the dual cone
$(\Sigma_{n,2d}\times\mathbb{R}_{\geq
0}^{2})^{*}=\Sigma_{n,2d}^{*}\times\mathbb{R}_{\geq 0}^{2}$. For
$s\in\Sigma_{n,2d}^{*}$ we used the barrier
$F(s)=-\log(\det(P^{\top}\operatorname{diag}(s)P))$ throughout the paper. This
barrier could easily extend via setting

$\overline{P}:=\begin{bmatrix}P&0&0\\\ 0&1&0\\\ 0&0&1\\\ \end{bmatrix}$

and setting the barrier for
$\overline{s}\in\Sigma_{n,2d}^{*}\times\mathbb{R}_{\geq 0}^{2}$ as

$\overline{F}(\overline{s})=-\log\det(\overline{P}^{\top}\operatorname{diag}(\overline{s})\overline{P})=-\log(\det(P^{\top}\operatorname{diag}(\overline{s}_{[:U]})P))-\log\overline{s}_{U+1}-\log\overline{s}_{U+2},$

which recovers the standard barrier for the product cone. Our algorithm
Algorithm 1 can now equally be run with $\overline{F}$ instead of $F$, where
then $P$ is replaced by $\overline{P}$.

##  2 Discussion on representation

In this paper we consider the dual formulation of (SOS), where
$A\in\mathbb{R}^{m\times U}$, $b\in\mathbb{R}^{m}$, $c\in\mathbb{R}^{U}$:

$\displaystyle\max$ $\displaystyle\langle y,b\rangle$ (7) $\displaystyle
A^{\top}y+s$ $\displaystyle=c$ $\displaystyle s$
$\displaystyle\in\Sigma_{n,2d}^{*}\,.$

Here $m$ is the number of constraints, and w.l.o.g. we assume $m\leq U$.

We remark that [py19] instead consider the following formulation:

$\displaystyle\min\;$ $\displaystyle\langle c,x\rangle$ $\displaystyle Ax$
$\displaystyle=b$ $\displaystyle x$ $\displaystyle\in\Sigma_{n,2d}^{*}\,.$

These two formulations are in fact equivalent, because Eq. (7) is equivalent
to

$\displaystyle\min\;$ $\displaystyle\langle\hat{c},x\rangle\quad$
$\displaystyle\hat{A}x$ $\displaystyle=\hat{b}$ $\displaystyle x$
$\displaystyle\in\Sigma_{n,2d}^{*}\,.$

Here we define $\hat{A}\in\mathbb{R}^{(U-m)\times U}$ to be a matrix that
satisfies $\ker(\hat{A})=\operatorname{Im}(A^{\top})$, i.e., the rows of
$\hat{A}$ are the orthogonal complement of the rows of $A$. Further we define
$\hat{c}\in\mathbb{R}^{U}$ to be any vector that satisfies $A\hat{c}=b$.
Finally, we define $\hat{b}:=\hat{A}c\in\mathbb{R}^{U-m}$.

###  2.1 SOS under the monomial basis

In [py19], the advantages and disadvantages of three different bases (monomial
basis, Chebyshev basis, and interpolant basis) are discussed. Further,
explicit expressions of $\Lambda$ (defined in Theorem 6.4) are given for all
three bases. Apart from the advantage of its computational efficiency, [py19]
chooses interpolant basis in their algorithm also because the interpolant
basis representation is numerically stable, and this is required for practical
algorithms. As we mainly focus on the theoretical running time of SOS
algorithms, in this section we further justify the choice of the interpolant
basis in our algorithm. We want to add to the exposition of [py19] and stress
that the standard monomial basis does not seem suitable even for theoretical
algorithms. It is unclear whether in the monomial basis an amortized runtime
faster than the naive $O(U^{\omega})$ can be achieved. (Similar argument holds
for the Chebyshev basis.)

The Hessian in the monomial basis is given by

$H_{ij}(s)=\sum_{a+b=i,k+\ell=j}[\Lambda(s)^{-1}]_{ak}[\Lambda(s)^{-1}]_{b\ell}.$

It is unclear how low-rank updates techniques could be applied here as low
rank updates to $\Lambda$ do not translate to low-rank updates to the Hessian
$H$. Further, the structure of $\Lambda$ itself in the multivariate case is
far less understood than its counterpart in the interpolant basis, which we
are going to elaborate on in the remainder of this section.

For $n$ variables and degree $d$ the monomial basis elements correspond to the
terms $x_{1}^{\alpha_{1}}\cdot\ldots\cdot x_{n}^{\alpha_{n}}$ for
$\alpha=(\alpha_{1},\ldots,\alpha_{n})\in\mathbb{N}^{n}$, $\|\alpha\|_{1}\leq
d$. When choosing both $\mathbf{p}$ and $\mathbf{q}$ to be the monomial basis,
$\Lambda$ has a special structure. For any basis elements
$p_{1},p_{2}\in\mathbf{p}$ we have that $p_{1}p_{2}$ itself is a monomial in
$\mathcal{V}_{n,2d}$. As such, the coefficients
$\lambda_{ij}\in\mathbb{R}^{U}$ have the special form that $\lambda_{ij}$ is
zero everywhere but in the coordinate that corresponds to the element in
$\mathbf{q}$ equalling $p_{1}p_{2}$. As $\Lambda$ is also uniquely defined by
the images of the elements in $\mathbf{q}$, we can write
$E_{u}\in\mathbb{R}^{L\times L}$ as $E_{u}:=\Lambda(e_{u})$ where $e_{u}$ is
the vector that is zero everywhere but in position $u$ for some $u\in[U]$. Let
$q_{u}\in\mathbf{q}$ be the associated basis polynomial. Then we see that
$(E_{u})_{ij}=1_{[p_{i}p_{j}=q_{u}]}$. It also follows that every matrix
$S\in\operatorname{Im}(\Lambda)$ has at most $U$ different entries and each
entry is uniquely defined by the corresponding basis element in $\mathbf{q}$.
While it is not known how this special structure could be exploited in
general, a speedup is known for the univariate case $n=1$. Here,
$\operatorname{Im}(\Lambda)$ are all Hankel matrices: We have
$\mathbf{q}=\\{1,x,\ldots,x^{2d}\\}$ and $\mathbf{p}=\\{1,x,\ldots,x^{d}\\}$
so for any vector $u\in\mathbb{R}^{2d+1}$ we have that

$\Lambda(v)=\begin{bmatrix}v_{0}&v_{1}&v_{2}&\ldots&\ldots&v_{d}\\\
v_{1}&v_{2}&&&&\vdots\\\ v_{2}&&&&&\vdots\\\ \vdots&&&&&v_{2d-2}\\\
\vdots&&&&v_{2d-2}&v_{2d-1}\\\
v_{d}&\ldots&\ldots&v_{2d-2}&v_{2d-1}&v_{2d}\end{bmatrix}.$

These highly structured matrices are known to be invertible in time
$\tilde{O}(d^{2})$ (see e.g. [Pan2001]).

As mentioned above even the bivariate case becomes far more complicated. Let
$n=2$ and $d=2$, and pick the ordered bases as

$\mathbf{q}=\\{1,x,y,x^{2},xy,y^{2},x^{3},x^{2}y,xy^{2},y^{3},x^{4},xy^{3},x^{2}y^{2},xy^{3},y^{4}\\}\,,\leavevmode\nobreak\
\leavevmode\nobreak\ \mathbf{p}=\\{1,x,y,x^{2},xy,y^{2}\\}\,.$

Then for $v\in\mathbb{R}^{n+2d\choose 2d}=\mathbb{R}^{15}$ we get for
$\mathbf{p}\mathbf{p}^{\top}$ and the corresponding matrix $\Lambda$ that

$\mathbf{p}\mathbf{p}^{\top}=\begin{bmatrix}1&x&y&x^{2}&xy&y^{2}\\\
x&x^{2}&xy&x^{3}&x^{2}y&xy^{2}\\\ y&xy&y^{2}&x^{2}y&xy^{2}&y^{3}\\\
x^{2}&x^{3}&x^{2}y&x^{4}&x^{3}y&x^{2}y^{2}\\\
xy&x^{2}y&xy^{2}&x^{3}y&x^{2}y^{2}&xy^{3}\\\
y^{2}&xy^{2}&y^{3}&x^{2}y^{2}&xy^{3}&y^{4}\end{bmatrix},\quad\Lambda(v)=\begin{bmatrix}v_{0}&v_{1}&v_{2}&v_{3}&v_{4}&v_{5}\\\
v_{1}&v_{3}&v_{4}&v_{6}&v_{7}&v_{8}\\\ v_{2}&v_{4}&v_{5}&v_{7}&v_{8}&v_{9}\\\
v_{3}&v_{6}&v_{7}&v_{10}&v_{11}&v_{7}\\\
v_{4}&v_{7}&v_{8}&v_{11}&v_{12}&v_{13}\\\
v_{5}&v_{8}&v_{9}&v_{7}&v_{13}&v_{14}\\\ \end{bmatrix}.$

While $\Lambda(v)$ still has some structure it is unclear how
$(\Lambda(v))^{-1}$ could be computed more efficiently than in matrix
multiplication time.

##  3 Proof of amortization lemma

We include the proof of Lemma 10.2 for completeness. The main difference
between this proof and that of [hjst21] is that we cut off at $U/L$ instead of
$L$.

Our proof makes use of the following two facts about $\omega$ and $\alpha$
(Lemma A.4 and Lemma A.5 of [cls19]).

######  Fact 3.1 (Relation of $\omega$ and $\alpha$).

$\omega\leq 3-\alpha$.

######  Fact 3.2 (Upper bound of ${\cal T}_{\mathrm{mat}}(n,n,r)$).

For any $r\leq n$, we have that ${\cal T}_{\mathrm{mat}}(n,n,r)\leq
n^{2+o(1)}+r^{\frac{\omega-2}{1-\alpha}}\cdot
n^{2-\frac{\alpha(\omega-2)}{(1-\alpha)}+o(1)}$.

######  Lemma 3.3 (Restatement of Lemma 10.2).

Let $t$ denote the total number of iterations. Let $r_{i}\in[L]$ be the rank
for the $i$-th iteration for $i\in[t]$. Assume $r_{i}$ satisfies the following
condition: for any vector $g\in\mathbb{R}_{+}^{L}$ which is non-increasing, we
have

$\displaystyle\sum_{i=1}^{t}r_{i}\cdot g_{r_{i}}\leq O(t\cdot\|g\|_{2}).$

If the cost in the $i$-th iteration is $O({\cal
T}_{\mathrm{mat}}(U,U,\min\\{Lr_{i},U\\}))$, when $\alpha\geq 5-2\omega$, the
amortized cost per iteration is

$U^{2+o(1)}+U^{\omega-1/2+o(1)}\cdot L^{1/2}.$

###### Proof.

For $r_{i}$ that satisfies $r_{i}\leq U/L$, we have

$\displaystyle{\cal T}_{\mathrm{mat}}(U,U,Lr_{i})\leq$
$\displaystyle\leavevmode\nobreak\
U^{2+o(1)}+(Lr_{i})^{\frac{\omega-2}{1-\alpha}}\cdot
U^{2-\frac{\alpha(\omega-2)}{1-\alpha}+o(1)}$ (8) $\displaystyle=$
$\displaystyle\leavevmode\nobreak\
U^{2+o(1)}+U^{2-\frac{\alpha(\omega-2)}{1-\alpha}+o(1)}\cdot
L^{\frac{\omega-2}{1-\alpha}}\cdot r_{i}^{\frac{\omega-2}{1-\alpha}},$

where the first step follows from Fact 3.2.

Define a sequence $g\in\mathbb{R}_{+}^{L}$ such that for $r\in[L]$,

$\displaystyle
g_{r}=\begin{cases}r^{\frac{\omega-2}{1-\alpha}-1}&\text{if\leavevmode\nobreak\
}r\leq U/L,\\\ (U/L)^{\frac{\omega-2}{1-\alpha}}\cdot
r^{-1}&\text{if\leavevmode\nobreak\ }r>U/L.\end{cases}$

Note that $g$ is non-increasing because $\frac{\omega-2}{1-\alpha}\leq 1$
(Fact 3.1). Then using the condition in the lemma statement, we have

$\displaystyle\sum_{i=1}^{t}\min\\{r_{i}^{\frac{\omega-2}{1-\alpha}},(U/L)^{\frac{\omega-2}{1-\alpha}}\\}=$
$\displaystyle\leavevmode\nobreak\ \sum_{i=1}^{t}r_{i}\cdot g_{r_{i}}$ (9)
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\ t\cdot\|g\|_{2}$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
t\cdot\Big{(}\int_{x=1}^{U/L}x^{\frac{2(\omega-2)}{1-\alpha}-2}\mathrm{d}x+(U/L)^{\frac{2(\omega-2)}{1-\alpha}}\cdot\int_{x=U/L}^{L}x^{-2}\mathrm{d}x\Big{)}^{1/2}$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
t\cdot\Big{(}c\cdot(U/L)^{\frac{2(\omega-2)}{1-\alpha}-1}+(U/L)^{2(\frac{\omega-2}{1-\alpha})}\cdot(U/L)^{-1}\Big{)}^{1/2}$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\ t\cdot
O((U/L)^{\frac{(\omega-2)}{1-\alpha}-1/2}),$

where the first step follows from the definition of $g\in\mathbb{R}^{L}$, the
second step follows from the assumption $\sum_{t=1}^{t}r_{i}\cdot
g_{r_{i}}\leq t\cdot\|g\|_{2}$ in the lemma statement, the third step follows
from upper bounding the $\ell_{2}$ norm
$\|g\|_{2}^{2}=\sum_{r=1}^{L}g_{r}^{2}$, and the fourth step follows
$\frac{2(\omega-2)}{1-\alpha}\geq 1$ when $\alpha\geq 5-2\omega$, so the
integral $\int_{x=1}^{U/L}x^{\frac{2(\omega-2)}{1-\alpha}-2}\mathrm{d}x=c\cdot
x^{\frac{2(\omega-2)}{1-\alpha}-1}\big{|}_{1}^{U/L}=O\big{(}(U/L)^{\frac{2(\omega-2)}{1-\alpha}-1}\big{)}$
where $c:=1/(\frac{2(\omega-2)}{1-\alpha}-1)$.

Thus we have

$\displaystyle\sum_{t=1}^{t}{\cal
T}_{\mathrm{mat}}(U,U,\min\\{Lr_{i},U\\})\leq$
$\displaystyle\leavevmode\nobreak\
\sum_{t=1}^{t}\Big{(}U^{2+o(1)}+U^{2-\frac{\alpha(\omega-2)}{1-\alpha}+o(1)}\cdot
L^{\frac{\omega-2}{1-\alpha}}\cdot\min\\{r_{i}^{\frac{\omega-2}{1-\alpha}},(U/L)^{\frac{\omega-2}{1-\alpha}}\\}\Big{)}$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\ t\cdot
U^{2+o(1)}+U^{2-\frac{\alpha(\omega-2)}{1-\alpha}+o(1)}\cdot
L^{\frac{\omega-2}{1-\alpha}}\cdot\sum_{t=1}^{t}\min\\{r_{i}^{\frac{\omega-2}{1-\alpha}},(U/L)^{\frac{\omega-2}{1-\alpha}}\\}$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\ t\cdot
U^{2+o(1)}+U^{2-\frac{\alpha(\omega-2)}{1-\alpha}+o(1)}\cdot
L^{\frac{\omega-2}{1-\alpha}}\cdot
t\cdot(U/L)^{\frac{(\omega-2)}{1-\alpha}-1/2}$ $\displaystyle=$
$\displaystyle\leavevmode\nobreak\ t\cdot(U^{2+o(1)}+U^{\omega-1/2+o(1)}\cdot
L^{1/2}),$

where the first step follows from Eq. (8) and ${\cal
T}_{\mathrm{mat}}(U,U,U)=U^{\omega}=U^{2-\frac{\alpha(\omega-2)}{1-\alpha}}\cdot
L^{\frac{\omega-2}{1-\alpha}}\cdot(U/L)^{\frac{\omega-2}{1-\alpha}}$, the
second step follows from moving summation inside, the third step follows from
Eq. (9), and the last step follows from adding the terms together. ∎

##  4 Introduction

Polynomial optimization is a fundamental problem in many areas of applied
mathematics, operations research, and theoretical computer science, including
combinatorial optimization [BS11, LasserreNing, BHKKMP19], statistical
estimation [HopkinsKPRSS17, Hopkins018], experimental design [Papp12], control
theory [HessHLP16], signal processing [RohDV07], power systems engineering
[GPF], discrete geometry [BV08, BBBCGKS09] and computational algebraic
geometry [AG09]. In the most basic formulation, we are given a collection of
$k$ real $n$-variate polynomials $g_{1},\cdots,g_{k}$ and an objective
function $f:\mathbb{R}^{n}\to\mathbb{R}$, and the goal is to minimize $f$ over
the set $\mathcal{S}:=\\{t\in\mathbb{R}^{n}\mid\forall
i\in\\{1,\cdots,k\\}:g_{i}(t)\geq 0\\}$, that is, to find

$\inf_{t\in\mathbb{R}^{n}}\\{f(t)\mid t\in\mathcal{S}\\},$ (10)

which is equivalent to checking polynomial nonnegativity
$\sup_{c\in\mathbb{R}}\\{c\mid f(t)-c\geq 0,\leavevmode\nobreak\ \forall
t\in\mathcal{S}\\}$. This is then equivalent to computing
$\sup_{c\in\mathbb{R}}\\{c\mid f-c\in\mathcal{K}(\mathcal{S})\\}$, where
$\mathcal{K}(\mathcal{S})$ denotes the convex cone of all polynomials of
degree at most $\deg(f)$ that are non-negative on the set $\mathcal{S}$. This
is an instance of the more general conic programming:

$\displaystyle\min_{x\in\mathbb{R}^{N}}$ $\displaystyle c^{\top}x\quad$ (11)
$\displaystyle Ax$ $\displaystyle=b$ $\displaystyle x$
$\displaystyle\in\mathcal{K},\ $

where $\mathcal{K}\subset\mathbb{R}^{N}$ is some convex cone111A subset
$\mathcal{K}\subset\mathbb{R}^{N}$ is a convex cone if
$\forall\;x,y\in\mathcal{K}$ and $\alpha,\beta\in\mathbb{R}_{+}$, $\alpha
x+\beta y\in\mathcal{K}$.. The conic optimization problem over the cone
$\mathcal{K}(\mathcal{S})$ is intractable in general because there is no
simple characterization of $\mathcal{K}(\mathcal{S})$. Nevertheless, there
always exists an increasing family of convex cones of _weighted sum-of-
squares_ polynomials that converges to any such cone
$\mathcal{K}(\mathcal{S})$.

We first introduce the notion of _sum-of-squares_ (SOS) polynomials: Denoting
by $\mathcal{V}_{n,d}$ the vector space of all $n$-variate polynomials of
(total) degree at most $d$, a polynomial $p\in\mathcal{V}_{n,2d}$ is said to
be _sum-of-squares_ (SOS) if it can be written as a finite sum of square
polynomials, i.e., there exist $q_{1},\cdots,q_{\ell}$ such that
$p=\sum_{i=1}^{\ell}q_{i}^{2}$. The set $\Sigma_{n,2d}$ of SOS polynomials of
degree at most $2d$ is a (proper) cone contained in $\mathcal{V}_{n,2d}$, of
dimension $U:=\dim(\mathcal{V}_{n,2d})={n+2d\choose d}$, as the vector space
$\mathcal{V}_{n,2d}$ is isomorphic to $\mathbb{R}^{U}$. If $p$ can be written
as $p=\sum_{i=1}^{k}g_{i}s_{i}$ for
$s_{1}\in\Sigma_{n,2d_{1}},\cdots,\Sigma_{n,2d_{k}}$ and $k$ nonzero
polynomials $\mathbf{f}:=(f_{1},\cdots,f_{k})$, then it is said to be
_weighted sum-of-squares_ (WSOS).

_Putinar’s Positivstellensatz_ [PUTINAR99] states that under mild conditions,
any polynomial $p$ that is non-negative on $\mathcal{S}$ can be written as a
WSOS polynomial $\sum_{i=1}^{k}g_{i}s_{i}$, albeit with (potentially)
_unbounded degree_ $s_{i}$’s. In WSOS optimization we consider sum-of-squares
polynomials $s_{i}$ with bounded degree, so the hierarchy of WSOS optimization
with increasing degree (known as the Lasserre hierarchy) can be viewed as a
tool for approximating general polynomial optimization. For more details of
this approximation scheme for polynomial optimization, we refer the readers to
the matextbooks [l15, bpt13].

This paper concerns algorithms for _(W)SOS optimization_ , which is the conic
optimization program (11) where the underlying cone $\mathcal{K}$ is the
(W)SOS cone:

$\displaystyle\min_{x\in\mathbb{R}^{U}}$ $\displaystyle c^{\top}x\quad$ (12)
$\displaystyle Ax$ $\displaystyle=b$ $\displaystyle x$
$\displaystyle\in\Sigma_{n,2d},\ $

where $x\in\mathcal{V}_{n,2d}$ is the vector of coefficients which encodes the
polynomial. Henceforth, we focus on the case where $\mathcal{K}=\Sigma_{n,2d}$
is the SOS cone, and discuss how to extend our algorithm for SOS optimization
to WSOS in Section 11.

The computational complexity of solving Problem 12 naturally depends on the
dimensions

$L:=\dim(\mathcal{V}_{n,d})={n+d\choose
d}\;\;\;\;\;,\;\;\;\;U:=\dim(\mathcal{V}_{n,2d})={n+2d\choose 2d}$ (13)

of the underlying vector spaces (Note that $L\leq U\leq L^{2}$). We now turn
to explain the previous approaches for SOS optimization solvers.

#### SOS Optimization as SDPs

A fundamental fact is that the _dual_ SOS cone is a _slice of the SDP cone_
[nesterov2000squared]. More formally, for any _fixed bases_
$\mathbf{p}=(p_{1},p_{2},\cdots,p_{L})$ and
$\mathbf{q}=(q_{1},q_{2},\cdots,q_{U})$ to $\mathcal{V}_{n,d}$ and
$\mathcal{V}_{n,2d}$ respectively, there exists a unique linear mapping
$\Lambda:\mathbb{R}^{U}\to\mathbb{R}^{L\times L}$ satisfying

$\Lambda(\mathbf{q}(t))=\mathbf{p}(t)\mathbf{p}(t)^{\top},\leavevmode\nobreak\
\leavevmode\nobreak\ \forall t\in\mathbb{R}^{n}.$ (14)

Here we define $\mathbf{p}(t)=(p_{1}(t),p_{2}(t),\cdots,p_{L}(t))^{\top}$ and
$\mathbf{q}(t)=(q_{1}(t),q_{2}(t),\cdots,q_{U}(t))^{\top}$. An equivalent way
to view the definition of $\Lambda$ in (14) is as follows: For polynomials
$p_{i},p_{j}\in\mathbf{p}$ there are unique coefficients $\lambda_{iju}$ such
that $p_{i}p_{j}=\sum_{u\in U}\lambda_{iju}q_{u}$. These $\lambda_{iju}$
define the mapping $\Lambda$ unambiguously.

This in turn implies that a polynomial $s\in\mathcal{V}_{n,2d}$ (we view $s$
as a vector in $\mathbb{R}^{U}$ that corresponds to its coefficients over the
basis $\mathbf{q}$) is in the dual SOS cone $\Sigma_{n,2d}^{*}$ if and only if
$\Lambda(s)$ is a _positive semidefinite_ (PSD) matrix (proved by
[nesterov2000squared], see Theorem 6.4 for details). As [py19] recently
observed, the choice of the bases $\mathbf{p},\mathbf{q}$ crucially affects
the complexity of the optimization problem, more on this below.

Equation (14) implies the well-known fact that optimization over SOS
polynomials (12) can be reduced to semidefinite programming

$\min_{X\succeq 0}\\{\langle
C,X\rangle\mid\operatorname{tr}(A_{i}X)=b_{i},\leavevmode\nobreak\ \forall
i\in[m]\\},$ (SDP)

and can thus be solved using off-the-shelve SDP solvers. However, despite
recent breakthroughs on the runtime of general SDP solvers via _interior-point
methods_ (IPMs) [jiang2020faster, hjst21], the SDP reformulation (SDP) of (12)
does not scale well for moderately large degrees, i.e., whenever $U\ll L^{2}$
in (13). This is because the SDP reformulation always incurs a factor of at
least $L^{2}$, even when $U\ll L^{2}$, as this is the SDP variable size (the
PSD matrix $X$ has size $L\times L$). Indeed, for current fast-matrix-
multiplication time $\omega\approx 2.37$ [l14, aw21], the running time of
state-of-the-art SDP solvers [jiang2020faster, hjst21] for SOS optimization
(Problem 12) is222We use $\widetilde{O}(\cdot)$ to hide $U^{o(1)}$ and
$\log(1/\delta)$ factors.

$\widetilde{O}\left(L^{0.5}\cdot\min\\{UL^{2}+U^{2.37},\leavevmode\nobreak\
L^{4.24}\\}\right).$ (15)

An alternative approach is to solve Problem 12 directly by designing an _ad-
hoc_ IPM for the dual SOS cone, avoiding the blowup in the SDP reformulation.
This was exactly the motivation of [py19]. In more detail, all aformentioned
SDP solvers are based on IPMs [nesterov1994interior], which iteratively
minimize the original objective function plus a barrier function via Newton
steps. When applied to the SOS Problem (12), _the choice of the specific bases
$\mathbf{p},\mathbf{q}$ crucially affects the structure of the (Hessian of
the) barrier function $F(s)=F(\Lambda(s))$_, and hence the cost-per-iteration
of the IPM. As such, choosing a “good” and efficient basis is key to a fast
algorithm for (12). One of the main contribution of [py19] is an efficient
basis for the SOS cone, which efficiently scales to multivariate SOS, yielding
an IPM whose total runtime is

$\widetilde{O}\left(L^{0.5}U^{\omega}\right)\approx\widetilde{O}\left(L^{0.5}U^{2.37}\right).$
(16)

Our main result is a polynomially faster IPM for Problem 12:

######  Theorem 4.1 (Main Result, Informal version of Theorem 9.3).

With current FMM exponent, there is an algorithm for solving Problem (12),
whose total running time is $\widetilde{O}\left(LU^{1.87}\right).$

Indeed, this runtime is polynomially faster than (16) and (15), as shown in
Figure 1. We now turn to elaborate on the technical approach for proving
Theorem 4.1.

$1$$1.2$$1.4$$1.6$$1.8$$2$$2.4$$2.6$$2.8$$3$$3.2$$3.4$$3.6$$\log_{L}U$$\log_{U}(T)$This
paperPY19SDPlower bound

Figure 1: Overview of current running times of recent solvers for SOS. The
lower bound bound stems from solving a linear system in $U$ variables, i.e.,
$T=\Omega(U^{\omega})$ where $\omega\approx 2.37$.

.

#### Faster IPMs via Inverse-Maintenance

Interior-Point Methods (IPMs [karmarkar1984new,
renegarMathematicalViewInteriorPoint2001]) are a powerful class of second-
order optimization algorithms for convex optimization, which essentially
reduce a conic optimization problem (11) to solving _a sequence of slowly-
changing linear systems_ (via Newton steps). Since their discovery in the mid
80’s, IPMs have emerged as the “gold-standard” of convex optimization, as they
are known to converge fast in _both theory and practice_ [s87]. The main
computational cost of IPMs is computing, in each iteration, the inverse of the
Hessian of the underlying barrier function $F(s)=F(\Lambda(s))$, which naively
costs at least $U^{\omega}$ time per iteration for the SOS optimization
problem [py19]. A recent influential line of work [cls19, jiang2020faster],
inspired by [v89_lp]’s seminal work, has demonstrated that _dynamically
maintaining_ the inverse of the Hessian matrix under _low-rank_ updates using
clever _data structures_ , can lead to much cheaper _cost-per-iteration_. All
of these results rely on a careful combination of dynamic data structures with
the geometry (e.g., spectral approximation) of the underlying optimization
method and barrier function. This paper extends this line of work to SOS
optimization.

#### Our Techniques

We follow the framework of [py19] which chooses the _polynomial interpolant
basis_ representation and the corresponding linear operator
$\Lambda:\mathbb{R}^{U}\to\mathbb{R}^{L\times L}$ is
$\Lambda(s)=P^{\top}\operatorname{diag}(s)P$, where $P\in\mathbb{R}^{U\times
L}$ is the matrix whose entries are the evaluation of the Lagrange
interpolation polynomials, through some unisolvent333Any set of points in
$\mathbb{R}^{n}$ for which the evaluation of a polynomial in
$\mathcal{V}_{n,d}$ on these points uniquely defines the polynomial. set of
points in $\mathcal{V}_{n,d}$ (see Section 6 for a formal definition). This
basis induces the aforementioned convenient form of $\Lambda$, and generalizes
to the multivariate case. The Hessian of the barrier function
$F(s)=-\log\det(\Lambda(s))$ is given by

$H(s)=\big{(}P(P^{\top}\operatorname{diag}(s)P)^{-1}P^{\top}\big{)}^{\circ
2}\in\mathbb{R}^{U\times U},$

where $A\circ B$ denotes the element-wise (Hadamard) product of two matrices.
The main bottleneck of each iteration of IPMs is to compute the Hessian
inverse $H(x)^{-1}$ of the Newton step, which naïvely takes $O(U^{\omega})$
time.

In IPM theory, it has long been known that it suffices to compute a spectral
approximation of the Hessian. We follow the “lazy update” framework in recent
developments of LP and SDP solvers [cls19, jiang2020faster], which batches
together low-rank updates to
$M:=P(P^{\top}\operatorname{diag}(s)P)^{-1}P^{\top}$, where
$\operatorname{rk}(M)=L$. In each iteration, we can compute a spectral
approximation $M^{\operatorname{new}}=M+UV^{\top}$, where $U,V$ are low rank
matrices with size $U\times r$ where $r\ll U$ is chosen to optimize the
runtime. Since $\widetilde{M}\approx M$ implies that $\widetilde{M}^{\circ
2}\approx M^{\circ 2}$, this also gives a spectral approximation of the
Hessian.

The main challenge here, compared to previous LP and SDP solvers [v89_lp,
ls14, cls19, jiang2020faster, hjst21], is that low-rank updates to $M$ do not
readily translate to a low-rank update to $(M^{\circ 2})^{-1}$, since
Hadamard-products can _increase_ the rank $\operatorname{rk}(A\circ
B)\leq\operatorname{rk}(A)\cdot\operatorname{rk}(B)$, in contrast to standard
matrix multiplication which does not increase the rank
$\operatorname{rk}(AB)\leq\max\\{\operatorname{rk}(A),\operatorname{rk}(B)\\}$.
This means that we cannot directly apply Woodbury’s identity to efficiently
update the inverse of the Hessian, which is the common approach in all
aforementioned works. Instead, we employ the following property which relates
rank-one Hadamard-product perturbations to standard matrix products

$M\circ(u\cdot v^{\top})=\operatorname{diag}(u)\cdot
M\cdot\operatorname{diag}(v),$

which means that we can translate the rank-$r$ update of $M$ into a rank-$Lr$
update of $M^{\circ 2}$ for $r\leq L$. With some further calculations,
applying Woodbury’s identity on the resulting matrix, implies that we can
compute $((M^{\operatorname{new}})^{\circ 2})^{-1}$ in time

$O\big{(}{\cal T}_{\mathrm{mat}}(U,U,Lr)\big{)},$

which is never worse than ${\cal T}_{\mathrm{mat}}(U,U,U)=U^{\omega}$ as long
as $r\leq U/L$. Modifying the amortization tools of [jiang2020faster] and
[hjst21], combined with basic spectral theory for Hadamard products, we show
that our amortized cost per iteration is bounded by

$O\big{(}U^{2}+U^{\omega-1/2}\cdot L^{1/2}\big{)},$

which becomes $O\big{(}U^{2}+U^{1.87}L^{0.5}\big{)}$ if we plug in the current
matrix multiplication exponent.

##  5 Preliminaries

In this section we provide the definitions and the tools that we will use. For
any integer $n>0$, we define $[n]=\\{1,2,\cdots,n\\}$. We use $\mathbb{R}_{+}$
and $\mathbb{R}_{\geq 0}$ to denote the set of positive and non-negative real
numbers respectively. We use $0_{n},1_{n}\in\mathbb{R}^{n}$ to denote the all-
zero and all-one vectors of size $n$.

Given a vector $v\in\mathbb{R}^{n}$, for any $m\leq n$, we use
$v_{[:m]}\in\mathbb{R}^{m}$ to denote the first $m$ entries of $v$. For a
vector $v\in\mathbb{R}^{n}$, we use
$\operatorname{diag}(v)\in\mathbb{R}^{n\times n}$ to denote the diagonal
matrix whose diagonal entries are $v$. For a square matrix
$A\in\mathbb{R}^{n\times n}$, we use $\operatorname{diag}(A)\in\mathbb{R}^{n}$
to denote the vector of the diagonal entries of $A$. We use
$\operatorname{rk}(A)$ to denote the rank of a matrix $A$. We use $\ker(A)$
and $\operatorname{Im}(A)$ to denote the kernel space and the column space of
$A$.

We say a matrix $A\in\mathbb{R}^{n\times n}$ is PSD (denoted as $A\succeq 0$)
if $A$ is symmetric and $x^{\top}Ax\geq 0$ for all $x\in\mathbb{R}^{n}$. We
use $\mathbb{S}^{n\times n}$ to denote the set of PSD matrices of size
$n\times n$. The spectral norm of a matrix $A\in\mathbb{R}^{n\times d}$ is
defined as $\|A\|_{2}=\max_{x\in\mathbb{R}^{d},\|x\|_{2}=1}\|Ax\|_{2}$. The
Frobenius norm of $A$ is defined as
$\|A\|_{F}=\sqrt{\sum_{i\in[n]}\sum_{j\in[d]}A_{i,j}^{2}}$. For any PSD matrix
$M\in\mathbb{S}^{n\times n}$, we define the $M$-norm as
$\|x\|_{M}=\sqrt{x^{\top}Mx}$, $\forall x\in\mathbb{R}^{n}$.

We use ${\cal T}_{\mathrm{mat}}(a,b,c)$ to denote the time to multiply two
matrices of sizes $a\times b$ and $b\times c$. A basic fact of fast matrix
multiplication is that ${\cal T}_{\mathrm{mat}}(a,b,c)={\cal
T}_{\mathrm{mat}}(b,c,a)={\cal T}_{\mathrm{mat}}(c,a,b)$ (see e.g. [b13]), and
we will use these three terms interchangeably.

######  Fact 5.1 (Woodbury identity).

Let $A\in\mathbb{R}^{n\times n},C\in\mathbb{R}^{k\times
k},U\in\mathbb{R}^{n\times k},V\in\mathbb{R}^{k\times n}$ where $A$ and $C$
are invertible, then

$(A+UCV)^{-1}=A^{-1}-A^{-1}U(C^{-1}+VA^{-1}U)^{-1}VA^{-1}.$

######  Definition 5.2 (Hadamard product).

For any two matrices $A,B\in\mathbb{R}^{m\times n}$, the Hadamard product
$A\circ B$ is defined as

$(A\circ B)_{i,j}=A_{i,j}\cdot B_{i,j},\leavevmode\nobreak\
\leavevmode\nobreak\ \forall i\in[m],j\in[n].$

We also use $A^{\circ 2}$ to denote $A\circ A$.

The Hadamard product has the following properties (the proofs are
straightforward).

######  Fact 5.3 (Properties of Hadamard product).

For matrices $A,B\in\mathbb{R}^{m\times n}$, and vectors $x\in\mathbb{R}^{m}$,
$y\in\mathbb{R}^{n}$, we have the following properties.

  1. 1.

$x^{\top}(A\circ
B)y=\operatorname{tr}[\operatorname{diag}(x)A\operatorname{diag}(y)B^{\top}]$,

  2. 2.

$A\circ(x\cdot y^{\top})=\operatorname{diag}(x)\cdot
A\cdot\operatorname{diag}(y)$.

######  Definition 5.4 (Spectral approximation).

For any two symmetric matrices $A,\widetilde{A}\in\mathbb{R}^{n\times n}$, any
parameter $\epsilon\in(0,1)$, we say $\widetilde{A}$ and $A$ are
$\epsilon$-spectral approximation of each other, denoted as
$\widetilde{A}\approx_{\epsilon}A$, if we have

$e^{-\epsilon}\cdot x^{\top}Ax\leq x^{\top}\widetilde{A}x\leq
e^{\epsilon}\cdot x^{\top}Ax,\leavevmode\nobreak\ \leavevmode\nobreak\ \forall
x\in\mathbb{R}^{n}.$

Spectral approximation has the following properties (for completeness we
include a proof in the appendix). [Properties of spectral approximation] For
any two PSD matrices $A,\widetilde{A}\in\mathbb{R}^{n\times n}$, any parameter
$\epsilon\in(0,1)$, if $\widetilde{A}\approx_{\epsilon}A$, then we have

  1. 1.

$B^{\top}AB\approx_{\epsilon}B^{\top}\widetilde{A}B$, for any matrix
$B\in\mathbb{R}^{n\times n}$.

  2. 2.

If both $A$ and $\widetilde{A}$ are invertible, then
$A^{-1}\approx_{\epsilon}\widetilde{A}^{-1}$.

  3. 3.

$e^{-\epsilon}\operatorname{tr}[A]\leq\operatorname{tr}[\widetilde{A}]\leq
e^{\epsilon}\operatorname{tr}[A]$.

  4. 4.

$\widetilde{A}^{\circ 2}\approx_{2\epsilon}A^{\circ 2}$.

###### Proof.

The proofs of the first three claims are straightforward. We only prove the
last claim.

For any vector $x\in\mathbb{R}^{n}$, we have

$\displaystyle x^{\top}A^{\circ 2}x=$ $\displaystyle\leavevmode\nobreak\
\operatorname{tr}[\operatorname{diag}(x)A\operatorname{diag}(x)A]$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\operatorname{tr}[A^{1/2}\operatorname{diag}(x)A\operatorname{diag}(x)A^{1/2}]$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
e^{\epsilon}\cdot\operatorname{tr}[A^{1/2}\operatorname{diag}(x)\widetilde{A}\operatorname{diag}(x)A^{1/2}]$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
e^{\epsilon}\cdot\operatorname{tr}[\widetilde{A}^{1/2}\operatorname{diag}(x)A\operatorname{diag}(x)\widetilde{A}^{1/2}]$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
e^{2\epsilon}\cdot\operatorname{tr}[\widetilde{A}^{1/2}\operatorname{diag}(x)\widetilde{A}\operatorname{diag}(x)\widetilde{A}^{1/2}]=e^{2\epsilon}\cdot
x^{\top}\widetilde{A}^{\circ 2}x$

where the first step follows from Fact 5.3, the second and the fourth steps
follow from the trace invariance under cyclic permutations and the fact that
$A^{1/2}$ exists when $A$ is PSD, the third and the fifth steps follow from
Part 3 of this fact.

Similarly we can prove $x^{\top}A^{\circ 2}x\geq e^{-2\epsilon}\cdot
x^{\top}\widetilde{A}^{\circ 2}x$. Thus we have $A^{\circ
2}\approx_{2\epsilon}\widetilde{A}^{\circ 2}$. ∎

##  6 Background of sum-of-squares optimization

In this section we provide the background of sum-of-squares optimization. We
refer the readers to [parrilo2020sum, py19] for more details.

######  Definition 6.1 (Polynomial space).

We use $\mathcal{V}_{n,d}$ to denote the set of $n$-variate polynomials over
the reals of degree at most $d$, where the degree means the total degree,
i.e., the degree of $x_{1}^{d_{1}}\cdots x_{n}^{d_{n}}$ is
$\sum_{i=1}^{n}d_{i}$.

######  Definition 6.2 (Degree of polynomial space).

We define $L:=\dim(\mathcal{V}_{n,d})=\binom{n+d}{n}$ and
$U:=\dim(\mathcal{V}_{n,2d})=\binom{n+2d}{n}$.

After fixing a basis $(p_{1},p_{2},\cdots,p_{L})$ of $\mathcal{V}_{n,d}$,
there exists a one-to-one correspondence between any polynomial
$p=\sum_{i=1}^{L}x_{i}\cdot p_{i}\in\mathcal{V}_{n,d}$ and the vector
$[x_{1},x_{2},\cdots,x_{L}]\in\mathbb{R}^{L}$. From now on when the basis is
clear from context, we will use $\mathcal{V}_{n,d}$ and $\mathbb{R}^{L}$
interchangeably, and similarly $\mathcal{V}_{n,2d}$ and $\mathbb{R}^{U}$
interchangeably.

######  Definition 6.3 (SOS polynomials).

A polynomial $p\in\mathcal{V}_{n,2d}$ is said to be a sum-of-squares $(SOS)$
polynomial if $p$ can be written as a sum of squares of polynomials, i.e.
$p=\sum_{i=1}^{M}q_{i}^{2}$ for some $M\in\mathbb{N}$ and polynomials
$q_{1},q_{2},\cdots,q_{M}\in\mathcal{V}_{n,d}$.

We use $\Sigma_{n,2d}$ to denote the set of $n$-variate SOS polynomials of
degree at most $2d$.

The set $\Sigma_{n,2d}$ is a closed convex and pointed cone in
$\mathcal{V}_{n,2d}$ with non-empty interior (Theorem 17.1 of
[nesterov2000squared]). The SOS optimization problem requires the variable
$x\in\mathbb{R}^{U}$ to be in the SOS cone, and it is a special case of conic
programming. Given a constraint matrix $A\in\mathbb{R}^{m\times U}$ where
$m\leq U$, and $b\in\mathbb{R}^{m}$ and $c\in\mathbb{R}^{U}$, the SOS
optimization can be written in the following primal-dual formulation:

$\begin{aligned} \text{Primal:\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ }\min\;&\langle c,x\rangle\quad\\\
\mathrm{s.t.}\leavevmode\nobreak\ \leavevmode\nobreak\ Ax&=b\\\
x&\in\Sigma_{n,2d}\,,\\\ \end{aligned}\quad\quad\quad\begin{aligned}
\text{Dual:\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
}\max\;&\langle y,b\rangle\\\ \mathrm{s.t.}\leavevmode\nobreak\
\leavevmode\nobreak\ A^{\top}y+s&=c\\\ s&\in\Sigma_{n,2d}^{*}\,.\\\
\end{aligned}$ (SOS)

Here $\Sigma_{n,2d}^{*}:=\\{s\in\mathbb{R}^{U}\mid s^{\top}x\geq
0,\leavevmode\nobreak\ \forall x\in\Sigma_{n,2d}\\}$ denotes the dual cone of
$\Sigma_{n,2d}$.

Nesterov in [nesterov2000squared] noted that the dual SOS cone allows the
following characterization.

######  Theorem 6.4 (Dual cone characterization, Theorem 17.1 of
[nesterov2000squared]).

For any ordered bases $\mathbf{p}=(p_{1},\ldots,p_{L})$ and
$\mathbf{q}=(q_{1},\ldots,q_{U})$ of $\mathcal{V}_{n,d}$ and
$\mathcal{V}_{n,2d}$, let $\Lambda:\mathbb{R}^{U}\to\mathbb{R}^{L\times L}$ be
the unique linear mapping satisfying
$\Lambda(\mathbf{q})=\mathbf{p}\mathbf{p}^{\top}$.444This equation means
$\forall t\in\mathbb{R}^{n}$,
$\Lambda([q_{1}(t),\cdots,q_{U}(t)]^{\top})=[p_{1}(t),\cdots,p_{L}(t)]^{\top}\cdot[p_{1}(t),\cdots,p_{L}(t)]$.
Then the dual cone $\Sigma_{n,2d}^{*}$ admits the characterization under the
bases $\mathbf{p}$ and $\mathbf{q}$:

$\Sigma_{n,2d}^{*}=\big{\\{}s\in\mathbb{R}^{U}\mid\Lambda(s)\succeq
0\big{\\}}.$ (17)

As barrier functions for the cone of positive semidefinite matrices are well-
known, this also gives rise to a barrier function for the dual SOS cone. With
the standard log-det barrier for the semidefinite cone, the following function
$F:\Sigma_{n,2d}^{*}\to\mathbb{R}$ is a barrier function for
$\Sigma_{n,2d}^{*}$:

$F(s)=-\log\det(\Lambda(s)).$

Furthermore, the barrier parameter $\nu_{F}$ of $F(s)$ is bounded by the
barrier parameter $L$ of the original log-det barrier function
([nesterov2000squared]).

#### Interpolant basis

The barrier function depends on the choice of the basis for both
$\mathcal{V}_{n,d}$ and $\mathcal{V}_{n,2d}$, as the linear map $\Lambda$
depends on these two bases. We follow the approach of [py19] and focus on the
so-called _interpolant_ bases, which generalises well to multivariate
polynomials and is numerically stable.

For the vector space $\mathcal{V}_{n,2d}$, consider a set of _unisolvent_
points $\mathcal{T}=\\{t_{1},t_{2},\cdots,t_{U}\\}\subseteq\mathbb{R}^{n}$,
which is a set points such that every polynomial in $\mathcal{V}_{n,2d}$ is
uniquely determined by its values on the points in $\mathcal{T}$. For
univariate polynomials any set of $U$ points suffices, but this does not hold
anymore for the multivariate case. To also ensure numerical stability, the so
called (approximate) Fekete points can be used as unisolvent points
[sommariva09, sommariva10].

The interpolant basis is defined as follows. Let us fix a set of unisolvent
points $\mathcal{T}=\\{t_{1},t_{2},\cdots,t_{U}\\}\subseteq\mathbb{R}^{n}$.
Now every $t_{u}\in\mathcal{T}$ implies a Lagrange polynomial $q_{u}$ which is
the unique polynomial that satisfies $q_{u}(t_{u})=1$ and $q_{u}(t_{v})=0$ for
all $t_{v}\neq t_{u}\in\mathcal{T}$. The Lagrange polynomials form a basis
$\mathbf{q}=(q_{1},\cdots,q_{U})$ of $\mathcal{V}_{n,2d}$. Choose any basis
$\mathbf{p}=(p_{1},\ldots,p_{L})$ of $\mathcal{V}_{n,d}$. Define the matrix
$P\in\mathbb{R}^{U\times L}$ as

$P_{u,\ell}=p_{\ell}(t_{u}),\leavevmode\nobreak\ \leavevmode\nobreak\ \forall
u\in[U],\ell\in[L].$

By the definition of the Lagrange polynomials,
$p_{i}p_{j}=\sum_{u=1}^{U}p_{i}(t_{u})p_{j}(t_{u})q_{u}$, so we have
$\mathbf{p}\mathbf{p}^{\top}=P^{\top}\operatorname{diag}(\mathbf{q})P$. Thus
under the bases $\mathbf{p}$ and $\mathbf{q}$, the linear map
$\Lambda:\mathbb{R}^{U}\to\mathbb{R}^{L\times L}$ takes on the following
convenient form:

$\Lambda(s)=P^{\top}\operatorname{diag}(s)P.$ (18)

##  7 Algorithm

Since in this paper we focus on the theoretical running time of the algorithm,
for simplicity we use the barrier method (see e.g.
[renegarMathematicalViewInteriorPoint2001, Chapter 2]) instead of the more
sophisticated Skajaa–Ye Algorithm used by [py19].

The dual formulation of (SOS) is equivalent to the following optimization
problem

$\min-b^{\top}y\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
s.t.\leavevmode\nobreak\ \leavevmode\nobreak\ y\in\overline{D}_{F},$

where with an abuse of the notation we define
$F:\mathbb{R}^{m}\to\mathbb{R}_{+}$ to be the barrier function

$F(y)=-\log\det(\Lambda(c-A^{\top}y))$ (19)

for $c\in\mathbb{R}^{U}$, and $A\in\mathbb{R}^{m\times U}$, and
$\Lambda(s)=P^{\top}\operatorname{diag}(s)P$ is the linear operator defined in
Eq (18). $D_{F}\subseteq\mathbb{R}^{m}$ is the domain of $F$, and
$\overline{D}_{F}$ is the closure of $D_{F}$.

The barrier parameter of the barrier function $F$ is $\nu_{F}=L$. The gradient
and the Hessian of the barrier function $F$ are (define $s:=c-A^{\top}y$):

$\displaystyle g(y)=$ $\displaystyle\leavevmode\nobreak\
A\cdot\operatorname{diag}\Big{(}P\big{(}P^{\top}\operatorname{diag}(s)P\big{)}^{-1}P^{\top}\Big{)},$
$\displaystyle H(y)=$ $\displaystyle\leavevmode\nobreak\
A\cdot\Big{(}P\big{(}P^{\top}\operatorname{diag}(s)P\big{)}^{-1}P^{\top}\Big{)}^{\circ
2}\cdot A^{\top}.$

For any $\eta>0$, define a function $F_{\eta}:\mathbb{R}^{m}\to\mathbb{R}$:

$F_{\eta}(y)=-\eta\cdot b^{\top}y+F(y).$

The gradient and the Hessian of $F_{\eta}(y)$ are:

$\displaystyle g_{\eta}(y)=$ $\displaystyle\leavevmode\nobreak\ -\eta\cdot
b+A\cdot\operatorname{diag}\Big{(}P\big{(}P^{\top}\operatorname{diag}(s)P\big{)}^{-1}P^{\top}\Big{)},$
$\displaystyle H_{\eta}(y)=$ $\displaystyle\leavevmode\nobreak\
A\cdot\Big{(}P\big{(}P^{\top}\operatorname{diag}(s)P\big{)}^{-1}P^{\top}\Big{)}^{\circ
2}\cdot A^{\top}.$

Note that $H_{\eta}(y)=H(y)$ for any $\eta$.

In each iteration the barrier method increases $\eta$ by a factor of
$1+\Theta(\frac{1}{\sqrt{L}})$, and it performs a Newton step

$y\leftarrow y-H_{\eta}(y)^{-1}\cdot g_{\eta}(y).$

By standard IPM theory it suffices to use a spectral approximation of the
Hessian matrix in the Newton step. For more details see e.g.
[renegarMathematicalViewInteriorPoint2001].

The main technical part of our algorithm is to efficiently maintain a matrix
$N$ that is the spectral approximation of the inverse of the Hessian matrix.
To do this, we maintain another matrix $\widetilde{S}$ that is a spectral
approximation of $S:=P^{\top}\operatorname{diag}(s)P$, and we use the
subroutine LowRankUpdate(Algorithm 3, Lemma 9.1) to update $\widetilde{S}$.
After $\widetilde{S}$ is updated, we use another subroutine UpdateHessianInv
(Algorithm 2, Lemma 8.1) to update $N$. A complete description of our
algorithm can be found in Algorithm 1.

Parameters : $\delta\in(0,1)$, $\epsilon_{N}\in(0,0.05)$,
$\alpha=\frac{\epsilon_{N}}{20\sqrt{L}}$,
$t=40\epsilon_{N}^{-1}\sqrt{L}\log(L/\delta)$.

Input : $A\in\mathbb{R}^{m\times U}$, $b\in\mathbb{R}^{m}$,
$c\in\mathbb{R}^{U}$

Output : A near feasible and optimal solution.

1

2Construct $P\in\mathbb{R}^{U\times L}$ of the interpolant basis. Convert
$A,b,c$ to the interpolant basis.

3Use Lemma 1.1 to obtain a modified dual SOS optimization problem which has an
initial solution $(y,s)\in\mathbb{R}^{m}\times\mathbb{R}^{U}$ that is optimal
for $F_{\eta}$, where $\eta=1$.

4$\widetilde{S}\leftarrow S\leftarrow P^{\top}\operatorname{diag}(s)P$ ;

5 // $\widetilde{S},S\in\mathbb{R}^{L\times L}$

6

7$T\leftarrow S^{-1}$ ;

8 // $T\in\mathbb{R}^{L\times L}$

9

10$N\leftarrow\big{(}A(PTP^{\top})^{\circ 2}A^{\top}\big{)}^{-1}$ ;

11 // $N\in\mathbb{R}^{m\times m}$

12

13$g\leftarrow-\eta\cdot
b+A\cdot\operatorname{diag}\big{(}P(P^{\top}\operatorname{diag}(s)P)^{-1}P^{\top}\big{)}$
;

14 // $g\in\mathbb{R}^{m}$

15

16for  _$i=1,2,\cdots,t$_ do

17 $\delta_{y}\leftarrow-N\cdot g$ ;

18 // $\delta_{y}\in\mathbb{R}^{m}$

19

20 $y^{\operatorname{new}}\leftarrow y+\delta_{y}$ ;

21 // $y^{\operatorname{new}}\in\mathbb{R}^{m}$

22

23 $s^{\operatorname{new}}\leftarrow c-A^{\top}y^{\operatorname{new}}$ ;

24 // $s^{\operatorname{new}}\in\mathbb{R}^{U}$

25

26 $\eta^{\operatorname{new}}\leftarrow\eta\cdot(1+\alpha)$;

27

28 $S^{\operatorname{new}}\leftarrow
P^{\top}\operatorname{diag}(s^{\operatorname{new}})P$ ;

29 // $S^{\operatorname{new}}\in\mathbb{R}^{L\times L}$

30

31
$\widetilde{S}^{\operatorname{new}},V_{1},V_{2}\leftarrow\textsc{LowRankUpdate}(S^{\operatorname{new}},\widetilde{S})$;

32

33 // Lemma 9.1, $\widetilde{S}^{\operatorname{new}}\in\mathbb{R}^{L\times
L}$, $V_{1},V_{2}\in\mathbb{R}^{L\times r_{i}}$ or $V_{1}=V_{2}=\mathtt{null}$

34

35 if  _$V_{1}=V_{2}=\mathtt{null}$_ then

36

37 $T^{\operatorname{new}}\leftarrow(\widetilde{S}^{\operatorname{new}})^{-1}$
;

38 // $T^{\operatorname{new}}\in\mathbb{R}^{L\times L}$

39

40
$N^{\operatorname{new}}\leftarrow\big{(}A\cdot(PT^{\operatorname{new}}P^{\top})^{\circ
2}\cdot A^{\top}\big{)}^{-1}$ ;

41 // $N^{\operatorname{new}}\in\mathbb{R}^{m\times m}$

42

43 else

44

45
$T^{\operatorname{new}},N^{\operatorname{new}}\leftarrow\textsc{UpdateHessianInv}(T,N,V_{1},V_{2})$;

46

47 // Lemma 8.1, $T^{\operatorname{new}}\in\mathbb{R}^{L\times
L},N^{\operatorname{new}}\in\mathbb{R}^{m\times m}$

48

49

50 $g^{\operatorname{new}}\leftarrow-\eta^{\operatorname{new}}\cdot
b+A\cdot\operatorname{diag}\Big{(}P\big{(}P^{\top}\operatorname{diag}(s^{\operatorname{new}})P\big{)}^{-1}P^{\top}\Big{)}$
;

51 // $g^{\operatorname{new}}\in\mathbb{R}^{m}$

52

53
$(\eta,y,s,\widetilde{S},T,N,g)\leftarrow(\eta^{\operatorname{new}},y^{\operatorname{new}},s^{\operatorname{new}},\widetilde{S}^{\operatorname{new}},T^{\operatorname{new}},N^{\operatorname{new}},g^{\operatorname{new}})$;

54

return _$(y,s)$_

Algorithm 1 Main SOS algorithm.

##  8 Updating Hessian inverse efficiently

In this section we prove how to update the Hessian inverse efficiently. We
present the algorithm UpdateHessianInv in Algorithm 2.

Input : $T\in\mathbb{R}^{L\times L}$, $N\in\mathbb{R}^{m\times m}$,
$V_{1},V_{2}\in\mathbb{R}^{L\times r}$

Output : $T^{\operatorname{new}}\in\mathbb{R}^{L\times L}$,
$N^{\operatorname{new}}\in\mathbb{R}^{m\times m}$

1

2// Step 1

3

4$\overline{V_{1}}\leftarrow-TV_{1}\cdot(I+V_{2}^{\top}TV_{1})^{-1}$ ;

5 // $\overline{V_{1}}\in\mathbb{R}^{L\times r}$

6

7$\overline{V_{2}}\leftarrow TV_{2}$ ;

8 // $\overline{V_{2}}\in\mathbb{R}^{L\times r}$

9

10$T^{\operatorname{new}}\leftarrow
T+\overline{V_{1}}\cdot\overline{V_{2}}^{\top}$ ;

11 // $T^{\operatorname{new}}\in\mathbb{R}^{L\times L}$

12

13// Step 2

14

15$Y^{\prime}\leftarrow[2PT,P\overline{V_{1}}]$ ;

16 // $Y^{\prime}\in\mathbb{R}^{U\times(L+r)}$

17

18$Z^{\prime}\leftarrow[P,P\overline{V_{2}}]$ ;

19 // $Z^{\prime}\in\mathbb{R}^{U\times(L+r)}$

20

21$Y\leftarrow[\operatorname{diag}(u_{1})Y^{\prime},\cdots,\operatorname{diag}(u_{r})Y^{\prime}]$,
$u_{i}$ is the $i$-th column of $P\overline{V_{1}}$ ;

22 // $Y\in\mathbb{R}^{U\times(L+r)r}$

23

24$Z\leftarrow[\operatorname{diag}(v_{1})Z^{\prime},\cdots,\operatorname{diag}(v_{r})Z^{\prime}]$,
$v_{i}$ is the $i$-th column of $P\overline{V_{2}}$ ;

25 // $Z\in\mathbb{R}^{U\times(L+r)r}$

26

27// Step 3

28

29$N^{\operatorname{new}}\leftarrow
N-N\cdot(AY)\cdot\big{(}I+(AZ)^{\top}N(AY)\big{)}^{-1}\cdot(AZ)^{\top}\cdot N$
;

30 // $N^{\operatorname{new}}\in\mathbb{R}^{m\times m}$

31

return _$T^{\operatorname{new}},N^{\operatorname{new}}$_

Algorithm 2 UpdateHessianInv

######  Lemma 8.1 (Hessian inverse update).

In the algorithm UpdateHessianInv (Algorithm 2), the inputs are the maintained
matrices $T,N$ and the updates $V_{1},V_{2}\in\mathbb{R}^{L\times r}$ where
$r$ satisfies $Lr\leq U$. The inputs satisfy that for some
$\widetilde{S}\in\mathbb{S}^{L\times L}$,

$\displaystyle T=\widetilde{S}^{-1}\in\mathbb{R}^{L\times L},$ $\displaystyle
N=\big{(}A\cdot(P\widetilde{S}^{-1}P^{\top})^{\circ 2}\cdot
A^{\top}\big{)}^{-1}\in\mathbb{R}^{m\times m},$

Let $\widetilde{S}^{\operatorname{new}}=\widetilde{S}+V_{1}V_{2}^{\top}$. The
algorithm outputs two matrices $T^{\operatorname{new}}$,
$N^{\operatorname{new}}$ such that

$\displaystyle
T^{\operatorname{new}}=(\widetilde{S}^{\operatorname{new}})^{-1}\in\mathbb{R}^{L\times
L},$ $\displaystyle
N^{\operatorname{new}}=\big{(}A\cdot(P(\widetilde{S}^{\operatorname{new}})^{-1}P^{\top})^{\circ
2}\cdot A^{\top}\big{)}^{-1}\in\mathbb{R}^{m\times m}.$

Furthermore, the algorithm takes $O({\cal T}_{\mathrm{mat}}(U,U,Lr))$ time.

###### Proof.

We first prove the correctness by analyzing each step of the algorithm.

Step 1. Compute $\overline{V_{1}},\overline{V_{2}}\in\mathbb{R}^{L\times r}$
and $T^{\operatorname{new}}\in\mathbb{R}^{L\times L}$.

$\displaystyle T^{\operatorname{new}}=$ $\displaystyle\leavevmode\nobreak\
T+\overline{V_{1}}\cdot\overline{V_{2}}^{\top}$ $\displaystyle=$
$\displaystyle\leavevmode\nobreak\
T-TV_{1}\cdot(I+V_{2}^{\top}TV_{1})^{-1}\cdot V_{2}^{\top}T^{\top}$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
(\widetilde{S}+V_{1}V_{2}^{\top})^{-1}=(\widetilde{S}^{\operatorname{new}})^{-1},$

where the first two steps follow from algorithm description, the third step
follows from the Woodbury identity (Fact 5.1) and $T=\widetilde{S}^{-1}$.

Thus $T^{\operatorname{new}}$ satisfies the requirement of the output.

Step 2. Compute $Y^{\prime},Z^{\prime}\in\mathbb{R}^{U\times(L+r)}$ and
$Y,Z\in\mathbb{R}^{U\times(L+r)r}$

We prove that $Y$ and $Z$ satisfy $(PTP^{\top})^{\circ 2}+Y\cdot
Z^{\top}=(PT^{\operatorname{new}}P^{\top})^{\circ 2}$:

$\displaystyle(PTP^{\top})^{\circ 2}+Y\cdot Z^{\top}=$
$\displaystyle\leavevmode\nobreak\ (PTP^{\top})^{\circ
2}+\sum_{i=1}^{r}\operatorname{diag}(u_{i})\cdot\big{(}Y^{\prime}\cdot(Z^{\prime})^{\top}\big{)}\cdot\operatorname{diag}(v_{i})$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\ (PTP^{\top})^{\circ
2}+\big{(}Y^{\prime}\cdot(Z^{\prime})^{\top}\big{)}\circ\big{(}\sum_{i=1}^{r}u_{i}\cdot
v_{i}^{\top}\big{)}$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\
(PTP^{\top})^{\circ
2}+\big{(}2PTP^{\top}+(P\overline{V_{1}})\cdot(P\overline{V_{2}})^{\top}\big{)}\circ\big{(}(P\overline{V_{1}})\cdot(P\overline{V_{2}})^{\top}\big{)}$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\big{(}PTP^{\top}+(P\overline{V_{1}})\cdot(P\overline{V_{2}})^{\top}\big{)}^{\circ
2}$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\
(PT^{\operatorname{new}}P^{\top})^{\circ 2},$

where the first step follows from the algorithm description of $Y$ and $Z$,
the second step follows from Part 2 of Fact 5.3 that
$\operatorname{diag}(x)\cdot A\cdot\operatorname{diag}(y)=A\circ(x\cdot
y^{\top})$, the third step follows from
$Y^{\prime}\cdot(Z^{\prime})^{\top}=2PTP^{\top}+(P\overline{V_{1}})\cdot(P\overline{V_{2}})$
and $(P\overline{V_{1}})\cdot(P\overline{V_{2}})=\sum_{i=1}^{r}u_{i}\cdot
v_{i}^{\top}$ (see algorithm description of $Y^{\prime}$ and $Z^{\prime}$),
the last step follows from
$T^{\operatorname{new}}=T+\overline{V_{1}}\cdot\overline{V_{2}}^{\top}$.

Step 3. Compute $N^{\operatorname{new}}\in\mathbb{R}^{m\times m}$.

$\displaystyle N^{\operatorname{new}}=$ $\displaystyle\leavevmode\nobreak\
N-N\cdot(AY)\cdot\big{(}I+(AZ)^{\top}N(AY)\big{)}^{-1}\cdot(AZ)^{\top}\cdot N$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\big{(}A\cdot(PTP^{\top})^{\circ 2}\cdot
A^{\top}+(AY)\cdot(AZ)^{\top}\big{)}^{-1}$ $\displaystyle=$
$\displaystyle\leavevmode\nobreak\
\big{(}A\cdot(PT^{\operatorname{new}}P^{\top})^{\circ 2}\cdot
A^{\top}\big{)}^{-1},$

where the first step follows from the algorithm description of
$N^{\operatorname{new}}$, the second step follows from
$N=\big{(}A\cdot(PTP^{\top})^{\circ 2}\cdot A^{\top}\big{)}^{-1}$ and the
Woodbury identity (Fact 5.1), and the last step follows from
$(PT^{\operatorname{new}}P^{\top})^{\circ 2}=(PTP^{\top})^{\circ 2}+Y\cdot
Z^{\top}$.

Thus $N^{\operatorname{new}}$ satisfies the requirement of the output.

Time complexity. It is easy to see that the most time-consuming step is to
compute $N^{\operatorname{new}}$ on Line 2, and in total this step takes
$O({\cal T}_{\mathrm{mat}}(m,U,Lr)+{\cal
T}_{\mathrm{mat}}(m,m,Lr)+(Lr)^{\omega})$ time.

Since $Lr\leq U$ and $m\leq U$, overall this algorithm takes at most $O({\cal
T}_{\mathrm{mat}}(U,U,Lr))$ time. ∎

##  9 Correctness

###  9.1 Standard results from IPM theory

We use the following two results of the barrier method that hold for any cone
with a barrier function. The proofs are standard, (see e.g.,
[renegarMathematicalViewInteriorPoint2001, Section 2.4]), and for completeness
we include a proof in the Appendix. [Invariance of Newton step,
[renegarMathematicalViewInteriorPoint2001]] Consider the following
optimization problem: $\min-b^{\top}y$ s.t. $y\in\overline{D}_{F}$, where
$F:\mathbb{R}^{m}\to\mathbb{R}_{+}$ is a barrier function with barrier
parameter $\nu_{F}$, $D_{F}\subseteq\mathbb{R}^{m}$ is the domain of $F$, and
$\overline{D}_{F}$ is the closure of $D_{F}$. For any $\eta\geq 1$, define
$F_{\eta}(y)=-\eta b^{\top}y+F(y)$. Let $g_{\eta}(y)\in\mathbb{R}^{m}$ and
$H(y)\in\mathbb{R}^{m\times m}$ denote the gradient and the Hessian of
$F_{\eta}$ at $y$.

Let $0<\epsilon_{N}\leq 0.05$ be a parameter. If a feasible solution $y\in
D_{F}$, a parameter $\eta>0$, and a positive definite matrix
$\widetilde{H}\in\mathbb{S}^{n\times n}$ satisfy the following:

$\|g_{\eta}(y)\|_{H(y)^{-1}}\leq\epsilon_{N},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \widetilde{H}\approx_{0.02}H(y).$

Then
$\eta^{\operatorname{new}}=\eta\cdot(1+\frac{\epsilon_{N}}{20\sqrt{\nu_{F}}})$,
$y^{\operatorname{new}}=y+\delta_{y}$ where
$\delta_{y}=-\widetilde{H}^{-1}g_{\eta^{\operatorname{new}}}(y)$ satisfy
$y^{\operatorname{new}}\in D_{F}$ and

$\|\delta_{y}\|_{H(y)}\leq 2\epsilon_{N},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\|g_{\eta^{\operatorname{new}}}(y^{\operatorname{new}})\|_{H(y^{\operatorname{new}})^{-1}}\leq\epsilon_{N}.$

###### Proof.

For clarity in the proof we use $\eta_{1}:=\eta$,
$\eta_{2}:=\eta^{\operatorname{new}}$, $y_{1}:=y$, and
$\widetilde{y}_{2}:=y^{\operatorname{new}}=y_{1}-\widetilde{H}^{-1}g_{\eta_{2}}(y_{1})$.
We also define $\alpha_{H}:=e^{0.02}\leq 1.03$. Note that
$\alpha_{H}^{-1}H(y_{1})\preceq\widetilde{H}\preceq\alpha_{H}H(y_{1})$.

We first introduce some notation:

$\displaystyle n_{\eta_{1}}(y_{1}):=$ $\displaystyle\leavevmode\nobreak\
-H(y_{1})^{-1}g_{\eta_{1}}(y_{1}),$ $\displaystyle n_{\eta_{2}}(y_{1}):=$
$\displaystyle\leavevmode\nobreak\ -H(y_{1})^{-1}g_{\eta_{2}}(y_{1}),$
$\displaystyle\widetilde{n}_{\eta_{2}}(y_{1}):=$
$\displaystyle\leavevmode\nobreak\
-\widetilde{H}(y_{1})^{-1}g_{\eta_{2}}(y_{1}),$ $\displaystyle
n_{\eta_{2}}(\widetilde{y}_{2}):=$ $\displaystyle\leavevmode\nobreak\
-H(\widetilde{y}_{2})^{-1}g_{\eta_{2}}(\widetilde{y}_{2})$

Note that $\widetilde{y}_{2}=y_{1}+\widetilde{n}_{\eta_{2}}(y_{1})$. Recall
that we define $\|x\|_{H(y)}^{2}=x^{\top}H(y)x,\forall x\in\mathbb{R}^{m}$ to
be the local norm at $y$. The local norm induces a matrix norm
$\|M\|_{H(y)}=\max_{x}\frac{\|Mx\|_{H(y)}}{\|x\|_{H(y)}}=\|H(Y)^{1/2}MH(Y)^{-1/2}\|_{2},\forall
M\in\mathbb{R}^{m\times m}$.

Step 1 (Bound $\|n_{\eta_{2}}(y_{1})\|_{H(y_{1})}$). We have

$\displaystyle n_{\eta_{2}}(y_{1})=$ $\displaystyle\leavevmode\nobreak\
-H(y_{1})^{-1}g_{\eta_{2}}(y_{1})$ $\displaystyle=$
$\displaystyle\leavevmode\nobreak\ -H(y_{1})^{-1}\cdot(-\eta_{2}\cdot
b+g(y_{1}))$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\frac{\eta_{2}}{\eta_{1}}\cdot
n_{\eta_{1}}(y_{1})+(\frac{\eta_{2}}{\eta_{1}}-1)\cdot H(y_{1})^{-1}g(y_{1}).$

Thus we have

$\displaystyle\|n_{\eta_{2}}(y_{1})\|_{H(y_{1})}\leq$
$\displaystyle\leavevmode\nobreak\
\frac{\eta_{2}}{\eta_{1}}\cdot\|n_{\eta_{1}}(y_{1})\|_{H(y_{1})}+|\frac{\eta_{2}}{\eta_{1}}-1|\cdot\|H(y_{1})^{-1}g(y_{1})\|_{H(y_{1})}$
(20) $\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
\frac{\eta_{2}}{\eta_{1}}\cdot\|n_{\eta_{1}}(y_{1})\|_{H(y_{1})}+|\frac{\eta_{2}}{\eta_{1}}-1|\cdot\sqrt{\nu_{F}}$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\ 1.1\epsilon_{N},$

where the second step follows from the definition of the barrier parameter
(Section 2.3.1 of [renegarMathematicalViewInteriorPoint2001]), and the third
step follows from
$\eta_{2}=\eta_{1}\cdot(1+\frac{\epsilon_{N}}{20\sqrt{\nu_{F}}})$ and
$\|n_{\eta_{1}}(y_{1})\|_{H(y_{1})}\leq\epsilon_{N}$.

Step 2 (Bound $\|\widetilde{n}_{\eta_{2}}(y_{1})\|_{H(y_{1})}$). We have

$\displaystyle\|\widetilde{n}_{\eta_{2}}(y_{1})\|_{H(y_{1})}^{2}=$
$\displaystyle\leavevmode\nobreak\
g_{\eta_{2}}(y_{1})^{\top}\widetilde{H}^{-1}H(y_{1})\widetilde{H}^{-1}g_{\eta_{2}}(y_{1})$
(21) $\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\ \alpha_{H}\cdot
g_{\eta_{2}}(y_{1})^{\top}\widetilde{H}^{-1}g_{\eta_{2}}(y_{1})$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\ \alpha_{H}^{2}\cdot
g_{\eta_{2}}(y_{1})^{\top}H(y_{1})^{-1}g_{\eta_{2}}(y_{1})$ $\displaystyle=$
$\displaystyle\leavevmode\nobreak\
\alpha_{H}^{2}\cdot\|n_{\eta_{2}}(y_{1})\|_{H(y_{1})}^{2}$ $\displaystyle\leq$
$\displaystyle\leavevmode\nobreak\ 4\epsilon_{N}^{2},$

where the second step follows from $H(y_{1})\preceq\alpha_{H}\widetilde{H}$,
the third step follows from
$\widetilde{H}^{-1}\preceq\alpha_{H}H(y_{1})^{-1}$, and the fifth step follows
from $\alpha_{H}\leq 1.03$ and $\|n_{\eta_{2}}(y_{1})\|_{H(y_{1})}\leq
1.1\epsilon_{N}$ (Eq. (20)).

Note that $\widetilde{y}_{2}\in D_{F}$ directly follows from this bound and
the definition of self-concordant functions. This proves the first inequality
that of $\|\delta_{y}\|_{H(y)}\leq 2\epsilon_{N}$.

Step 3 (Bound $\|H(\widetilde{y}_{2})^{-1}H(y_{1})\|_{H(y_{1})}$). We have

$\displaystyle\|H(\widetilde{y}_{2})^{-1}H(y_{1})\|_{H(y_{1})}\leq$
$\displaystyle\leavevmode\nobreak\
\frac{1}{(1-\|\widetilde{y}_{2}-y_{1}\|_{H(y_{1})})^{2}}$ (22)
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\frac{1}{(1-\|\widetilde{n}_{\eta_{2}}(y_{1})\|_{H(y_{1})})^{2}}$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\ 1/(0.9)^{2}$

where the first step follows from Theorem 2.2.1 of
[renegarMathematicalViewInteriorPoint2001], the second step follows from
$\widetilde{y}_{2}=y_{1}+\widetilde{n}_{\eta_{2}}(y_{1})$, the third step
follows from $\|\widetilde{n}_{\eta_{2}}(y_{1})\|_{H(y_{1})}\leq
2\epsilon_{N}$ (Eq. (21)) and $\epsilon_{N}\leq 0.05$.

Step 4 (Bound $\|H(y_{1})^{-1}g_{\eta_{2}}(\widetilde{y}_{2})\|_{H(y_{1})}$).
We have

$\displaystyle\|H(y_{1})^{-1}g_{\eta_{2}}(\widetilde{y}_{2})\|_{H(y_{1})}\leq$
$\displaystyle\leavevmode\nobreak\
\underbrace{\|H(y_{1})^{-1}\cdot\big{(}g_{\eta_{2}}(\widetilde{y}_{2})-g_{\eta_{2}}(y_{1})-H(y_{1})\widetilde{n}_{\eta_{2}}(y_{1})\big{)}\|_{H(y_{1})}}_{a_{1}}$
(23) $\displaystyle\leavevmode\nobreak\
+\underbrace{\|H(y_{1})^{-1}g_{\eta_{2}}(y_{1})+\widetilde{n}_{\eta_{2}}(y_{1})\|_{H(y_{1})}}_{a_{2}}.$

Next we bound these two terms separately. For the first term, we have

$\displaystyle a_{1}=$ $\displaystyle\leavevmode\nobreak\
\|H(y_{1})^{-1}\cdot\int_{0}^{1}[H(y_{1}+t\widetilde{n}_{\eta_{2}}(y_{1}))-H(y_{1})]\cdot\widetilde{n}_{\eta_{2}}(y_{1})\leavevmode\nobreak\
\mathrm{d}t\|_{H(y_{1})}$ (24) $\displaystyle\leq$
$\displaystyle\leavevmode\nobreak\
\|\widetilde{n}_{\eta_{2}}(y_{1})\|_{H(y_{1})}\cdot\int_{0}^{1}\|I-H(y_{1})^{-1}\cdot
H(y_{1}+t\widetilde{n}_{\eta_{2}}(y_{1}))\|_{H(y_{1})}\leavevmode\nobreak\
\mathrm{d}t$ $\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
\|\widetilde{n}_{\eta_{2}}(y_{1})\|_{H(y_{1})}\cdot\int_{0}^{1}\Big{(}\frac{1}{(1-t\|\widetilde{n}_{\eta_{2}}(y_{1})\|_{H(y_{1})})^{2}}-1\Big{)}\mathrm{d}t$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\frac{\|\widetilde{n}_{\eta_{2}}(y_{1})\|_{H(y_{1})}^{2}}{1-\|\widetilde{n}_{\eta_{2}}(y_{1})\|_{H(y_{1})}}$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\ 5\epsilon_{N}^{2},$

where the first step follows from
$\widetilde{y}_{2}-y_{1}=\widetilde{n}_{\eta_{2}}(y_{1})$ and Proposition
1.5.7 of [renegarMathematicalViewInteriorPoint2001], the third step follows
from Theorem 2.2.1 of [renegarMathematicalViewInteriorPoint2001], and the
fifth step follows from $\|\widetilde{n}_{\eta_{2}}(y_{1})\|_{H(y_{1})}\leq
2\epsilon_{N}$ (Eq. (21)) and $\epsilon_{N}\leq 0.05$.

For the second term, we have

$\displaystyle a_{2}^{2}=$ $\displaystyle\leavevmode\nobreak\
\Big{|}g_{\eta_{2}}(y_{1})^{\top}\cdot\Big{(}H(y_{1})^{-1}-\widetilde{H}^{-1}H(y_{1})\widetilde{H}^{-1}\Big{)}\cdot
g_{\eta_{2}}(y_{1})\Big{|}$ (25) $\displaystyle\leq$
$\displaystyle\leavevmode\nobreak\
\max\\{\alpha_{H}^{2}-1,1-\alpha_{H}^{-2}\\}\cdot
g_{\eta_{2}}(y_{1})^{\top}\cdot H(y_{1})^{-1}\cdot g_{\eta_{2}}(y_{1})$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\ 0.09\epsilon_{N}^{2}$

where the first step follows from
$\widetilde{n}_{\eta_{2}}(y_{1})=-\widetilde{H}^{-1}\cdot
g_{\eta_{2}}(y_{1})$, the second step follows from
$\alpha_{H}^{-1}H(y_{1})\preceq\widetilde{H}\preceq\alpha_{H}H(y_{1})$, the
third step follows from $\|n_{\eta_{2}}(y_{1})\|_{H(y_{1})}\leq
1.1\epsilon_{N}$ (Eq. (20)) and $\alpha_{N}\leq 1.03$.

Plug in Eq. (24) and (25) into Eq. (23), we have

$\|H(y_{1})^{-1}g_{\eta_{2}}(\widetilde{y}_{2})\|_{H(y_{1})}\leq
5\epsilon_{N}^{2}+0.3\epsilon_{N}\leq 0.6\epsilon_{N},$ (26)

where the last step follows from $\epsilon_{N}\leq 0.05$.

Step 5 (Bound $\|n_{\eta_{2}}(\widetilde{y}_{2})\|_{H(\widetilde{y}_{2})}$).
We are finally ready to bound
$\|n_{\eta_{2}}(\widetilde{y}_{2})\|_{H(\widetilde{y}_{2})}$. We have

$\displaystyle\|n_{\eta_{2}}(\widetilde{y}_{2})\|_{H(\widetilde{y}_{2})}^{2}=$
$\displaystyle\leavevmode\nobreak\
g_{\eta_{2}}(\widetilde{y}_{2})^{\top}H(\widetilde{y}_{2})^{-1}g_{\eta_{2}}(\widetilde{y}_{2})$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
\|H(\widetilde{y}_{2})^{-1}H(y_{1})\|_{H(y_{1})}^{2}\cdot\|H(y_{1})^{-1}g_{\eta_{2}}(\widetilde{y}_{2})\|_{H(y_{1})}^{2}$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
1/(0.9)^{4}\cdot(0.6\epsilon_{N})^{2}\leq\epsilon_{N}^{2},$

where the third step follows from Eq. (22) and Eq. (26). This proves the
second inequality that
$\|g_{\eta^{\operatorname{new}}}(y^{\operatorname{new}})\|_{H(y^{\operatorname{new}})^{-1}}\leq\epsilon_{N}$.
∎

[Approximate optimality, [renegarMathematicalViewInteriorPoint2001]] Consider
the following optimization problem: $\min-b^{\top}y$ s.t.
$y\in\overline{D}_{F}$, where $F:\mathbb{R}^{m}\to\mathbb{R}_{+}$ is a barrier
function with barrier parameter $\nu_{F}$, $D_{F}\subseteq\mathbb{R}^{m}$ is
the domain of $F$, and $\overline{D}_{F}$ is the closure of $D_{F}$. Let
$\mathrm{OPT}$ be the optimal objective value of this optimization problem.
For any $\eta\geq 1$, define $F_{\eta}(y)=-\eta b^{\top}y+F(y)$. Let
$g_{\eta}(y)\in\mathbb{R}^{m}$ and $H(y)\in\mathbb{R}^{m\times m}$ denote the
gradient and the Hessian of $F_{\eta}$ at $y$.

Let $0<\epsilon_{N}\leq 0.05$. If a feasible solution $y\in D_{F}$ satisfies
$\|g_{\eta}(y)\|_{H(y)^{-1}}\leq\epsilon_{N}$, then we have
$-b^{\top}y\leq\mathrm{OPT}+\frac{\nu_{F}}{\eta}\cdot(1+2\epsilon_{N})$.

###### Proof.

We use $g(y)\in\mathbb{R}^{m}$ to denote the gradient of $F$. Note that the
Hessian of $F$ and $F_{\eta}$ are the same. Let $z(\eta)\in D_{F}$ be the
minimizer of $F_{\eta}(y)=-\eta\langle b,y\rangle+F(y)$. Let
$y^{*}\in\overline{D}_{F}$ denote the optimal solution of $\min-\langle
b,y\rangle$ s.t. $y\in\overline{D}_{F}$.

Step 1. We have

$\displaystyle-\langle b,z(\eta)\rangle-\mathrm{OPT}=$
$\displaystyle\leavevmode\nobreak\ -\langle b,z(\eta)\rangle+\langle
b,y^{*}\rangle$ (27) $\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\frac{1}{\eta}\langle\eta b,y^{*}-z(\eta)\rangle$ $\displaystyle=$
$\displaystyle\leavevmode\nobreak\ \frac{1}{\eta}\langle
g(z(\eta)),y^{*}-z(\eta)\rangle<\frac{1}{\eta}\nu_{F},$

where the second step follows from $g_{\eta}(z(\eta))=-\eta b+g(z(\eta))=0$
since $z(\eta)$ is the minimizer of $F_{\eta}$, the third step follows from
Theorem 2.3.3 of [renegarMathematicalViewInteriorPoint2001] that $\langle
g(x),y-x\rangle<\nu_{F}$ for any $x,y\in D_{F}$.

Step 2. By the definition of self-concordant functions (Section 2.2.1 of
[renegarMathematicalViewInteriorPoint2001]), for any $z\in D_{F}$, we have
that $z+t\cdot H(z)^{-1}b\in D_{F}$ if
$t\in[0,\frac{1}{\|H(z)^{-1}b\|_{H(z)}})$. So we have

$-\langle b,z+t\cdot H(z)^{-1}b\rangle\geq\mathrm{OPT},$

which implies that

$\langle b,t\cdot H(z)^{-1}b\rangle\leq-\langle b,z\rangle-\mathrm{OPT},$

and as $t$ approaches
$\frac{1}{\|H(z)^{-1}b\|_{H(z)}}=\frac{1}{\|b\|_{H(z)^{-1}}}$, we have

$\|b\|_{H(z)^{-1}}\leq-\langle b,z\rangle-\mathrm{OPT}.$ (28)

Step 3. Hence for all $y,z\in D_{F}$, we have

$\displaystyle\frac{-\langle b,y\rangle-\mathrm{OPT}}{-\langle
b,z\rangle-\mathrm{OPT}}=$ $\displaystyle\leavevmode\nobreak\ 1+\frac{-\langle
b,y-z\rangle}{-\langle b,z\rangle-\mathrm{OPT}}$ (29) $\displaystyle\leq$
$\displaystyle\leavevmode\nobreak\
1+\frac{\|b\|_{H(z)^{-1}}\cdot\|y-z\|_{H(z)}}{-\langle
b,z\rangle-\mathrm{OPT}}$ $\displaystyle\leq$
$\displaystyle\leavevmode\nobreak\ 1+\|y-z\|_{H(z)},$

where the last step follows from Eq. (28).

Step 4. Plugging in $z=z(\eta)$ to Eq. (29), we have

$\displaystyle-\langle b,y\rangle\leq$ $\displaystyle\leavevmode\nobreak\
\mathrm{OPT}+(-\langle
b,z(\eta)\rangle-\mathrm{OPT})\cdot(1+\|y-z(\eta))\|_{H(z(\eta))})$ (30)
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
\mathrm{OPT}+\frac{\nu_{F}}{\eta}\cdot(1+\|y-z(\eta))\|_{H(z(\eta))}),$

where the second step follows from Eq. (27).

Step 5. Finally, for simplicity denote the Newton step of $F_{\eta}$ at $y$ as
$n(y):=H(y)^{-1}g_{\eta}(y)\in\mathbb{R}^{m}$. We have
$\|n(y)\|_{H(y)}\leq\epsilon_{N}$ where $\epsilon_{N}\leq 0.05$, and using
Theorem 2.2.5 of [renegarMathematicalViewInteriorPoint2001], we have

$\|y-z(\eta)\|_{H(y)}\leq\|n(y)\|_{H(y)}+\frac{3\|n(y)\|_{H(y)}^{2}}{(1-\|n(y)\|_{H(y)})^{3}}\leq\epsilon_{N}\cdot(1+\frac{3\epsilon_{N}}{(1-\epsilon_{N})^{3}})\leq
1.5\epsilon_{N}.$

Then by the definition of self-concordance (Section 2.2.1 of
[renegarMathematicalViewInteriorPoint2001]), we have

$\|y-z(\eta)\|_{H(z(\eta))}\leq\frac{\|y-z(\eta)\|_{H(y)}}{1-\|y-z(\eta)\|_{H(y)}}\leq\frac{1.5\epsilon_{N}}{1-1.5\epsilon_{N}}\leq
2\epsilon_{N}.$

Thus, from Eq. (30) we have

$-\langle
b,y\rangle\leq\mathrm{OPT}+\frac{\nu_{F}}{\eta}\cdot(1+2\epsilon_{N}).\qed$

###  9.2 Low rank update

We use the following low rank update procedure of [jiang2020faster] and
[hjst21], which we modify by using a cutoff when $r\geq U/L$. The proof of the
following lemma can be found in [hjst21, Theorem 10.8].

Parameters : A real number $\epsilon_{S}<0.01$.

Input : New exact matrix $S^{\operatorname{new}}\in\mathbb{R}^{L\times L}$,
old approximate matrix $\widetilde{S}\in\mathbb{R}^{L\times L}$,

Output : New approximate matrix
$\widetilde{S}^{\operatorname{new}}\in\mathbb{R}^{L\times L}$, update matrices
$V_{1},V_{2}\in\mathbb{R}^{L\times r}$.

1

2$Z_{\operatorname{mid}}\leftarrow(S^{\operatorname{new}})^{-1/2}\widetilde{S}(S^{\operatorname{new}})^{-1/2}-I$

3Compute spectral decomposition
$Z_{\operatorname{mid}}=X\operatorname{diag}(\lambda)X^{\top}$

4Let $\pi:[L]\to[L]$ be a sorting permutation such that
$|\lambda_{\pi(i)}|\geq|\lambda_{\pi(i+1)}|$.

5if  _$|\lambda_{\pi(1)}|\leq\epsilon_{S}$_ then

6 $\widetilde{S}_{\operatorname{new}}\leftarrow\widetilde{S}$;

7

8 return _$(\widetilde{S}^{\operatorname{new}},0,0)$_

9else

10 $r\leftarrow 1$;

11

12 while  _$2r\leq U/L$ _and_ $(|\lambda_{\pi(2r)}|>\epsilon_{S}$ _or_
$|\lambda_{\pi(2r)}|>(1-1/\log L)|\lambda_{\pi(r)}|)$_ do

13 $r\leftarrow r+1$;

14

15 $r\leftarrow 2r$;

16

17 if  _$r\geq U/L$_ then

18 $\widetilde{S}^{\operatorname{new}}\leftarrow S^{\operatorname{new}}$ ;

19 // Here we deviate from [jiang2020faster]

20

21 return _$(\widetilde{S}^{\operatorname{new}},\mathtt{null},\mathtt{null})$_

22 else

23 $\lambda^{\operatorname{new}}_{\pi(i)}\leftarrow\begin{cases}0&\text{if
}i=1,2,\ldots,r\\\ \lambda_{\pi(i)}&\text{else }\end{cases}$

24 $\Omega\leftarrow\mathrm{supp}(\lambda^{\operatorname{new}}-\lambda)$ ;

25 // $|\Omega|=r$

26

27 $V_{1}\leftarrow((S^{\operatorname{new}})^{1/2}\cdot
X\cdot\operatorname{diag}(\lambda^{\operatorname{new}}-\lambda))_{:,\Omega}$ ;

28 // $V_{1}\in\mathbb{R}^{L\times r}$

29

30 $V_{2}\leftarrow((S^{\operatorname{new}})^{1/2}\cdot X)_{:,\Omega}$ ;

31 // $V_{2}\in\mathbb{R}^{L\times r}$

32

33
$\widetilde{S}^{\operatorname{new}}\leftarrow\widetilde{S}+(S^{\operatorname{new}})^{1/2}X\operatorname{diag}(\lambda^{\operatorname{new}}-\lambda)X^{\top}(S^{\operatorname{new}})^{1/2}$
;

34 //
$\widetilde{S}^{\operatorname{new}}=\widetilde{S}+V_{1}V_{2}^{\top}\in\mathbb{R}^{L\times
L}$

35

36 return _$(\widetilde{S}^{\operatorname{new}},V_{1},V_{2})$_ ;

37

38

Algorithm 3 LowRankUpdate of [jiang2020faster]

######  Lemma 9.1 (Low rank update).

The algorithm LowRankUpdate (Algorithm 3) has the following properties:

  1. (i)

The output matrix
$\widetilde{S}^{\operatorname{new}}=\widetilde{S}+V_{1}V_{2}^{\top}$ is a
spectral approximation of the input matrix:

$\widetilde{S}^{\operatorname{new}}\approx_{\epsilon_{S}}S^{\operatorname{new}}.$

  2. (ii)

Consider $t$ iterations of LowRankUpdate. Initially
$\widetilde{S}^{(0)}=S^{(0)}$, and we use $(S^{(i)},\widetilde{S}^{(i-1)})$
and $(\widetilde{S}^{(i)},V_{1}^{(i)},V_{2}^{(i)})$ to denote the input and
the output of the $i$-th iteration. We define the rank $r_{i}$ to be the rank
of $V_{1}^{(i)}$ if $V_{1}^{(i)}\neq\mathtt{null}$, and otherwise we define
$r_{i}=U/L$.

If the input exact matrices
$S^{(0)},S^{(1)},\cdots,S^{(t)}\in\mathbb{R}^{L\times L}$ satisfy

$\|(S^{(i-1)})^{-1/2}S^{(i)}(S^{(i-1)})^{-1/2}-I\|_{F}\leq
0.02,\leavevmode\nobreak\ \forall i\in[t].$ (31)

Then for any non-increasing sequence $g\in\mathbb{R}^{L}_{+}$, the ranks
$r_{i}$ satisfy

$\sum_{i=1}^{t}r_{i}\cdot g_{r_{i}}\leq O(t\cdot\|g\|_{2}\cdot\log L).$

Furthermore, the algorithm LowRankUpdate takes $O(L^{\omega})$ time.

###  9.3 Slowly moving guarantee

In SOS optimization, the matrix $S=P^{\top}\operatorname{diag}(s)P$
corresponds to the slack matrix of the SDP. The following lemma proves similar
to SDP, in SOS the matrix $S$ is changing slowly. The proof is deferred to the
appendix. Using this lemma we will prove that the requirement Eq. (31) of
Lemma 9.1 is satisfied, which means we can approximate the change to the slack
by a low-rank matrix.

######  Lemma 9.2 (Slowly moving guarantee).

Let $c\in\mathbb{R}^{U}$ and $A\in\mathbb{R}^{m\times U}$ be the input to the
optimization problem. Let $P\in\mathbb{R}^{U\times L}$ be the matrix of the
interpolant basis.

For any $y\in\mathbb{R}^{m}$ and
$y^{\operatorname{new}}=y+\delta_{y}\in\mathbb{R}^{m}$, let
$s=c-A^{\top}y\in\mathbb{R}^{U}$ and
$S=P^{\top}\operatorname{diag}(s)P\in\mathbb{R}^{L\times L}$. Similarly define
$s^{\operatorname{new}}$ and $S^{\operatorname{new}}$ from
$y^{\operatorname{new}}$. Let
$H(y)=A\cdot\big{(}P\big{(}P^{\top}\operatorname{diag}(s)P\big{)}^{-1}P^{\top}\big{)}^{\circ
2}\cdot A^{\top}\in\mathbb{R}^{m\times m}$. If
$s,s^{\operatorname{new}}\in\Sigma_{n,2d}^{*}$, then $S$ and
$S^{\operatorname{new}}$ are both PSD, and we have

$\|S^{-1/2}S^{\operatorname{new}}S^{-1/2}-I\|_{F}=\|\delta_{y}\|_{H(y)}.$

###### Proof.

Note that if $s,s^{\operatorname{new}}\in\Sigma_{n,2d}^{*}$, then by the dual
cone characterization (Theorem 6.4) $S$ and $S^{\operatorname{new}}$ are both
PSD.

For convenience we define $M=PS^{-1}P^{\top}\in\mathbb{R}^{U\times U}$. Note
that $H(y)=A\cdot M^{\circ 2}\cdot A^{\top}$. We also define
$\delta_{s}=s^{\operatorname{new}}-s=-A^{\top}\delta_{y}$. $\forall u\in[U]$,
we use $p_{u}\in\mathbb{R}^{L}$ to denote the $u$-th row of $P$.

$\displaystyle\|S^{-1/2}S^{\operatorname{new}}S^{-1/2}-I\|_{F}^{2}$
$\displaystyle=\|S^{-1/2}\big{(}S^{\operatorname{new}}-S\big{)}S^{-1/2}\|_{F}^{2}$
(32)
$\displaystyle=\|S^{-1/2}\big{(}P\operatorname{diag}(\delta_{s})P^{\top}\big{)}S^{-1/2}\|_{F}^{2}$
$\displaystyle=\operatorname{tr}\big{(}S^{-1}(P^{\top}\operatorname{diag}(\delta_{s})P)S^{-1}(P^{\top}\operatorname{diag}(\delta_{s})P)\big{)}$
$\displaystyle=\operatorname{tr}\big{(}S^{-1}(\sum_{u\in
U}(\delta_{s})_{u}\cdot p_{u}p_{u}^{\top})S^{-1}(\sum_{v\in
U}(\delta_{s})_{v}\cdot p_{v}p_{v}^{\top})\big{)}$ $\displaystyle=\sum_{u,v\in
U}(\delta_{s})_{u}(\delta_{s})_{v}\cdot\operatorname{tr}\big{(}S^{-1}p_{u}p_{u}^{\top}S^{-1}p_{v}p_{v}^{\top}\big{)}$
$\displaystyle=\sum_{u,v\in
U}(\delta_{s})_{u}(\delta_{s})_{v}\cdot(p_{v}^{\top}S^{-1}p_{u})^{2}$
$\displaystyle=\sum_{u,v\in U}(\delta_{s})_{u}(\delta_{s})_{v}\cdot
M_{uv}^{2}\leavevmode\nobreak\ =\leavevmode\nobreak\ \|\delta_{s}\|_{M^{\circ
2}}^{2},$

where the third step follows from $\|A\|_{F}^{2}=\operatorname{tr}(A^{\top}A)$
and the cyclic property of trace, and the sixth step again follows from the
cyclic property of trace.

Since $\delta_{s}=-A^{\top}\delta_{y}$, we have

$\|\delta_{s}\|_{M^{\circ 2}}^{2}=\delta_{s}^{\top}M^{\circ
2}\delta_{s}=\delta_{y}^{\top}AM^{\circ
2}A^{\top}\delta_{y}=\delta_{y}^{\top}H(y)\delta_{y}=\|\delta_{y}\|_{H(y)}^{2}.$
(33)

Combining Eq. (32) and (33) we get the bound in the lemma statement. ∎

###  9.4 Proof of correctness

Finally we are ready to prove the correctness of Algorithm 1.

######  Theorem 9.3 (Correctness of Algorithm 1).

Consider the following optimization problem with $A\in\mathbb{R}^{m\times U}$,
$b\in\mathbb{R}^{m}$, and $c\in\mathbb{R}^{U}$:

$\begin{aligned} \text{Primal:\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ }\min\;&\langle c,x\rangle\quad\\\
\mathrm{s.t.}\leavevmode\nobreak\ \leavevmode\nobreak\ Ax&=b\\\
x&\in\Sigma_{n,2d}\,,\\\ \end{aligned}\quad\quad\quad\begin{aligned}
\text{Dual:\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
}\max\;&\langle y,b\rangle\\\ \mathrm{s.t.}\leavevmode\nobreak\
\leavevmode\nobreak\ A^{\top}y+s&=c\\\ s&\in\Sigma_{n,2d}^{*}\,.\\\
\end{aligned}$

Let $\mathrm{OPT}$ denote the optimal objective value of this optimization
problem. Assume Slater’s condition and that any primal feasible
$x\in\Sigma_{n,2d}$ satisfies $\|x\|_{1}\leq R$.

Then for any error parameters $\delta\in(0,1)$, $\epsilon_{S}\leq 0.01$, and
$\epsilon_{N}\leq 0.05$, Algorithm 1 outputs $x\in\Sigma_{n,2d}$ that
satisfies

$\langle c,x\rangle\leq\mathrm{OPT}+\delta\cdot
R\|c\|_{\infty}\quad\text{and}\quad\|Ax-b\|_{1}\leq 8\delta
L\cdot(LR\|A\|_{\infty}+\|b\|_{1}).$

###### Proof.

We consider the optimization problem $\min-b^{\top}y$ s.t.
$y\in\overline{D}_{F}$, where $F:\mathbb{R}^{m}\to\mathbb{R}_{+}$ is the
barrier function defined in Eq. (19), and $\overline{D}_{F}$ is the closure of
the domain of $F$. The barrier parameter of $F$ is $\nu_{F}=L$. This
optimization problem is equivalent to the dual formulation and its optimal
value is $-\mathrm{OPT}$. For any $\eta$, let $F_{\eta}(y)=-\eta
b^{\top}y+F(y)$.

In the beginning Algorithm 1 first uses Lemma 1.1 to convert the optimization
problem to another form which has an initial feasible solution $y$ that is
close to the optimal solution of $F_{\eta}$ with $\eta=1$. The initial $y$
satisfies $\|g_{\eta}(y)\|_{H(y)^{-1}}\leq\epsilon_{N}$ by Lemma 1.1.
Initially we also have $\widetilde{S}=S=P^{\top}\operatorname{diag}(s)P$ (Line
1 in Algorithm 1).

Next we prove the correctness of Algorithm 1 inductively. At each iteration,
we assume the following induction hypothesis is satisfied: (1)
$\|g_{\eta}(y)\|_{H(y)^{-1}}\leq\epsilon_{N}$, (2)
$\widetilde{S}\approx_{\epsilon_{S}}S$. We aim to prove that the updated
$y^{\operatorname{new}}$, $\eta^{\operatorname{new}}$,
$S^{\operatorname{new}}$, and $\widetilde{S}^{\operatorname{new}}$ still
satisfy these two conditions.

In Lemma 8.1 we have proved that in Algorithm 1 we always maintain
$N=\big{(}A\cdot(P\widetilde{S}^{-1}P^{\top})^{\circ 2}\cdot
A^{\top}\big{)}^{-1}$. Let $\widetilde{H}=N^{-1}$, we have

$\widetilde{H}=A\cdot\big{(}P\widetilde{S}^{-1}P^{\top}\big{)}^{\circ 2}\cdot
A^{\top}\approx_{2\epsilon_{S}}A\cdot\big{(}PS^{-1}P^{\top}\big{)}^{\circ
2}\cdot A^{\top}=H(y),$

where in the second step we use the induction hypothesis that
$\widetilde{S}\approx_{\epsilon_{S}}S$, and by Fact 5.4 we have
$\widetilde{S}^{-1}\approx_{\epsilon_{S}}S^{-1}$, and hence
$P\widetilde{S}^{-1}P^{\top}\approx_{\epsilon_{S}}PS^{-1}P^{\top}$, and hence
$(P\widetilde{S}^{-1}P^{\top})^{\circ
2}\approx_{2\epsilon_{S}}(PS^{-1}P^{\top})^{\circ 2}$.

The new vector $y^{\operatorname{new}}$ is computed as
$y^{\operatorname{new}}=y+\delta_{y}$ where
$\delta_{y}=-\widetilde{H}^{-1}g_{\eta}(y)$ (Line 1 and 1 of Algorithm 1). And
$\eta$ is updated to
$\eta^{\operatorname{new}}=\eta\cdot(1+\frac{\epsilon_{N}}{20\sqrt{L}})$ (Line
1 of Algorithm 1). Since $\|g_{\eta}(y)\|_{H(y)^{-1}}\leq\epsilon_{N}$, and
$\widetilde{H}\approx_{2\epsilon_{S}}H(y)$ where $2\epsilon_{S}\leq 0.02$ by
its definition in Algorithm 3, the requirements of Lemma 9.1 are satisfied, so
we have

$\|g_{\eta^{\operatorname{new}}}(y^{\operatorname{new}})\|_{H(y^{\operatorname{new}})^{-1}}\leq\epsilon_{N},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \text{and}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\|\delta_{y}\|_{H(y)}\leq 2\epsilon_{N}.$

This proves the first induction hypothesis.

Then using Lemma 9.2 and since $\epsilon_{N}\leq 0.01$ by its definition in
Algorithm 1, we have

$\|S^{-1/2}(S^{\operatorname{new}})S^{-1/2}-I\|_{F}\leq\|\delta_{y}\|_{H(y)}\leq
2\epsilon_{N}\leq 0.02.$

Thus the input matrix $S^{\operatorname{new}}$ to LowRankUpdate satisfies the
requirement of Eq. (31) of Lemma 9.1, and we have that
$\widetilde{S}^{\operatorname{new}}\approx_{\epsilon_{S}}S^{\operatorname{new}}$.
This proves the second induction hypothesis.

Finally, we know that after $t=40\epsilon_{N}^{-1}\sqrt{L}\log(L/\delta)$
iterations, $\eta$ becomes $(1+\frac{\epsilon_{N}}{20\sqrt{L}})^{t}\geq
2L/\delta^{2}$, so using Lemma 9.1, we have

$b^{\top}y\geq\mathrm{OPT}-\frac{\nu_{F}}{\eta}\cdot(1+2\epsilon_{N})\geq\mathrm{OPT}-\delta^{2}.$

Thus the initialization lemma (Lemma 1.1) ensures that we have a solution
$x\in\Sigma_{n,2d}$ to the original primal optimization problem which
satisfies

$\langle c,x\rangle\leq\mathrm{OPT}+\delta\cdot
R\|c\|_{\infty},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \|Ax-b\|_{1}\leq 8\delta
L\cdot(LR\|A\|_{\infty}+\|b\|_{1}).\qed$

##  10 Time complexity

###  10.1 Worst case time

We first bound the worst case running time of Algorithm 1. The running time of
the $i$-th iteration depends on the updated rank $r_{i}$ of LowRankUpdate,
which is defined to be the size of $V_{1}^{(i)}$ if
$V_{1}^{(i)}\neq\mathtt{null}$, and $U/L$ otherwise (see Lemma 9.1).

######  Lemma 10.1 (Worst case time of Algorithm 1).

In Algorithm 1, the initialization time is $O(U^{\omega})$, and the running
time in the $i$-th iteration is $O({\cal
T}_{\mathrm{mat}}(U,U,\min\\{Lr_{i},U\\}))$.

###### Proof.

Initialization time. The most time-consuming step of initialization is Line 1,
where computing $N=\big{(}A(PTP^{\top})^{\circ 2}A^{\top}\big{)}^{-1}$ takes
$O({\cal T}_{\mathrm{mat}}(U,U,L)+{\cal T}_{\mathrm{mat}}(U,U,m))$ time. This
is bounded by $O(U^{\omega})$ since $L,m\leq U$.

Time per iteration. In each iteration the most time-consuming steps are (1)
computing $S^{\operatorname{new}}$ and calling LowRankUpdate on Line 1-1, (2)
executing the if-clause on Line 1-1, and (3) computing
$g^{\operatorname{new}}$ on Line 1.

  1. 1.

Computing $S^{\operatorname{new}}$ on Line 1 takes ${\cal
T}_{\mathrm{mat}}(U,L,L)$ time. Calling LowRankUpdate on Line 1 takes
$O(L^{\omega})$ time by Lemma 9.1.

  2. 2.

In the if-clause on Line 1-1, if $V_{1}=V_{2}=\mathtt{null}$, then $Lr_{i}\geq
U$, and we compute
$N^{\operatorname{new}}=\big{(}A\cdot(PT^{\operatorname{new}}P^{\top})^{\circ
2}\cdot A^{\top}\big{)}^{-1}$, which takes $O({\cal T}_{\mathrm{mat}}(U,U,U))$
time. Otherwise we call UpdateHessianInv on Line 1, which takes $O({\cal
T}_{\mathrm{mat}}(U,U,Lr_{i}))$ time by Lemma 8.1. In total the if-clause has
running time $O({\cal T}_{\mathrm{mat}}(U,U,\min\\{Lr_{i},U\\}))$.

  3. 3.

Computing the gradient $g=-\eta^{\operatorname{new}}\cdot
b+A\cdot\operatorname{diag}\Big{(}P\big{(}P^{\top}\operatorname{diag}(s^{\operatorname{new}})P\big{)}^{-1}P^{\top}\Big{)}$
on Line 1 takes $O({\cal T}_{\mathrm{mat}}(U,U,L))$ time since $m\leq U$.

Thus the total time per iteration is $O({\cal
T}_{\mathrm{mat}}(U,U,\min\\{Lr_{i},U\\}))$. ∎

###  10.2 Amortized time

In this section we bound the amortized running time of Algorithm 1.

Let $\omega$ be the matrix multiplication exponent, let $\alpha$ be the dual
matrix multiplication exponent. The current best values are $\omega\approx
2.373$ and $\alpha\approx 0.314$ [l14, gu18, aw21]. Note that the current best
values of $\omega$ and $\alpha$ satisfies that $\alpha\geq 5-2\omega$. We use
the following modified lemma from [hjst21]:

######  Lemma 10.2 (Helpful lemma for amortization, modified version of Lemma
10.13 of [hjst21]).

Let $t$ denote the total number of iterations. Let $r_{i}\in[L]$ be the rank
for the $i$-th iteration for $i\in[t]$. Assume $r_{i}$ satisfies the following
condition: for any vector $g\in\mathbb{R}_{+}^{L}$ which is non-increasing, we
have $\sum_{i=1}^{t}r_{i}\cdot g_{r_{i}}\leq O(t\cdot\|g\|_{2})$.

If the cost in the $i$-th iteration is $O({\cal
T}_{\mathrm{mat}}(U,U,\min\\{Lr_{i},U\\}))$, when $\alpha\geq 5-2\omega$, the
amortized cost per iteration is $U^{2+o(1)}+U^{\omega-1/2+o(1)}\cdot L^{1/2}$.

For completeness we provide a proof of this lemma in Section 3.

######  Theorem 10.3 (Time of Algorithm 1).

When $\alpha\geq 5-2\omega$, the running time of Algorithm 1 is

$(U^{2}\cdot L^{1/2}+U^{\omega-1/2}\cdot L)\cdot(\log(1/\delta)+U^{o(1)}).$

###### Proof.

Using Lemma 10.1 the initialization time is $O(U^{\omega})\leq
O(U^{\omega-1/2}\cdot L)$ since $U\leq L^{2}$.

Using Lemma 9.1 we know that the ranks $r_{i}$ indeed satisfy the requirement
of Lemma 10.2, and since the worst case time per iteration is $O({\cal
T}_{\mathrm{mat}}(U,U,\min\\{Lr_{i},U\\}))$ (Lemma 10.1), using Lemma 10.2 the
time per iteration is $U^{2+o(1)}+U^{\omega-1/2+o(1)}\cdot L^{1/2}$. Since
there are in total $t=40\epsilon_{N}^{-1}\sqrt{L}\log(L/\delta)$ iterations,
we get the total running time claimed in the lemma statement. ∎

###  10.3 Comparison with previous results

In this section we compare the running time of [py19], [jiang2020faster,
hjst21], and our result. We assume that $m=\Theta(U)$ when making the
comparisons.

Ignoring $\log(1/\delta)$ and $U^{o(1)}$ factors, and since $L\leq U\leq
L^{2}$, the running times are

$\displaystyle\text{\cite[cite]{[\@@bibref{}{py19}{}{}]} (SOS)}:$
$\displaystyle\leavevmode\nobreak\ L^{0.5}\cdot U^{\omega},$
$\displaystyle\text{\cite[cite]{[\@@bibref{}{jiang2020faster, hjst21}{}{}]}
(SDP)}:$ $\displaystyle\leavevmode\nobreak\
L^{0.5}\cdot\min\\{UL^{2}+U^{\omega},\leavevmode\nobreak\
L^{4}+L^{2\omega-0.5}\\},$ $\displaystyle\text{Ours (SOS)}:$
$\displaystyle\leavevmode\nobreak\ L^{0.5}\cdot(U^{2}+U^{\omega-0.5}\cdot
L^{0.5}).$

#### Current $\omega$ and $\alpha$.

Plugging in the current best values $\omega\approx 2.373$ and $\alpha\approx
0.314$, we have

$\displaystyle\text{\cite[cite]{[\@@bibref{}{py19}{}{}]} (SOS)}:$
$\displaystyle\leavevmode\nobreak\ L^{0.5}\cdot U^{2.373},$
$\displaystyle\text{\cite[cite]{[\@@bibref{}{jiang2020faster, hjst21}{}{}]}
(SDP)}:$ $\displaystyle\leavevmode\nobreak\
L^{0.5}\cdot\min\\{UL^{2}+U^{2.373},\leavevmode\nobreak\ L^{4.246}\\}$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
L^{0.5}\cdot\begin{cases}UL^{2}&\text{when\leavevmode\nobreak\
}U\in(L,L^{1.457}],\\\ U^{2.373}&\text{when\leavevmode\nobreak\
}U\in(L^{1.457},L^{1.789}],\\\ L^{4.246}&\text{when\leavevmode\nobreak\
}U\in(L^{1.789},L^{2}),\end{cases}$ $\displaystyle\text{Ours (SOS)}:$
$\displaystyle\leavevmode\nobreak\ L^{0.5}\cdot(U^{2}+U^{1.873}L^{0.5}).$

Note that our running time is always better than the previous results, and for
several values of $L$ and $U$ we improve by a polynomial factor. See Figure 1
for an illustration.

##  11 Weighted SOS

In this section we provide the background of weighted sum-of-squares (WSOS)
optimization, following the notation of [py19].

Recall the motivation for sum-of-squares to solve polynomial optimization in
(10). Given $k$ non-zero polynomials $f_{1},\ldots,f_{k}$, we aim to optimize
over the polynomials that are non-negative on
$\mathcal{S}:=\\{x\in\mathbb{R}^{n}\mid f_{i}(x)\geq 0,\forall i\in[k]\\}$.
While polynomial optimization is hard in general, we can relax the problem to
optimize over the following set of non-negative polynomials over
$\mathcal{S}$.

######  Definition 11.1 (WSOS polynomials).

Let $\mathbf{f}=(f_{1},\ldots,f_{k})$ where $f_{i}$ are non-zero polynomials,
and let $\mathbf{d}=(d_{1},\ldots,d_{k})$ where $d_{i}>0$. Define
$\mathcal{V}_{n,2\mathbf{d}}^{\mathbf{f}}:=\\{\sum_{i=1}^{k}f_{i}s_{i}\mid
s_{i}\in\mathcal{V}_{n,2d_{i}},\forall i\in[k]\\}$. The WSOS polynomials are
defined as

$\Sigma_{n,2\mathbf{d}}^{\mathbf{f}}:=\Big{\\{}\sum_{i=1}^{k}f_{i}s_{i}\mid
s_{i}\in\Sigma_{n,2d_{i}},\forall
i\in[k]\Big{\\}}\subset\mathcal{V}_{n,2\mathbf{d}}^{\mathbf{f}}.$

We assume that the set $\mathcal{S}$ is unisolvent for
$\mathcal{V}_{n,2\mathbf{d}}^{\mathbf{f}}$, so that
$\Sigma_{n,2\mathbf{d}}^{\mathbf{f}}$ is a pointed and closed cone in
$\mathcal{V}_{n,2\mathbf{d}}^{\mathbf{f}}$. (See Proposition 6.1 of [py19].)
Note that w.l.o.g. one can add $f_{k+1}=1$ to $\mathbf{f}$ to capture the SOS
polynomials. Define

$U:=\dim(\mathcal{V}_{n,2\mathbf{d}}^{\mathbf{f}}),\leavevmode\nobreak\
\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\ \leavevmode\nobreak\
L_{i}:=\dim(\mathcal{V}_{n,d_{i}})={n+d_{i}\choose n},\forall i\in[k].$

Note that $U\geq L_{i}$ for all $i\in[k]$. We also define $L=\max_{i}L_{i}$.
Similar to SOS optimization, after fixing a basis $(q_{1},q_{2},\cdots,q_{U})$
of $\mathcal{V}_{n,2\mathbf{d}}^{\mathbf{f}}$, a polynomial in
$\mathcal{V}_{n,2\mathbf{d}}^{\mathbf{f}}$ corresponds to the vector of
coefficients in $\mathbb{R}^{U}$.

We consider the following optimization problem over the WSOS cone
($A\in\mathbb{R}^{m\times U}$, $c\in\mathbb{R}^{U}$, $b\in\mathbb{R}^{m}$):

$\begin{aligned} \text{Primal:\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ }\min\;&\langle c,x\rangle\quad\\\
\mathrm{s.t.}\leavevmode\nobreak\ \leavevmode\nobreak\ Ax&=b\\\
x&\in\Sigma_{n,2\mathbf{d}}^{\mathbf{f}}\,,\\\
\end{aligned}\quad\quad\quad\begin{aligned} \text{Dual:\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ }\max\;&\langle y,b\rangle\\\
\mathrm{s.t.}\leavevmode\nobreak\ \leavevmode\nobreak\ A^{\top}y+s&=c\\\
s&\in\Sigma_{n,2\mathbf{d}}^{\mathbf{f}\;*}\,.\\\ \end{aligned}$ (WSOS)

Here $\Sigma_{n,2\mathbf{d}}^{\mathbf{f}\;*}$ denotes the dual cone of
$\Sigma_{n,2\mathbf{d}}^{\mathbf{f}}$.

Similar to Theorem 6.4, the dual cone $\Sigma_{n,2\mathbf{d}}^{\mathbf{f}\;*}$
again admits a semidefinite characterisation. For any basis
$\mathbf{q}=(q_{1},q_{2},\cdots,q_{U})$ of
$\mathcal{V}_{n,2\mathbf{d}}^{\mathbf{f}}$, and any basis
$\mathbf{p}_{i}=(p_{i,1},\ldots,p_{i,L_{i}})$ of $\mathcal{V}_{n,d_{i}}$ for
$i\in[k]$, let $\Lambda_{i}:\mathbb{R}^{U}\to\mathbb{R}^{L_{i}\times L_{i}}$
be the unique linear mapping satisfying
$\Lambda_{i}(\mathbf{q})=\mathbf{p}_{i}\mathbf{p}_{i}^{\top}$. Then the dual
cone admits the characterization

$\Sigma_{n,2\mathbf{d}}^{\mathbf{f}\;*}=\\{s\in\mathbb{R}^{U}\mid\Lambda_{i}(s)\succeq
0,\forall i\in[k]\\}.$ (34)

For more details see Theorem 17.6 of [nesterov2000squared].

#### Barrier function under interpolant basis

[py19] extends the interpolant basis to WSOS. Let
$\mathcal{T}=\\{t_{1},\cdots,t_{U}\\}\subset\mathbb{R}^{n}$ be a unisolvent
set for $\mathcal{V}_{n,2\mathbf{d}}^{\mathbf{f}}$. Let
$\mathbf{q}=(q_{1},q_{2},\cdots,q_{U})$ be the Lagrange polynomials
corresponding to $\mathcal{T}$, and $\mathbf{q}$ forms a basis of
$\mathcal{V}_{n,2\mathbf{d}}^{\mathbf{f}}$. For all $i\in[k]$, choose any
basis $\mathbf{p}_{i}=(p_{i,1},\ldots,p_{i,L_{i}})$ of
$\mathcal{V}_{n,d_{i}}$. It was shown in [py19] that the unique linear mapping
$\Lambda_{i}:\mathbb{R}^{U}\to\mathbb{R}^{L_{i}\times L_{i}}$ satisfying
$\Lambda_{i}(\mathbf{q})=\mathbf{p}_{i}\mathbf{p}_{i}^{\top}$ is

$\Lambda_{i}(s):=P_{i}^{\top}\operatorname{diag}(\mathbf{f}_{i}\circ s)P_{i}.$
(35)

where $P_{i}\in\mathbb{R}^{U\times L_{i}}$ is defined as
$(P_{i})_{u,\ell}:=p_{i,\ell}(t_{u})$, and
$\mathbf{f}_{i}:=(f_{i}(t_{1}),\ldots,f_{i}(t_{U}))\in\mathbb{R}^{U}$.

Define $F_{i}:\mathbb{R}^{m}\to\mathbb{R}$ as
$F_{i}(y)=-\log\det(\Lambda_{i}(c-A^{\top}y))$. We get the following barrier
with barrier parameter at most $\sum_{i=1}^{k}L_{i}$:

$F^{\mathbf{f}}(y)=\sum_{i=1}^{k}F_{i}(y)=-\sum_{i=1}^{k}\log\det(\Lambda_{i}(c-A^{\top}y)).$

For any $\eta>0$ we define a function
$F_{\eta}^{\mathbf{f}}:\mathbb{R}^{m}\to\mathbb{R}$ as $F_{\eta}(y)=-\eta\cdot
b^{\top}y+F^{\mathbf{f}}(y)$.

The corresponding gradients and Hessians of the $F_{i}$’s are:

$\displaystyle g_{i}(y)$
$\displaystyle=A\cdot\Big{(}\mathbf{f}_{i}\circ\operatorname{diag}(P_{i}(P_{i}^{\top}\operatorname{diag}(\mathbf{f}_{i}\circ
s)P_{i})^{-1}P_{i}^{\top})\Big{)},$ $\displaystyle H_{i}(y)$
$\displaystyle=A\cdot\Big{(}(\mathbf{f}_{i}\mathbf{f}_{i}^{\top})\circ(P_{i}(P_{i}^{\top}\operatorname{diag}(\mathbf{f}_{i}\circ
s)P_{i})^{-1}P_{i}^{\top})^{\circ 2}\Big{)}\cdot A^{\top},$

and so gradient and Hessian of $F^{\mathbf{f}}_{\eta}(y)$ are

$g^{\mathbf{f}}_{\eta}(y)=-\eta b+\sum_{i=1}^{k}g_{i}(y),\leavevmode\nobreak\
\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\ \leavevmode\nobreak\
H^{\mathbf{f}}(y)=\sum_{i=1}^{k}H_{i}(y).$ (36)

Note that we omit the subscript $\eta$ for the Hessian since it is independent
of $\eta$.

###  11.1 Extension of our SOS algorithm to WSOS

We next show how to extend our SOS algorithm (Algorithm 1) to WSOS. The main
adjustments are the following.

  1. 1.

Maintain $k$ approximate slack. For each $i\in[k]$, we maintain a matrix
$\widetilde{S}_{i}\in\mathbb{R}^{L_{i}\times L_{i}}$ to approximate the exact
matrix
$S_{i}:=\Lambda_{i}(s)=P_{i}^{\top}\operatorname{diag}(\mathbf{f}_{i}\circ
s)P_{i}\in\mathbb{R}^{L_{i}\times L_{i}}$. We also maintain
$T_{i}=\widetilde{S}_{i}^{-1}$. In each iteration we compute
$S_{i}^{\operatorname{new}}=P^{\top}\operatorname{diag}(\mathbf{f}_{i}\circ
s^{\operatorname{new}})P_{i}\in\mathbb{R}^{L_{i}\times L_{i}}$

  2. 2.

LowRankUpdate on diagonal block matrices. In the $j$-th iteration, to compute
the low-rank update to the $\widetilde{S}_{i}$’s, we call LowRankUpdate
(Algorithm 3) with the two block diagonal matrices
$\operatorname{diag}(S_{1}^{\operatorname{new}},\cdots,S_{k}^{\operatorname{new}})$
and $\operatorname{diag}(\widetilde{S}_{1},\cdots,\widetilde{S}_{k})$ as
inputs. We do not need to explicitly construct these two block diagonal
matrices. It’s easy to check that in LowRankUpdate all the computations
maintain the block diagonal structure. Thus the outputs of LowRankUpdate are
also block diagonal matrices, and for all $i\in[k]$ we read off
$\widetilde{S}_{i}^{\operatorname{new}}\in\mathbb{R}^{L_{i}\times L_{i}}$, and
the low-rank updates $V_{i,1},V_{i,2}\in\mathbb{R}^{L_{i}\times r_{i,j}}$.

  3. 3.

Extend UpdateHessianInv to WSOS. We extend UpdateHessianInv (Algorithm 2), as
shown in UpdateHessianInvWSOS (Algorithm 4). Now we need to stack the $k$
updates of the $H_{i}$’s together, and perform a rank-$\sum_{i=1}^{k}r_{i,j}$
update. Using a similar proof as that of Lemma 8.1, we can show that

$\displaystyle
T_{i}^{\operatorname{new}}=(\widetilde{S}_{i}^{\operatorname{new}})^{-1}\in\mathbb{R}^{L_{i}\times
L_{i}},\leavevmode\nobreak\ \leavevmode\nobreak\ \forall i\in[k],$
$\displaystyle
N^{\operatorname{new}}=\Big{(}A\cdot\big{(}\sum_{i=1}^{k}(\mathbf{f}_{i}\mathbf{f}_{i}^{\top})\circ(P_{i}(\widetilde{S}_{i}^{\operatorname{new}})^{-1}P_{i}^{\top})^{\circ
2}\big{)}\cdot A^{\top}\Big{)}^{-1}\in\mathbb{R}^{m\times m}.$

Input : $T_{i}\in\mathbb{R}^{L_{i}\times L_{i}}$ for all $i\in[k]$,
$N\in\mathbb{R}^{m\times m}$, $V_{i,1},V_{i,2}\in\mathbb{R}^{L_{i}\times
r_{i}}$ for all $i\in[k]$

Output : $T_{i}^{\operatorname{new}}\in\mathbb{R}^{L_{i}\times L_{i}}$ for all
$i\in[k]$, $N^{\operatorname{new}}\in\mathbb{R}^{m\times m}$

1

2// Step 1

3

4for  _$i=1,2,\cdots,k$_ do

5 $\overline{V_{i,1}}\leftarrow-TV_{i,1}\cdot(I+V_{i,2}^{\top}TV_{i,1})^{-1}$
;

6 // $\overline{V_{i,1}}\in\mathbb{R}^{L_{i}\times r_{i}}$

7

8 $\overline{V_{i,2}}\leftarrow TV_{i,2}$ ;

9 // $\overline{V_{i,2}}\in\mathbb{R}^{L_{i}\times r_{i}}$

10

11 $T_{i}^{\operatorname{new}}\leftarrow
T_{i}+\overline{V_{i,1}}\cdot\overline{V_{i,2}}^{\top}$ ;

12 // $T_{i}^{\operatorname{new}}\in\mathbb{R}^{L_{i}\times L_{i}}$

13

14

15// Step 2

16

17for  _$i=1,2,\cdots,k$_ do

18
$Y^{\prime}_{i}\leftarrow\operatorname{diag}(\mathbf{f}_{i})\cdot[2P_{i}T_{i},P_{i}\overline{V_{i,1}}]$
;

19 // $Y^{\prime}_{i}\in\mathbb{R}^{U\times(L_{i}+r_{i})}$

20

21
$Z^{\prime}_{i}\leftarrow\operatorname{diag}(\mathbf{f}_{i})\cdot[P_{i},P_{i}\overline{V_{i,2}}]$
;

22 // $Z^{\prime}_{i}\in\mathbb{R}^{U\times(L_{i}+r_{i})}$

23

24
$Y_{i}\leftarrow[\operatorname{diag}(u_{i,1})Y^{\prime}_{i},\cdots,\operatorname{diag}(u_{i,r})Y^{\prime}_{i}]$,
$u_{i,j}$ is $j$-th col. of $P_{i}\overline{V_{i,1}}$ ;

25 // $Y_{i}\in\mathbb{R}^{U\times(L_{i}+r_{i})r_{i}}$

26

27
$Z_{i}\leftarrow[\operatorname{diag}(v_{i,1})Z^{\prime}_{i},\cdots,\operatorname{diag}(v_{i,r})Z^{\prime}_{i}]$,
$v_{i,j}$ is $j$-th col. of $P_{i}\overline{V_{i,2}}$ ;

28 // $Z_{i}\in\mathbb{R}^{U\times(L_{i}+r_{i})r_{i}}$

29

30

31$Y\leftarrow[Y_{1},\cdots,Y_{k}]$, $Z\leftarrow[Z_{1},\cdots,Z_{k}]$ ;

32 // $Y,Z\in\mathbb{R}^{U\times\sum_{i=1}^{k}(L_{i}+r_{i})r_{i}}$

33

34// Step 3

35

36$N^{\operatorname{new}}\leftarrow
N-N\cdot(AY)\cdot\big{(}I+(AZ)^{\top}N(AY)\big{)}^{-1}\cdot(AZ)^{\top}\cdot N$
;

37 // $N^{\operatorname{new}}\in\mathbb{R}^{m\times m}$

38

return _$T^{\operatorname{new}},N^{\operatorname{new}}$_

Algorithm 4 UpdateHessianInvWSOS

###  11.2 Correctness

In this section we prove the correctness of the WSOS algorithm described in
the previous section. In the previous section we have shown that an analogue
of Lemma 8.1 holds for WSOS. Lemma 9.1 and Lemma 9.1 directly hold for our
WSOS algorithm since they can be applied to any barrier function. It remains
to prove the following analogue of Lemma 9.2.

######  Lemma 11.2 (Slowly moving guarantee for WSOS).

Let $s,s^{\operatorname{new}}\in\mathbb{R}^{U}$ be the current and the updated
slack variables of the WSOS algorithm. For each $i\in[k]$, define
$S_{i}=\Lambda_{i}(s)=P_{i}^{\top}\operatorname{diag}(\mathbf{f}_{i}\circ
s)P_{i}$, and
$S_{i}^{\operatorname{new}}=\Lambda_{i}(s^{\operatorname{new}})=P_{i}^{\top}\operatorname{diag}(\mathbf{f}_{i}\circ
s^{\operatorname{new}})P_{i}$. Let
$\overline{S}:=\operatorname{diag}(S_{1},\ldots,S_{k})$, and
$\overline{S}^{\operatorname{new}}:=\operatorname{diag}(S_{1}^{\operatorname{new}},\ldots,S_{k}^{\operatorname{new}})$.
Then we have

$\|\overline{S}^{-1/2}\overline{S}^{\operatorname{new}}\overline{S}^{-1/2}-I\|_{F}^{2}=\|\delta_{y}\|_{H^{\mathbf{f}}(y)}^{2},$

where $H^{\mathbf{f}}(y)$ is the Hessian of the barrier function
$F^{\mathbf{f}}_{\eta}(y)$, as defined in Eq. (36).

###### Proof.

For each $i\in[k]$, we define $M_{i}:=P_{i}S^{-1}P_{i}^{\top}$, and let
$p_{i,u}$ denote the $u$-th row of $P_{i}$ for each $u\in[U]$. Let
$\delta_{s}=-A^{\top}\delta y$, and note that
$s^{\operatorname{new}}-s=\delta_{s}$. We have that

$\displaystyle\|S_{i}^{-1/2}S_{i}^{\operatorname{new}}S_{i}^{-1/2}-I\|_{F}^{2}$
$\displaystyle=\|S_{i}^{-1/2}\big{(}S_{i}^{\operatorname{new}}-S_{i}\big{)}S_{i}^{-1/2}\|_{F}^{2}$
$\displaystyle=\|S_{i}^{-1/2}\big{(}P_{i}^{\top}\operatorname{diag}(\mathbf{f}_{i}\circ\delta_{s})P_{i}\big{)}S_{i}^{-1/2}\|_{F}^{2}$
$\displaystyle=\operatorname{tr}\Big{(}S_{i}^{-1}\big{(}P_{i}^{\top}\operatorname{diag}(\mathbf{f}_{i}\circ\delta_{s})P\big{)}S_{i}^{-1}\big{(}P_{i}^{\top}\operatorname{diag}(\mathbf{f}_{i}\circ\delta_{s})P_{i}\big{)}\Big{)}$
$\displaystyle=\operatorname{tr}\Big{(}S_{i}^{-1}\big{(}\sum_{u\in
U}(\mathbf{f}_{i})_{u}(\delta_{s})_{u}\cdot
p_{i,u}p_{i,u}^{\top}\big{)}S_{i}^{-1}\big{(}\sum_{v\in
U}(\mathbf{f}_{i})_{v}(\delta_{s})_{v}\cdot
p_{i,v}p_{i,v}^{\top}\big{)}\Big{)}$ $\displaystyle=\sum_{u,v\in
U}(\mathbf{f}_{i})_{u}(\mathbf{f}_{i})_{v}(\delta_{s})_{u}(\delta_{s})_{v}\cdot\operatorname{tr}\big{(}S_{i}^{-1}p_{i,u}p_{i,u}^{\top}S_{i}^{-1}p_{i,v}p_{i,v}^{\top}\big{)}$
$\displaystyle=\sum_{u,v\in
U}(\mathbf{f}_{i})_{u}(\mathbf{f}_{i})_{v}(\delta_{s})_{u}(\delta_{s})_{v}\cdot(p_{i,v}^{\top}S_{i}^{-1}p_{i,u})^{2}$
$\displaystyle=\sum_{u,v\in
U}(\mathbf{f}_{i})_{u}(\mathbf{f}_{i})_{v}(\delta_{s})_{u}(\delta_{s})_{v}\cdot(M_{i})_{uv}^{2}\leavevmode\nobreak\
=\leavevmode\nobreak\
\|\delta_{s}\|_{(\mathbf{f}_{i}\mathbf{f}_{i}^{\top})\circ{M_{i}^{\circ
2}}}^{2},$

where the third step follows from $\|A\|_{F}^{2}=\operatorname{tr}(A^{\top}A)$
and the cyclic property of trace, and the sixth step again follows from the
cyclic property of trace.

Since $\delta_{s}=-A^{\top}\delta y$ and
$H_{i}(y)=A\big{(}(\mathbf{f}_{i}\mathbf{f}_{i}^{\top})\circ{M_{i}^{\circ
2}}\big{)}A^{\top}$, we have

$\|\delta_{s}\|_{(\mathbf{f}_{i}\mathbf{f}_{i}^{\top})\circ{M_{i}^{\circ
2}}}=\delta_{s}^{\top}\big{(}(\mathbf{f}_{i}\mathbf{f}_{i}^{\top})\circ{M_{i}^{\circ
2}}\big{)}\delta_{s}=\delta_{y}^{\top}A\big{(}(\mathbf{f}_{i}\mathbf{f}_{i}^{\top})\circ{M_{i}^{\circ
2}}\big{)}A^{\top}\delta_{y}=\|\delta_{y}\|^{2}_{H_{i}(y)},$

and so

$\|S_{i}^{-1/2}S_{i}^{\operatorname{new}}S_{i}^{-1/2}-I\|_{F}^{2}=\|\delta_{y}\|_{H_{i}(y)}^{2}.$

Since $\overline{S}:=\operatorname{diag}(S_{1},\ldots,S_{k})$, and
$\overline{S}^{\operatorname{new}}:=\operatorname{diag}(S_{1}^{\operatorname{new}},\ldots,S_{k}^{\operatorname{new}})$,
we have

$\displaystyle\|\overline{S}^{-1/2}\overline{S}^{\operatorname{new}}\overline{S}^{-1/2}-I\|_{F}=$
$\displaystyle\leavevmode\nobreak\
\Big{(}\sum_{i=1}^{k}\|S_{i}^{-1/2}S_{i}^{\operatorname{new}}S_{i}^{-1/2}-I\|_{F}^{2}\Big{)}^{1/2}$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\sum_{i=1}^{k}\|\delta_{y}\|_{H_{i}(y)}^{2}=\|\delta_{y}\|_{H^{\mathbf{f}}(y)}^{2},$

where the last step follows from $H^{\mathbf{f}}(y)=\sum_{i=1}^{k}H_{i}(y)$. ∎

Using Section 9.1 with the WSOS barrier $F^{\mathbf{f}}$, we can bound
$\|\delta_{y}\|_{H^{\mathbf{f}}(y)}\leq 2\epsilon_{N}$ throughout the
algorithm. Combining with Lemma 11.2, the requirement of Lemma 9.1 is
satisfied for the WSOS algorithm where the inputs to LowRankUpdate are block
diagonal matrices
$\operatorname{diag}(S_{1}^{\operatorname{new}},\cdots,S_{k}^{\operatorname{new}})$
and $\operatorname{diag}(\widetilde{S}_{1},\cdots,\widetilde{S}_{k})$. In the
$j$-th iteration let $r_{i,j}$ denote the rank of the update to
$\widetilde{S}_{i}$, and let $r_{j}=\sum_{i=1}^{k}r_{i,j}$. Then for any
vector $g\in\mathbb{R}_{+}^{n}$ which is non-increasing, we have

$\sum_{j=1}^{t}r_{j}g_{r_{j}}\leq O(T\cdot\|g\|_{2}\cdot\log L).$ (37)

Now we have proved all the analogous lemmas for the WSOS algorithm. The
correctness of the WSOS algorithm follows from a similar argument at that of
Theorem 9.3. We summarize this in the following theorem.

######  Theorem 11.3 (Correctness of WSOS Algorithm).

Consider the following optimization problem with $A\in\mathbb{R}^{m\times U}$,
$b\in\mathbb{R}^{m}$, and $c\in\mathbb{R}^{U}$:

$\begin{aligned} \text{Primal:\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ }\min\;&\langle c,x\rangle\quad\\\
\mathrm{s.t.}\leavevmode\nobreak\ \leavevmode\nobreak\ Ax&=b\\\
x&\in\Sigma_{n,2\mathbf{d}}^{\mathbf{f}}\,,\\\
\end{aligned}\quad\quad\quad\begin{aligned} \text{Dual:\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ }\max\;&\langle y,b\rangle\\\
\mathrm{s.t.}\leavevmode\nobreak\ \leavevmode\nobreak\ A^{\top}y+s&=c\\\
s&\in\Sigma_{n,2\mathbf{d}}^{\mathbf{f}\;*}\,.\\\ \end{aligned}$

Let $\mathrm{OPT}$ denote the optimal objective value of this optimization
problem. Assume that any primal feasible
$x\in\Sigma_{n,2\mathbf{d}}^{\mathbf{f}}$ satisfies $\|x\|_{1}\leq R$.

Then for any error parameters $\delta\in(0,1)$ and
$\epsilon_{S},\epsilon_{N}\leq 0.01$, the variant of Algorithm 1 described in
Section 11.1 outputs $x\in\Sigma_{n,2\mathbf{d}}^{\mathbf{f}}$ that satisfies

$\langle c,x\rangle\leq\mathrm{OPT}+\delta\cdot
R\|c\|_{\infty}\quad\text{and}\quad\|Ax-b\|_{1}\leq 8\delta
L\cdot(LR\|A\|_{\infty}+\|b\|_{1}).$

###  11.3 Time analysis

In this section we analyse the running time of the WSOS algorithm described in
Section 11.1. Recall that we define $L=\max_{i}L_{i}$.

#### Worst-case time

The worst-case time complexity of the WSOS algorithm can be analyzed in a
similar way as that of Lemma 10.1. We only mention the bottlenecks.

  * •

In initialization computing the Hessian inverse takes $O(k{\cal
T}_{\mathrm{mat}}(U,U,L)+U^{\omega})$ time.

  * •

Similar to the SOS algorithm, in each iteration, there are two main
bottlenecks:

    1. 1.

The call to UpdateHessianInvWSOS, which takes ${\cal
T}_{\mathrm{mat}}(U,U,\min(U,Lr_{j}))$ time in iteration $j\in[t]$

    2. 2.

Computing all $P_{i}(P_{i}^{\top}\operatorname{diag}(\mathbf{f}_{i}\circ
s^{\operatorname{new}})P_{i})^{-1}P_{i}^{\top}$ when computing the gradient,
which takes $k{\cal T}_{\mathrm{mat}}(U,U,L)$ time.

Thus the total time per iteration is $O({\cal
T}_{\mathrm{mat}}(U,U,\min(Lr_{j},U))+k{\cal T}_{\mathrm{mat}}(U,U,L))$.

#### Amortized time

We can bound the amortized time of our WSOS algorithm using a similar argument
as that of the SOS algorithm in Section 10.2. Again assume that $\alpha\geq
5-2\omega$. Using Lemma 10.2 and Equation 37 that we proved in the previous
section, the amortized iteration complexity is

$O\big{(}U^{2}+U^{\omega-1/2}L^{1/2}+k{\cal
T}_{\mathrm{mat}}(U,U,L)\big{)}\cdot(\log(1/\delta)+U^{o(1)}).$

Overall we have $\sqrt{kL}$ iterations (since the barrier parameter is at most
$\sum_{i=1}^{k}L_{i}\leq kL$), resulting in total complexity of

$O\big{(}(kL)^{0.5}\cdot(U^{2}+U^{\omega-1/2}L^{1/2}+k{\cal
T}_{\mathrm{mat}}(U,U,L))+U^{\omega}\big{)}\cdot(\log(1/\delta)+U^{o(1)}).$

#### Comparison with previous results

Again assume that the number of constraints $m=\Theta(U)$. The running times
of [py19], [jiang2020faster], [hjst21] and our algorithm are

$\displaystyle\text{\cite[cite]{[\@@bibref{}{py19}{}{}]} (WSOS)}:$
$\displaystyle\leavevmode\nobreak\
\tilde{O}\big{(}(kL)^{0.5}\cdot(U^{\omega}+k{\cal
T}_{\mathrm{mat}}(U,U,L))\big{)},$
$\displaystyle\text{\cite[cite]{[\@@bibref{}{jiang2020faster, hjst21}{}{}]}
(SDP)}:$ $\displaystyle\leavevmode\nobreak\
\tilde{O}\big{(}(kL)^{0.5}\cdot\min\\{kUL^{2}+U^{\omega}+(kL)^{\omega},\leavevmode\nobreak\
U^{2}+(kL)^{4}+U^{\omega}/(kL)^{0.5}+(kL)^{2\omega-0.5}\\}\big{)},$
$\displaystyle\text{Ours (WSOS)}:$ $\displaystyle\leavevmode\nobreak\
\tilde{O}\big{(}(kL)^{0.5}\cdot(U^{2}+U^{\omega-1/2}L^{1/2}+k{\cal
T}_{\mathrm{mat}}(U,U,L))+U^{\omega}).$

#### Reformulation as SDP

We briefly explain the running time of WSOS when reformulated as SDP. Let
$\mathbf{q}=(q_{1},q_{2},\cdots,q_{U})$ be a basis of
$\mathcal{V}_{n,2\mathbf{d}}^{\mathbf{f}}$. For each $i\in[k]$, let
$\mathbf{p}_{i}=(p_{i,1},\ldots,p_{i,L_{i}})$ be a basis of
$\mathcal{V}_{n,d_{i}}$. Let
$\Lambda_{i}:\mathbb{R}^{U}\to\mathbb{R}^{L_{i}\times L_{i}}$ be the unique
linear mapping satisfying
$\Lambda_{i}(\mathbf{q})=\mathbf{p}_{i}\mathbf{p}_{i}^{\top}$. Let
$\Lambda_{i}^{*}:\mathbb{R}^{L_{i}\times L_{i}}\to\mathbb{R}^{U}$ be the
adjoint of $\Lambda$, i.e., $\Lambda_{i}^{*}$ is the unique linear operator
that satisfies $\langle\Lambda_{i}(x),V\rangle=\langle
x,\Lambda_{i}^{*}(V)\rangle$ for all $x\in\mathbb{R}^{U}$ and
$V\in\mathbb{R}^{L_{i}\times L_{i}}$. Apart from the dual characterization
described in Eq. (34), the WSOS cone also has a primal characterization (also
proved in Theorem 17.6 of [nesterov2000squared]):

$\Sigma_{n,2\mathbf{d}}^{\mathbf{f}}=\big{\\{}x\in\mathbb{R}^{U}\mid
x=\sum_{i=1}^{k}\Lambda_{i}^{*}(V_{i}),\leavevmode\nobreak\
\leavevmode\nobreak\ V_{i}\succeq 0\in\mathbb{R}^{L_{i}\times L_{i}}\text{ for
}i\in[k]\big{\\}}.$

We can therefore reformulate the primal of (WSOS) as the following SDP (we use
$A_{j,*}\in\mathbb{R}^{U}$ to denote the $j$-th row of
$A\in\mathbb{R}^{m\times U}$):

$\displaystyle\min$
$\displaystyle\sum_{i=1}^{k}\langle\Lambda_{i}(c),V_{i}\rangle\quad$ (WSOS -
SDP) $\displaystyle\mathrm{s.t.}\leavevmode\nobreak\ \leavevmode\nobreak$
$\displaystyle\sum_{i=1}^{k}\langle\Lambda_{i}(A_{j,*}),V_{i}\rangle=b_{j},\quad\forall
j\in[m]$ $\displaystyle V_{i}\succeq 0,\quad\forall i\in[k].$

Note that the number of constraints is still $m$, unchanged from
$\eqref{WSOS_program_primal_dual}$, but the number of variables is now
$\sum_{i=1}^{k}L_{i}^{2}=O(kL^{2})$, so just processing the constraint matrix
in every iteration of an IPM algorithm costs up to $O(kUL^{2})$. The SDP
solvers [jiang2020faster] and [hjst21] can profit from the block-diagonal
structure of the PSD matrices, which leads to the runtimes claimed in the
previous paragraph.

##  12 Discussion of Bottleneck

The dominating term in the cost-per-iteration of our SOS algorithm is ${\cal
T}_{\mathrm{mat}}(U,U,Lr)$. A natural question is whether this can be improved
to ${\cal T}_{\mathrm{mat}}(U,U,L)$.

The bottleneck for this potential improvement is the following self-contained
batch matrix-product problem, which is interesting in its own right:

  * •

Input : invertible matrix $H\in\mathbb{R}^{U\times U}$, and
$Z=[\operatorname{diag}(u_{1})Y,\cdots,\operatorname{diag}(u_{r})Y]\in\mathbb{R}^{U\times
Lr}$ for some matrix $Y\in\mathbb{R}^{U\times L}$ and $r$ vectors
$u_{1},\cdots,u_{r}\in\mathbb{R}^{U}$.

  * •

Output : The inverse $(H+ZZ^{\top})^{-1}-H^{-1}$ under the assumption that
$(H+ZZ^{\top})^{-1}$ exists.

Naively computing this inverse using Woodbury’s identity via

$(H+ZZ^{\top})^{-1}-H^{-1}=H^{-1}Z(I+Z^{\top}Z)^{-1}Z^{\top}H^{-1}$ (38)

takes ${\cal T}_{\mathrm{mat}}(U,U,Lr)$ time. However, since the update matrix
$Z$ only depends on $Ur+UL$ variables ($Y$ and $u_{1},\cdots,u_{r}$), there is
no known lower bound that prohibits a better running time of ${\cal
T}_{\mathrm{mat}}(U,U,L)+{\cal T}_{\mathrm{mat}}(U,U,r)$.

Computing $H^{-1}Z$ in (38) reduces to the following fundamental problem:
Given matrices $A,C\in\mathbb{R}^{n\times n}$ and vectors $b_{1},\ldots
b_{n}\in\mathbb{R}^{n}$ the task is to compute
$D_{k}=A\operatorname{diag}(b_{k})C$ for all $k\in[n]$. Let
$d_{ijk}:=(D_{k})_{i,j}$. Then using $b_{k\ell}:=(b_{k})_{\ell}$ we can write

$d_{ijk}=\sum_{\ell=1}^{n}a_{i\ell}b_{k\ell}c_{\ell j}.$

For symmetry let us harmlessly permute the entries in $c$ to obtain

$d_{ijk}=\sum_{\ell=1}^{n}a_{i\ell}b_{k\ell}c_{j\ell}.$

Note that this is exactly the term for matrix multiplication extended to three
matrices. For a matrix $F=AB^{\top}$ we would have
$f_{ij}=\sum_{\ell}a_{i\ell}b_{j\ell}$. From the paragraphs above it is clear
that we can compute all $d_{ijk}$ in time ${\cal
T}_{\mathrm{mat}}(n,n,n^{2})$. But can this be done faster? Interestingly, the
problem of finding a 4-clique in a graph can be reduced to this problem,
further proving its significance. Note that the currently best known running
time for this problem matches ${\cal T}_{\mathrm{mat}}(n,n,n^{2})$ for
$n$-vertex graphs [eisenbrand2004complexity].

#### Acknowledgments

The second author would like to thank Vissarion Fisikopoulos and Elias
Tsigaridas for introducing him from a practical perspective to Sum-of-Square
optimization under the interpolant basis.

## References

  * [AW21] Josh Alman and Virginia Vassilevska Williams.  A refined laser method and faster matrix multiplication.  In Proceedings of the 2021 ACM-SIAM Symposium on Discrete Algorithms (SODA), pages 522–539. SIAM, 2021.
  * [BBC+09] Brandon Ballinger, Grigoriy Blekherman, Henry Cohn, Noah Giansiracusa, Elizabeth Kelly, and Achill Schürmann.  Experimental study of energy-minimizing point configurations on spheres.  Experimental Mathematics, 18(3):257–283, 2009.
  * [BDMSV10] L. Bos, S. De Marchi, A. Sommariva, and M. Vianello.  Computing multivariate fekete and leja points by numerical linear algebra.  SIAM Journal on Numerical Analysis, 48(5):1984–1999, 2010.
  * [BHK+19] Boaz Barak, Samuel B. Hopkins, Jonathan A. Kelner, Pravesh K. Kothari, Ankur Moitra, and Aaron Potechin.  A nearly tight sum-of-squares lower bound for the planted clique problem.  SIAM J. Comput., 48(2):687–735, 2019.
  * [Blä13] Markus Bläser.  Fast matrix multiplication.  Theory of Computing, pages 1–60, 2013.
  * [BPT12] Grigoriy Blekherman, Pablo A Parrilo, and Rekha R Thomas.  Semidefinite optimization and convex algebraic geometry.  SIAM, 2012.
  * [BRS11] Boaz Barak, Prasad Raghavendra, and David Steurer.  Rounding semidefinite programming hierarchies via global correlation, 2011\.
  * [BV06] Christine Bachoc and Frank Vallentin.  New upper bounds for kissing numbers from semidefinite programming.  Technical report, Journal of the American Mathematical Society, 2006.
  * [CLS19] Michael B Cohen, Yin Tat Lee, and Zhao Song.  Solving linear programs in the current matrix multiplication time.  In Proceedings of the 51st Annual ACM Symposium on Theory of Computing (STOC), 2019.
  * [EG04] Friedrich Eisenbrand and Fabrizio Grandoni.  On the complexity of fixed parameter clique and dominating set.  Theoretical Computer Science, 326(1-3):57–67, 2004.
  * [GMM16] Bissan Ghaddar, Jakub Marecek, and M. Mevissen.  Optimal power flow as a polynomial optimization problem.  IEEE Transactions on Power Systems, 31:539–546, 2016.
  * [GU18] François Le Gall and Florent Urrutia.  Improved rectangular matrix multiplication using powers of the coppersmith-winograd tensor.  In Proceedings of the 2018 ACM-SIAM Symposium on Discrete Algorithms (SODA), pages 1029–1046. SIAM, 2018.
  * [HHLP16] Roxana Heß, Didier Henrion, Jean-Bernard Lasserre, and Tien Son Pham.  Semidefinite approximations of the polynomial abscissa.  SIAM J. Control. Optim., 54(3):1633–1656, 2016.
  * [HJST21] Baihe Huang, Shunhua Jiang, Zhao Song, and Runzhou Tao.  Solving tall dense sdps in the current matrix multiplication time, 2021\.
  * [HKP+17] Samuel B. Hopkins, Pravesh K. Kothari, Aaron Potechin, Prasad Raghavendra, Tselil Schramm, and David Steurer.  The power of sum-of-squares for detecting hidden structures.  In Chris Umans, editor, 58th IEEE Annual Symposium on Foundations of Computer Science, FOCS 2017, Berkeley, CA, USA, October 15-17, 2017, pages 720–731. IEEE Computer Society, 2017.
  * [HL18] Samuel B. Hopkins and Jerry Li.  Mixture models, robustness, and sum of squares proofs.  In Ilias Diakonikolas, David Kempe, and Monika Henzinger, editors, Proceedings of the 50th Annual ACM SIGACT Symposium on Theory of Computing, STOC 2018, Los Angeles, CA, USA, June 25-29, 2018, pages 1021–1034. ACM, 2018.
  * [JKL+20] Haotian Jiang, Tarun Kathuria, Yin Tat Lee, Swati Padmanabhan, and Zhao Song.  A faster interior point method for semidefinite programming, 2020.
  * [Kar84] Narendra Karmarkar.  A new polynomial-time algorithm for linear programming.  In Proceedings of the sixteenth annual ACM symposium on Theory of computing, pages 302–311, 1984.
  * [Las15] Jean Bernard Lasserre.  An Introduction to Polynomial and Semi-Algebraic Optimization.  Cambridge Texts in Applied Mathematics. Cambridge University Press, 2015\.
  * [Lau09] M. Laurent.  Sums of squares, moment matrices and optimization over polynomials, pages 155–270.  Number 149 in The IMA Volumes in Mathematics and its Applications Series. Springer Verlag, Germany, 2009.
  * [LG14] François Le Gall.  Powers of tensors and fast matrix multiplication.  In Proceedings of the 39th international symposium on symbolic and algebraic computation, pages 296–303, 2014.
  * [LS19] Yin Tat Lee and Aaron Sidford.  Solving linear programs with sqrt(rank) linear system solves.  CoRR, abs/1910.08033, 2019.
  * [LSZ19] Yin Tat Lee, Zhao Song, and Qiuyi Zhang.  Solving empirical risk minimization in the current matrix multiplication time.  CoRR, abs/1905.04447, 2019.
  * [Nes00] Yuri Nesterov.  Squared functional systems and optimization problems.  In High performance optimization, pages 405–440. Springer, 2000\.
  * [NN94] Y. Nesterov and A. Nemirovskii.  Interior-point polynomial algorithms in convex programming.  Studies in Applied Mathematics. Society for Industrial and Applied Mathematics, 1994.
  * [Pan01] Victor Y. Pan.  Structured Matrices and Polynomials.  Birkhäuser Boston, 2001.
  * [Pap12] Dávid Papp.  Optimal designs for rational function regression.  Journal of the American Statistical Association, 107(497):400–411, 2012.
  * [Par20] Pablo Parrilo.  Sum of squares : theory and applications : AMS short course, sum of squares : theory and applications, January 14-15, 2019, Baltimore, Maryland.  American Mathematical Society, Providence, Rhode Island, 2020.
  * [PV99] Mihai Putinar and Florian-Horia Vasilescu.  Positive polynomials on semi-algebraic sets.  Comptes Rendus de l’Académie des Sciences - Series I - Mathematics, 328(7):585–589, 1999.
  * [PY19] Dávid Papp and Sercan Yildiz.  Sum-of-squares optimization without semidefinite programming.  SIAM Journal on Optimization, 29(1):822–851, 2019.
  * [RDV07] Tae Roh, Bogdan Dumitrescu, and Lieven Vandenberghe.  Multidimensional FIR filter design via trigonometric sum-of-squares optimization.  J. Sel. Topics Signal Processing, 1(4):641–650, 2007.
  * [Ren01] James Renegar.  A Mathematical View of Interior-Point Methods in Convex Optimization.  Society for Industrial and Applied Mathematics, January 2001.
  * [Str87] Gilbert Strang.  Karmarkar’s algorithm and its place in applied mathematics.  The Mathematical Intelligencer, 9(2):4–10, 1987.
  * [SV09] Alvise Sommariva and Marco Vianello.  Computing approximate fekete points by qr factorizations of vandermonde matrices.  Computers & Mathematics with Applications, 57(8):1324–1336, 2009\.
  * [Tan15] Ning Tan.  On the Power of Lasserre SDP Hierarchy.  PhD thesis, EECS Department, University of California, Berkeley, Dec 2015\.
  * [Vai89] Pravin M Vaidya.  Speeding-up linear programming using fast matrix multiplication.  In 30th Annual Symposium on Foundations of Computer Science (FOCS), pages 332–337. IEEE, 1989.
  * [YTM94] Yinyu Ye, Michael J Todd, and Shinji Mizuno.  An $o\sqrt{nL}$-iteration homogeneous and self-dual linear programming algorithm.  Mathematics of operations research, 19(1):53–67, 1994.