$W_{\mathrm{iso}}$ be the set of isolated vertices (vertices not incident to
any hyperedges). We argue that
$|V(\kappa)|-|S_{\min}|+|W_{\mathrm{iso}}|<2q+1+\lceil\sum_{i\in[2q]}a_{i}(\kappa)/2\rceil$.

To understand this, we can start from the extreme case with a shape
$\kappa_{0}$ where all vertices in
$\widehat{\alpha}_{1},\dots,\widehat{\alpha}_{2q}$ and
$\widehat{w}_{1},\dots,\widehat{w}_{2q}$ are distinct and included in each
multiset with cardinality at most 1. Then
$|V(\kappa_{0})|=2q+\sum_{i=1}^{2q}a_{i}(\kappa_{0})$ where the $2q$ is from
the middle vertices $\widehat{w}_{1},\dots,\widehat{w}_{2q}$ and
$S_{\mathrm{iso}}$ is empty. By the definition of single-flipped graph matrix
shape (and taking $a_{i}:=a_{i}(\kappa_{0})$), the $i$th hyperedge has $\lceil
a_{i}/2+1/2\rceil$ or $\lfloor a_{i}/2\rfloor$ or $\lfloor a_{i}/2-1/2\rfloor$
vertices in $U_{\kappa_{0}}$ and $a_{i}-|e_{i}\cap U_{\kappa_{0}}|$ in
$V_{\kappa_{0}}$. Therefore $\min(|e_{i}\cap U_{\kappa_{0}}|,|e_{i}\cap
V_{\kappa_{0}}|)\geqslant\lfloor a_{i}/2-1/2\rfloor$ and so since a vertex
separator must separate $e_{i}\cap U_{\kappa_{0}}$ from $e_{i}\cap
V_{\kappa_{0}}$, at least $\lfloor a_{i}/2-1/2\rfloor$ of the vertices in
$e_{i}$ are part of any minimum vertex separator $S_{\min}$. Furthermore, if
$|e_{i}\cap U_{\kappa_{0}}|\neq|e_{i}\cap V_{\kappa_{0}}|\neq a_{i}/2$, then
the next hyperedge $e_{j}$ in the cycle that has an imbalance in $|e_{i}\cap
U_{\kappa_{0}}|$ vs $|e_{i}\cap V_{\kappa_{0}}|$ has that imbalance in the
opposite direction, forcing either a middle vertex in $W_{\kappa_{0}}$ to also
be part of this minimum vertex separator or for the vertex separator to stay
on the same side (either all $U_{\kappa_{0}}$ or all $V_{\kappa_{0}}$). The
latter will be cheaper, at the average cost of up to $1/2$ of a vertex per
hyperedge (except for the last one, if the total number of modes is odd).
Therefore, $|S_{\min}|\geqslant\lfloor\sum_{i}a_{i}/2\rfloor$. Overall,
$|V(\kappa_{0})|-|S_{\min}|+|W_{\mathrm{iso}}|\leqslant
2q+\lceil\sum_{i\in[2q]}a_{i}/2\rceil$, which is the bound we wanted.

The remaining shapes are formed from these extreme $\kappa_{0}$ by merging
vertices together and removing even numbers of duplicate copies of identical
hyperedges. Merging two vertices together has the effect of reducing
$|V(\kappa)|$ by 1, and it cannot reduce $S_{\min}$ by more than 1, so the
merging operation itself cannot cause
$|V(\kappa)|-|S_{\min}|+|W_{\mathrm{iso}}|$ to increase. The identical-
hyperedge-removing operation can cause $|W_{\mathrm{iso}}|$ to increase
without affecting $|V(\kappa)|$ or $|S_{\min}|$, but we argue that for each
element of $W_{\mathrm{iso}}$, there must have been at least one merging
operation that decreased $|V(\kappa)|$ without changing $|S_{\min}|$. Since
our construction of a minimum vertex separator did not use any vertices in
$W$, any middle vertex that is merged would decrease $|V(\kappa)|$ but not
decrease $|S_{\min}|$.

We can lower-bound the number of vertices in $W$ that were merged by looking
at the connected components of the graph when it is restricted to only the
vertices in $W$. Any merging/removal of hyperedges in the original graph
matrix shape imply the merging/removing of the corresponding edges in this
restricted graph $(V,E)$ with $V=W$. Each removal of a pair of identical edges
can only increase the number $C$ of connected components in this restricted
graph by $1$, so $C-1\leqslant q-|E|/2$. Unless there are no edges at all in
the graph (which is impossible by Item 7), $|W_{\mathrm{iso}}|\leqslant C-1$.
And finally, since there is at least one component with edges and every vertex
has even degree so at least one component is not a tree,
$|E|\geqslant|V|-C+1$. Therefore, $|V|\leqslant|E|+C-1\leqslant
2q-C+1\leqslant 2q-|W_{\mathrm{iso}}|$. Therefore, at least one vertex of $W$
must have been merged for every element of $|W_{\mathrm{iso}}|$ and the
inequality $|V(\kappa)|-|S_{\min}|+|W_{\mathrm{iso}}|\leqslant
2q+\lceil\sum_{i\in[2q]}a_{i}/2\rceil$ holds for all shapes $\kappa\in K$.

By Theorem 3.51, we then have with probability at least
$1-\varepsilon^{\prime}$ for each $\kappa$,

$\left\lVert M_{\kappa}\right\rVert_{{}_{\mathsf{op}}}\leqslant
O(1)n^{q+\lceil\sum_{i\in[2q]}a_{i}(\kappa)/2\rceil/2}\left(O(1)\left(|S_{\min}|\log
n-\log\varepsilon^{\prime}\right)\right)^{2qd_{H}}\,.$

By a union bound, Equation 6.5, and taking
$\varepsilon^{\prime}=e^{-n^{\lambda}}$, we have with probability
$1-O(e^{n^{\lambda}})$

$\displaystyle\operatorname*{\mathbb{E}}_{\sigma}\mathsf{noise}^{2p}$
$\displaystyle\leqslant O(1)\sum_{\kappa\in
K}n^{p-2p\lceil\sum_{i\in[2q]}a_{i}(\kappa)/2\rceil/2}\left\lVert\frac{1}{n^{q}}M_{\kappa}\right\rVert_{{}_{\mathsf{op}}}^{2p}$
$\displaystyle\leqslant O(n^{p+4pqd_{H}\lambda})\,.$

∎

###  6.2 Certifiable Taylor expansion

We state an elementary combinatorial lemma using Hölder’s inequality that
transparently shows that Lemma 6.1 implies an upper bound on the certified
value of the hessian under the maximum entropy constraint.

######  Lemma 6.2 (Hölder’s inequality for Hessian certificate).

Given the certificate in Lemma 6.1 the following holds by a sum-of-squares
proof:

$\displaystyle\operatorname*{\mathbb{E}}_{\sigma}\left(\frac{1}{k}\left\langle\frac{1}{n}\nabla^{2}H(\sigma),\operatorname*{\mathbb{E}}_{v}\left[vv^{\mathsf{T}}\mid\sigma\right]\right\rangle\right)^{2pq(p-1)}\leqslant\left((1-\varepsilon+O(\varepsilon^{2}))2\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{1/2}\frac{1}{k}\right)^{2pq(p-1)}+o_{n}(1)\,.$

###### Proof.

The proof follows by a simple combinatorial argument to choose $p$ large
enough, so that $C_{p}\approx(1-\varepsilon+O(\varepsilon^{2}))$, followed by
an application of Hölder’s inequality. First, note that,

$\displaystyle(C_{p/2}\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{p/2})^{2q(p-1)}=\left(\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})C_{p/2}^{1/p}\right)^{2pq(p-1)}$
$\displaystyle\sim_{\text{\lx@cref{creftypecap~refnum}{lem:catalan-
stirling}}}\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{2pq(p-1)/2}\left(\left(\frac{4^{p/2}}{\sqrt{\pi}\left(\frac{p}{2}\right)^{3/2}}\right)^{1/p}\right)^{2pq(p-1)}$
$\displaystyle=\left(2\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{1/2}\right)^{2pq(p-1)}\left(\left(\frac{2^{3}}{p^{3}(\pi)}\right)^{1/2p}\right)^{2pq(p-1)}$

By the maximum entropy constraint,

$\left\lVert\operatorname*{\mathbb{E}}_{v}[vv^{\mathsf{T}}\mid\sigma]\right\rVert^{q}_{q}\leqslant\left(\frac{1+o_{n}(1)}{\delta
n}\right)^{q-1}\,,$

and using Hölder’s inequality in conjunction with Lemma C.10 then yields,

$\displaystyle\operatorname*{\mathbb{E}}_{\sigma}\left(\frac{1}{k}\left\langle\frac{1}{n}\nabla^{2}H(\sigma),\operatorname*{\mathbb{E}}_{v}\left[vv^{\mathsf{T}}\mid\sigma\right]\right\rangle\right)^{2pq(p-1)}=\frac{1}{k^{2pq(p-1)}}\operatorname*{\mathbb{E}}_{\sigma}\left(\left[\left\langle\frac{1}{n}\nabla^{2}H(\sigma),\operatorname*{\mathbb{E}}_{v}\left[vv^{\mathsf{T}}\mid\sigma\right]\right\rangle\right]\right)^{2pq(p-1)}$
$\displaystyle\leqslant$
$\displaystyle\frac{1}{k^{2pq(p-1)}}\operatorname*{\mathbb{E}}_{\sigma}\left\lVert\frac{1}{n}\nabla^{2}H(\sigma)\right\rVert^{2pq(p-1)}_{p}\left\lVert\operatorname*{\mathbb{E}}_{v}\left[vv^{\mathsf{T}}\mid\sigma\right]\right\rVert^{2pq(p-1)}_{q}$
$\displaystyle\leqslant$
$\displaystyle\frac{1}{k^{2pq(p-1)}}\operatorname*{\mathbb{E}}_{\sigma}\left\lVert\frac{1}{n}\nabla^{2}H(\sigma)\right\rVert^{2pq(p-1)}_{p}\left(\frac{1+o_{n}(1)}{\delta
n}\right)^{2p(q-1)(p-1)}$ $\displaystyle\leqslant$
$\displaystyle\frac{1}{k^{2pq(p-1)}}\left((n\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{p/2}C_{p/2})^{2q(p-1)}+o(n^{2q(p-1)})\right)\left(\frac{1+o_{n}(1)}{\delta
n}\right)^{2p(q-1)(p-1)}$ $\displaystyle\leqslant$
$\displaystyle\frac{1}{k^{2pq(p-1)}}\frac{n^{2q+2p-2q}}{n^{2pq-2q}}\left(\frac{\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{pq(p-1)}(C_{p/2})^{2q(p-1)}}{\delta^{2p(q-1)(p-1)}}\right)+o_{n}(1)$
$\displaystyle=$
${}_{\frac{1}{p}+\frac{1}{q}=1}\frac{1}{k^{2pq(p-1)}}\left(2\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{1/2}\right)^{2pq(p-1)}\left(\frac{1}{(\delta\sqrt{\pi})^{1/p}(p/2)^{3/(2p)}}\right)^{2pq(p-1)}+o_{n}(1)$

We first invoke a series argument (ignoring contributions second-order and
higher as $p$ will be large) to expand the denominator terms as,

$\displaystyle(\delta\sqrt{\pi})^{1/p}=1+\frac{\log(\delta\sqrt{\pi})}{p}+O\left(\frac{1}{p^{2}}\right)\,,$
$\displaystyle(p/2)^{3/2p}=1+\frac{3\log(p/2)}{2p}+O\left(\frac{1}{p^{2}}\right)\,.$

We would like the approximation term to be
$1-\varepsilon+O(\varepsilon)^{2}\approx\frac{1}{1+\varepsilon+O(\varepsilon^{2})}$.
Equating the denominator (and ignoring contributions larger than second-
order),

$\displaystyle
1+\varepsilon+O\left(\frac{1}{p^{2}}\right)=1+\frac{3\log(p/2)}{2p}+\frac{\log(\delta\sqrt{\pi})}{p}+O\left(\frac{\log(p/2)\log(\delta\sqrt{\pi})}{p^{2}}\right)\,.$

The first-order solution to the above equation is given by the so-called
Lambert-W function [weisstein2002lambert] and is,

$p(\varepsilon,\delta)=\frac{-3W\left(\frac{-4}{3\pi^{1/3}}\cdot\frac{\varepsilon}{\delta^{2/3}}\right)}{2\varepsilon}\,.$

The choice of $\delta=A(\varepsilon)=O(\varepsilon)$ comes from Corollary 6.7
yielding a final value for $p$ as,

$p(\varepsilon)=-\frac{3W\left(\frac{-4}{3\pi^{1/3}}\cdot
O(\varepsilon^{1/3})\right)}{2\varepsilon}\,.$

Note that as $\varepsilon<<1$ we can do a Taylor expansion of the Lambert-W
function around $0$,

$W_{0}(z)=\sum_{n=1}^{\infty}\frac{(-n)^{n-1}}{n!}z^{n}=z-z^{2}+\frac{3}{2}z^{3}-\frac{8}{3}z^{4}+\dots\,,$

yielding,

$p(\varepsilon)=-\frac{3}{2\varepsilon}W\left(\frac{-4}{3\pi^{1/3}}\cdot
O(\varepsilon^{1/3})\right)=\frac{2}{\pi^{1/3}}O(\varepsilon^{-2/3})-O(\varepsilon^{-1/3})+O(1)\,.$

Close to $0^{-}$, the Lambert function is linear in the first-order and,
therefore, monotone. This implies that one can scale by a constant factor
(smaller than $O(\varepsilon^{-2/3})$) so that $p(\varepsilon)$ is an integer
and assert that,

$\frac{2}{\pi^{1/3}}O(\varepsilon^{-2/3})-O(\varepsilon^{-1/3})+O(1)=p(\varepsilon)\geqslant
2\,,$

for a sufficiently small $\varepsilon$. This allows for the conclusion that
the Hölder conjugate $q\geqslant 1$, and that,

$\displaystyle\left(\frac{1}{(\delta\sqrt{q})^{1/p}(p/2)^{3/2p}}\right)^{2pq}$
$\displaystyle=1-\varepsilon+O(\varepsilon^{2})\implies$
$\displaystyle\operatorname*{\mathbb{E}}_{\sigma}\left(\frac{1}{k}\left\langle\frac{1}{n}\nabla^{2}H(\sigma),\operatorname*{\mathbb{E}}_{v}\left[vv^{\mathsf{T}}\mid\sigma\right]\right\rangle\right)^{2pq}\leqslant\left((1-O(\varepsilon))2\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{1/2}\frac{1}{k}\right)^{2pq}+o_{n}(1)\,.$

∎

While it is possible to achieve the upper bound above with _some_ choice of
$\operatorname*{\mathbb{E}}[vv^{\mathsf{T}}\mid\sigma]$, as our constraint
system restricts this to be a PSD matrix, we must write a small modification
of the lemma above to work for an “adjustment” of the hessian (where it is
shifted to be PSD with high probability).

######  Corollary 6.3 (Large correlation of conditional covariance with
“shifted” Hessian).

Given the certificate in Lemma 6.1, the following holds by a sum-of-squares
proof:

$\operatorname*{\mathbb{E}}_{\sigma}\left(\frac{1}{k}\left\langle\left(2\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{1/2}\operatorname{\mathrm{Id}}\right)+\frac{1}{n}\nabla^{2}H(\sigma),\operatorname*{\mathbb{E}}[vv^{\mathsf{T}}\mid\sigma]\right\rangle\right)\leqslant\frac{1}{k}2\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{1/2}\left(2-\varepsilon+O(\varepsilon^{2})\right)+o_{n}(1)\,.$

###### Proof.

We do a binomial expansion on the shifted Hessian and use one more application
of Stirling’s approximation of the factorial to get a uniform upper bound on
the Catalan numbers, leading to a _slightly_ different parameterization of the
Lambert-W function to obtain the right choice of
$(p(\varepsilon),q(\varepsilon))$-Hölder conjugates that correspond to the
desired approximation.

We simplify the proof used above by using the radical
$A=\operatorname*{\mathbb{E}}[vv^{\mathsf{T}}\mid\sigma]^{1/(p-1)}$ and using
it with an application of Lemma D.4 (with
$B=\frac{1}{n}\nabla^{2}H(\sigma)+2\nu^{\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{1/2}\operatorname{\mathrm{Id}}_{n}$)
and obtain,

$\displaystyle\left(\operatorname*{\mathbb{E}}_{\sigma}\left(\left\langle\left(2\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{1/2}\operatorname{\mathrm{Id}}\right)+\frac{1}{n}\nabla^{2}H(\sigma),\operatorname*{\mathbb{E}}[vv^{\mathsf{T}}\mid\sigma]\right\rangle\right)\right)^{p}$
$\displaystyle\qquad\qquad\leqslant_{\text{Jensen with
$p$=even}}\operatorname*{\mathbb{E}}_{\sigma}\left(\operatorname{Tr}[BA^{p-1}]\right)^{p}\leqslant_{\text{\lx@cref{creftypecap~refnum}{lem:matrix-
p-q-
holders}}}\left(\operatorname{Tr}[A^{p}]\right)^{p-1}\operatorname*{\mathbb{E}}_{\sigma}\operatorname{Tr}[B^{p}]$
$\displaystyle=(\operatorname{Tr}[\operatorname*{\mathbb{E}}[vv^{\mathsf{T}}\mid\sigma]^{p/(p-1)}])^{p-1}\operatorname*{\mathbb{E}}_{\sigma}\operatorname{Tr}[B^{p}]=_{q=p/(p-1)}\left(\operatorname{Tr}[\operatorname*{\mathbb{E}}[vv^{\mathsf{T}}\mid\sigma]^{q}]\right)^{p-1}\operatorname*{\mathbb{E}}_{\sigma}\operatorname{Tr}[B^{p}]$
$\displaystyle=\left(\left\lVert\operatorname*{\mathbb{E}}[vv^{\mathsf{T}}\mid\sigma]\right\rVert_{q}^{q}\right)^{p-1}\operatorname*{\mathbb{E}}_{\sigma}\operatorname{Tr}[B^{p}]\leqslant_{\text{\lx@cref{creftypecap~refnum}{prop:q-norm-
cond-cov}}}\left(\frac{1}{(\delta
n)^{q-1}}\right)^{p-1}\operatorname*{\mathbb{E}}_{\sigma}\operatorname{Tr}[B^{p}]=_{1/p+1/q=1}\frac{1}{\delta
n}\operatorname*{\mathbb{E}}_{\sigma}\operatorname{Tr}[B^{p}]$
$\displaystyle\leqslant_{\text{Binomial expansion}}\frac{1}{\delta
n}\left(\sum_{k=2,\text{even}}^{p}\binom{p}{k}(2\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{1/2})^{p-k}\operatorname*{\mathbb{E}}_{\sigma}\operatorname{Tr}\left[\left(\frac{1}{n}\nabla^{2}H(\sigma)\right)^{k}\right]+(2\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{1/2})^{p}n+o(n)\right)$
$\displaystyle\leqslant_{\text{\lx@cref{creftypecap~refnum}{lem:hessian-
schatten-norm}}}\frac{1}{\delta
n}\left(\sum_{k\in\\{2,4,\dots,p\\}}\binom{p}{k}(2\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{1/2})^{p-k}\left(\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{k/2}C_{k/2}+o(n)\right)+(2\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{1/2})^{p}n+o(n)\right)$
$\displaystyle=(2\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{1/2})^{p}\frac{1}{\delta}\left(\left\\{\sum_{k\in\\{2,4,\dots,p\\}}\binom{p}{k}2^{-k}C_{k/2}\right\\}+1\right)+o_{n}(1)$
$\displaystyle=_{\text{\lx@cref{creftypecap~refnum}{lem:catalan-
stirling}}}(2\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{1/2})^{p}\frac{1}{\delta}\left(\left\\{{\frac{2\sqrt{2}}{\sqrt{\pi}}}\sum_{k\in\\{2,4,\dots,p\\}}\binom{p}{k}\frac{1}{k^{3/2}}\left(1+O\left(\frac{1}{k}\right)\right)\right\\}+1\right)+o_{n}(1)\,,$

where we have used semi-algebraic degree reduction (one-sided) of the form as
shown in Lemma C.7 in the calculation above.
By a calculation similar to [spencer2014asymptopia, Section 5.4], using
estimates similar to those in [spencer2014asymptopia, Eq. 5.45 Eq. 5.49 & Eq.
5.50] we have that, for $k^{\prime}=\frac{p+100\sqrt{\ln(p)}}{2}$,

$\sum_{s\in\\{k^{\prime},p\\}}\binom{p}{s}\frac{1}{s^{3/2}}\leqslant\Theta\left(\frac{2^{p}}{p^{5001}}\right)\,,$

using the fact that for the range of $s$ above, $\frac{1}{s^{3/2}}\leqslant
O\left(\frac{1}{p^{3/2}}\right)$. Similarly, by symmetry,

$\sum_{s\in\\{0,p-k^{\prime}\\}}\binom{p}{s}\frac{1}{s^{3/2}}\leqslant\Theta\left(\frac{2^{p}}{p^{4999}}\right)\,,$

with the bound of $\frac{1}{s^{3/2}}\leqslant 1$ for the desired range. Since
the largest contributions come from the middle terms in the binomial sum,
those are bounded as,

$\displaystyle\sum_{s\in\\{p-k^{\prime},k^{\prime}\\}}\binom{p}{s}\frac{1}{s^{3/2}}$
$\displaystyle\leqslant\left(\frac{2\sqrt{2}}{p^{3/2}}+o(1/p^{3/2})\right)\sum_{s\in\\{p-k^{\prime},k^{\prime}\\}}\binom{p}{k}\sim\frac{2\sqrt{2}}{p^{3/2}}\binom{p}{p/2}\sum_{s=\frac{p+c\sqrt{\ln(p)}}{2}\in\\{k^{\prime},p-k^{\prime}\\}}e^{-(c^{2}/2)}$
$\displaystyle\leqslant\frac{2\sqrt{2}}{p^{3/2}}\binom{p}{p/2}\sqrt{p}\int_{-\infty}^{\infty}\frac{e^{-c^{2}/2}}{2}dc=\frac{2\sqrt{2}}{p^{3/2}}\binom{p}{p/2}\sqrt{p}\sqrt{\frac{\pi}{2}}\,$

where we used the fact that for $s\in\\{p-k^{\prime},k^{\prime}\\}$ it is the
case that,

$\displaystyle\frac{1}{s^{3/2}}$
$\displaystyle\leqslant\frac{1}{\left(\frac{p-100\sqrt{\ln(p)}}{2}\right)^{3/2}}=\frac{2\sqrt{2}}{(p-100\sqrt{\ln(p)})^{3/2}}$
$\displaystyle=\frac{2\sqrt{2}}{p^{3/2}\left(1-\frac{100\sqrt{\ln(p)}}{p^{3/2}}\right)}=\frac{2\sqrt{2}}{p^{3/2}}\left(1+O\left(\frac{\sqrt{\ln(p)}}{p^{3/2}}\right)\right)$
$\displaystyle=\frac{2\sqrt{2}}{p^{3/2}}+O\left(\frac{\sqrt{\ln(p)}}{p^{3}}\right)\,,$

for every $p$. Using the Stirling approximation for $\binom{p}{p/2}$ again
leads to,

$\binom{p}{p/2}\sqrt{p}\sqrt{\frac{\pi}{2}}\sim\frac{2^{p}}{\sqrt{p}}\sqrt{\frac{2}{\pi}}\sqrt{p}\sqrt{\frac{\pi}{2}}=2^{p}\,.$

A similar analysis can be used to show that the terms above strictly dominate
the remaining $o\left(\frac{1}{k^{3/2}}\right)$ terms, which continue to
contribute on the same order and $\sim$ is used to denote precise equality _up
to first order_ terms in the approximation.
Putting the above bounds altogether leads to the final bounds on the inner-
product as,

$\displaystyle\left(\operatorname*{\mathbb{E}}_{\sigma}\left(\left\langle\left(2\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{1/2}\operatorname{\mathrm{Id}}\right)+\frac{1}{n}\nabla^{2}H(\sigma),\operatorname*{\mathbb{E}}[vv^{\mathsf{T}}\mid\sigma]\right\rangle\right)\right)^{p}\leqslant(2\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{1/2})^{p}\frac{1}{\delta}\left(\frac{8}{\sqrt{\pi}}\frac{2^{p}}{p^{3/2}}+1\right)+o_{n}(1)$
$\displaystyle\implies_{\text{\lx@cref{creftypecap~refnum}{lem:ideal-
reduction-2p}}}\operatorname*{\mathbb{E}}_{\sigma}\left\langle\left(2\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{1/2}\operatorname{\mathrm{Id}}\right)+\frac{1}{n}\nabla^{2}H(\sigma),\operatorname*{\mathbb{E}}[vv^{\mathsf{T}}\mid\sigma]\right\rangle$
$\displaystyle\qquad\qquad\qquad\leqslant
2\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{1/2}\left(\frac{1}{\delta}\left(\frac{8}{\sqrt{\pi}}\frac{2^{p}}{p^{3/2}}+1\right)+o_{n}(1)\right)^{1/p}$
$\displaystyle\qquad\qquad\qquad\leqslant
2\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{1/2}\left(\frac{1}{\delta}\left(\frac{8}{\sqrt{\pi}}\frac{2^{p}}{p^{3/2}}+1\right)\right)^{1/p}+o_{n}(1)\,,$

and the goal is to set $p(\varepsilon)$ so that,

$\left(\frac{1}{\delta}\left(\frac{8}{\sqrt{\pi}}\frac{2^{p}}{p^{3/2}}+1\right)\right)^{1/p}=(2-\varepsilon+O(\varepsilon)^{2})\,.$

Note that,

$\displaystyle\left(\frac{1}{\delta}\right)^{1/p}\left(1+\left(\frac{2C^{1/p}}{p^{3/2p}}\right)^{p}\right)^{1/p}\leqslant\left(\frac{1}{\delta}\right)^{1/p}\left(2\left(\frac{2C^{1/p}}{p^{3/2p}}\right)^{p}\right)^{1/p}=2\left(\frac{2}{\delta}\right)^{1/p}\left(\frac{C}{p^{3/2}}\right)^{1/p}=2\left(\frac{(2C)^{1/p}}{\delta^{1/p}p^{3/2p}}\right)\,,$

and so the final desired ratio is,

$1+\frac{\varepsilon}{2}+O(\varepsilon^{2})=(\delta)^{1/p}\left(\frac{p}{2C}\right)^{3/2p}\,,$

at which point an analysis for the Lambert-W function identical to the one in
Lemma 6.2 gives the desired choice of $p(\varepsilon)=O(\varepsilon^{-2/3})$.
∎

######  Lemma 6.4 (Taylor Series expansion under maximum-entropy).

Given a program iterate $\sigma$ and increment $v$ from a polynomial grammar
satisfying the axioms of the HES SoS relaxation (up to an appropriate constant
degree), the following statement holds by a SoS proof:

$\operatorname*{\mathbb{E}}_{v,\sigma}\left(\frac{1}{n}\left(H\left(\sigma+\frac{v}{\sqrt{k}}\right)-H(\sigma)\right)\right)\leqslant(1-O(\varepsilon)-o_{n}(1))\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{1/2}\frac{1}{k}+O\left(\frac{1}{k^{3}}\right)+o_{n}(1)\,,$
(6.6)

with high probability over the randomness of the input.

###### Proof.

The proof for this follows by combining the following three facts:

  1. 1.

The low-degree $p$-Schatten norm SoS certificate of the hessian
$\nabla^{2}H(\sigma)$ with an appropriately large choice of $p$ (Lemma 6.1).
This conjoined with a Holder’s inequality yields the desired bound.

  2. 2.

The gradient term is zero since the first cumulant is zero.

  3. 3.

The low-degree SoS certificates that, for all $p\in[3,d_{H}]$, the
contributions of
$\frac{1}{j!k^{j/2}n}\langle\nabla^{r}H(\sigma_{j}),v^{\otimes r}\rangle$
(Theorem 5.9) are $\frac{O(\varepsilon^{2})}{k}$. In particular, when
$\left\lVert P\right\rVert_{{}_{\mathsf{op}}}=n^{\lceil d_{p}/2+1\rceil/2}$
(as is the case when each of its coefficients/entries is random, independent,
and on the order of $\sqrt{n}$), $k=1$ and $\nu=\frac{1}{\delta n}$, we have

$\operatorname*{\mathbb{E}}_{\sigma,u}p(\sigma,u)\leqslant
O(d_{p})^{O(d_{p}+i+d)}(\nu n)^{d_{p}/2}n\,.$

This immediately implies, with the appropriate choice of $\nu$ (and a choice
of step precision $\varepsilon$ that is _exponentially_ small in the desired
global precision) that,
$\frac{1}{n}\operatorname*{\mathbb{E}}_{\sigma,u}[p(\sigma,u)]\leqslant\frac{O(\varepsilon^{2})}{k}$.

The Taylor expansion of the $H(\sigma+\frac{v}{\sqrt{k}})$ around $\sigma$ is
given as,

$\displaystyle\frac{1}{n}H\left(\sigma+\frac{v}{\sqrt{k}}\right)=\frac{1}{n}H(\sigma)+\frac{1}{\sqrt{k}n}\left\langle\nabla
H(\sigma),v\right\rangle+\frac{1}{2kn}\left\langle\nabla^{2}H(\sigma),vv^{\mathsf{T}}\right\rangle+\sum_{r=3}^{d_{H}}\frac{1}{k^{r/2}r!n}\langle\nabla^{r}H(\sigma),v^{\otimes
r}\rangle\,.$

Applying the conditional expectation operators as well,

$\displaystyle\operatorname*{\mathbb{E}}_{\sigma,v}\left(\frac{1}{n}H\left(\sigma+\frac{v}{\sqrt{k}}\right)-\frac{1}{n}H(\sigma)\right)=\sum_{j=2}^{d_{H}}\frac{1}{j!k^{j/2}n}\operatorname*{\widetilde{\mathbb{E}}}_{\sigma,v}\langle\nabla^{r}H(\sigma),v^{\otimes
r}\rangle\,,$

where we used the fact that the gradient is zero. To isolate the contribution
of the derivatives, we will apply Lemma C.7 in conjunction with a pseudo-
Jensen inequality for even powers. Since there are only $k$ steps, only the
hessian term will survive while the others will contribute no more than
$O(1/k^{3/2})$ per step.

$\displaystyle\sum_{j=2}^{d_{H}}\frac{1}{j!k^{j/2}n}\operatorname*{\mathbb{E}}_{\sigma,v}\langle\nabla^{r}H(\sigma),v^{\otimes
r}\rangle$
$\displaystyle=\frac{1}{2k}\operatorname*{\mathbb{E}}_{\sigma}\left\langle\frac{1}{n}\nabla^{2}H(\sigma),\operatorname*{\mathbb{E}}_{v}\left[vv^{\mathsf{T}}\mid\sigma\right]\right\rangle+O\left(\frac{1}{k^{3/2}}\right)+o_{n}(1)\,.$

We now bound the hessian term.

$\displaystyle\left(\frac{1}{k}\operatorname*{\mathbb{E}}_{\sigma}\left[\left\langle\frac{1}{n}\nabla^{2}H(\sigma),\operatorname*{\mathbb{E}}_{v}\left[vv^{\mathsf{T}}\mid\sigma\right]\right\rangle\right]\right)^{2pq(p-1)}$
$\displaystyle\leqslant_{\text{\lx@cref{creftypecap~refnum}{lem:pseudo-jensen-
even}}}\operatorname*{\mathbb{E}}_{\sigma}\left(\frac{1}{k}\left\langle\frac{1}{n}\nabla^{2}H(\sigma),\operatorname*{\mathbb{E}}_{v}\left[vv^{\mathsf{T}}\mid\sigma\right]\right\rangle\right)^{2pq(p-1)}$
$\displaystyle\leqslant\left((1-O(\varepsilon))2\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{1/2}\frac{1}{k}\right)^{2pq(p-1)}+o_{n}(1)\,.$

Applying Lemma C.7 yields,

$\displaystyle\frac{1}{k}\left(\operatorname*{\mathbb{E}}_{\sigma}\left[\left\langle\frac{1}{n}\nabla^{2}H(\sigma),\operatorname*{\mathbb{E}}_{v}\left[vv^{\mathsf{T}}\mid\sigma\right]\right\rangle\right]\right)\leqslant\left(\left((1-O(\varepsilon))2\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{1/2}\frac{1}{k}\right)^{2pq(p-1)}+o_{n}(1)\right)^{1/(2pq(p-1))}$
$\displaystyle=(1-O(\varepsilon))2\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{1/2}\frac{1}{k}\left(1+\frac{o_{n}(1)}{\left((1-O(\varepsilon))2\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{1/2}\frac{1}{k}\right)^{2pq(p-1)}}\right)^{1/(2pq(p-1))}$
$\displaystyle=_{\text{Taylor
expansion}}(1-O(\varepsilon))2\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{1/2}\frac{1}{k}\left(1+\frac{o_{n}(1)}{2pq(-1)\left((1-O(\varepsilon))2\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{1/2}\frac{1}{k}\right)}+o_{n}(1)\right)$
$\displaystyle=(1-O(\varepsilon))2\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{1/2}\frac{1}{k}\left(1+o_{n}(1)\right)=\left(1-O(\varepsilon)-o_{n}(1)\right)2\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{1/2}\frac{1}{k}\,.$

A similar calculation as above for the higher-order derivative certificates
provided by Theorem 5.9 concludes the argument. ∎

We demonstrate that the global optimization stated in the axioms of the
program (which is feasible as shown in Section 7.1) in conjunction with the
upper bound on the (pseudo-expected) increase of the energy at every step
concludes a lower bound as well.

######  Lemma 6.5 (Local maximization via maximum-entropy and global
maximization).

Given the maximum entropy constraints and the (feasible) constraint of global
(hamiltonian) maximization, the following lower bound holds by a SoS proof:

$\operatorname*{\mathbb{E}}_{k,\dots,i}\left[H\left(\sum_{j=1}^{i}\frac{v_{j}}{\sqrt{k}}\right)-H\left(\sum_{j=1}^{i-1}\frac{v_{j}}{\sqrt{k}}\right)\right]\geqslant\frac{1}{k}(1-\varepsilon)\nu^{\prime\prime}(i/k)^{1/2}-O(\varepsilon^{2})-O\left(\frac{1}{k^{3/2}}\right)-o_{n}(1)\,.$

###### Proof.

Note that the bounds provided in Proposition C.9 are dual sided, and the
certificate for the hessian in Lemma 6.1 controls the deviation. Therefore,
combining Proposition C.9 and Lemma 6.1 immediately yields that,

$n\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{k/2}C_{k/2}-o(n)\leqslant\operatorname*{\mathbb{E}}_{\sigma}\left\lVert\frac{1}{n}\nabla^{2}H(\sigma)\right\rVert^{k}_{k}\leqslant
n\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{k/2}C_{k/2}+o(n)\,.$

Given the assumption of global optimality (under moments) which is made
possible because of feasibility, that is,

$\operatorname*{\mathbb{E}}\left[H\left(\sum_{i=1}^{k}\frac{v_{i}}{\sqrt{k}}\right)\right]\geqslant(1-\varepsilon)\max_{\sigma\in\mathcal{S}^{n-1}(1)}H(\sigma)\sim_{\text{by~{}\lx@cref{creftypecap~refnum}{lem:taylor}}}(1-\varepsilon)\int_{0}^{1}\nu^{\prime\prime}(q)^{1/2}dq+O(\varepsilon)^{2}+o_{n}(1)\,,$

and the bound on the higher-order derivatives certified by sum-of-squares to
be $O(1/k^{3/2})+o_{n}(1)$ (Theorem 5.9) allows one to conclude by simple
linear manipulation that,

$\operatorname*{\mathbb{E}}\left[H\left(\sum_{j=1}^{k}\frac{v_{j}}{\sqrt{k}}\right)\right]-\operatorname*{\mathbb{E}}\left[H\left(\sum_{j=1}^{k-1}\frac{v_{j}}{\sqrt{k}}\right)\right]\geqslant\frac{1}{k}(1-\varepsilon)\left(\nu^{\prime\prime}(1)^{1/2}\right)-O(\varepsilon^{2})-O\left(\frac{1}{k^{3/2}}\right)-o_{n}(1)\,,$

where the contribution is as is because one subtracts all steps but the last
one to get the relevant contribution from the Hessian. ∎

###  6.3 Matrix Bernstein approximation bound via analytic sparsity

Note that the _true_ Hessian is distributed as
$\nu^{\prime\prime}(\left\lVert\sigma\right\rVert)^{1/2}\mathsf{GOE}(n)$ with
high probability. Consequently, it is traceless. However, in order to analyze
the rounding scheme introduced in [ss24rounding, Section 2], it is critical
show that the Bernstein polynomial has trace $\approx 1$. To do this, we need
to take the Bernstein polynomial of the _shifted_ Hessian
$2\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{1/2}\operatorname{\mathrm{Id}}_{n}+\frac{1}{n}\nabla^{2}H(\sigma)=p\left(\frac{1}{n}\nabla^{2}H(\sigma)\right)$
and invoke Lemma 6.6 (proved below). This allows one to conclude that, for an
appropriate choice of the ramp (fixed in Corollary 6.7), the trace of the
(appropriately) normalized (divided by $1/(\delta n)$) Bernstein polynomial
can be made $\approx 1$.

######  Lemma 6.6 (Equivalence with semicircle moments).

Let
$X(\sigma)=M(\sigma)/(2\nu^{\prime\prime}(\left\lVert\sigma\right\rVert_{2}^{2})^{1/2})$
for a symmetric matrix $M(\sigma)$ with entries that are low-degree
polynomials in $\sigma$. Let $\mathcal{D}$ be a pseudo-distribution over high-
entropy step distributions over $\sigma$ satisfying the conclusions of Lemma
6.1 and let $p$ be a degree-$d$ polynomial. Then, for all positive integer
$r$,

$\operatorname*{\mathbb{E}}_{\sigma}\left(\left[\operatorname{Tr}\frac{1}{n}p(X(\sigma))\right]-n\operatorname*{\mathbb{E}}_{x\sim\mu_{sc}}\left[p(x)\right]\right)^{2r}\leqslant
o(n^{2r})\,\left\lVert p\right\rVert_{L^{2}}^{2r}$

where $\mu_{sc}$ is the Wigner semicircle law given by
$\mu_{sc}(x)=\frac{2}{\pi}\sqrt{1-x^{2}}$. This has an SoS proof.

###### Proof.

We write $p(\cdot)$ in the monomial basis, and then apply a Cauchy-Schwarz
inequality on the summands to separate the constants from the SoS program
variables (which are functions of $\sigma$). This is followed by an
application of Lemma 6.1 in concurrence with Lemma C.4. We will crucially use
the facts that the moment variable $r=O_{\varepsilon}(1)$ and that the degree
of the polynomial $d=O_{\varepsilon}(1)$.

First, we define

$\displaystyle
y_{i}(\sigma):=\operatorname{Tr}\left[\left(\frac{1}{n}X(\sigma)\right)^{i}\right]-n\operatorname*{\mathbb{E}}_{x\sim\mu_{sc}}\left[x^{i}\right]\,.$

Then, for all bounded $i$ and bounded even $j$,

$\operatorname*{\mathbb{E}}\left[y_{i}(\sigma)^{j}\right]\leqslant
o(n^{j})\,,$ (6.7)

which follows directly from Lemma 6.1.

Expressing $p(x)=\sum_{i=0}^{d}c_{i}x^{i}$ and using linearity of the trace
and expectation immediately gives

$\displaystyle\operatorname*{\mathbb{E}}\left[\left(\frac{1}{n}\operatorname{Tr}[p(X(\sigma))]-\operatorname*{\mathbb{E}}_{x\sim\mu_{sc}}[p(x)]\right)^{2r}\right]$
$\displaystyle=\operatorname*{\mathbb{E}}\left[\left(\sum_{i=0}^{d}c_{i}\operatorname{Tr}\left[(\frac{1}{n}X(\sigma))^{i}\right]-\sum_{i=0}^{d}c_{i}\operatorname*{\mathbb{E}}_{x\sim\mu_{sc}}[x^{i}]\right)^{2r}\right]$
$\displaystyle=\operatorname*{\mathbb{E}}\left[\left(\sum_{i}c_{i}y_{i}(\sigma)\right)^{2r}\right]\,.$

Applying the strategy of first isolating the terms containing the program
variables $\left(\\{y_{i}(\sigma))\\}_{i=0}^{d}\right)$ from the constant
terms that come from the coefficients $\\{c_{i}\\}_{i=0}^{d}$ via an
application of the Cauchy-Schwarz (CS) inequality yields,

$\displaystyle\operatorname*{\mathbb{E}}_{\sigma\sim\mathcal{D}}\left[\left(\sum_{i=0}^{d}c_{i}y_{i}(\sigma)\right)^{2r}\right]=\operatorname*{\mathbb{E}}_{\sigma\sim\mathcal{D}}\left[\left(\left\langle\mathbf{c},\mathbf{y(\sigma)}\right\rangle^{2}\right)^{r}\right]\overset{\mathsf{CS}}{\leqslant}\operatorname*{\widetilde{\mathbb{E}}}_{\sigma\sim\mathcal{D}}\left[\left(\sum_{i=0}^{d}c_{i}^{2}\right)^{r}\left(\sum_{i=0}^{d}y^{2}_{i}(\sigma)\right)^{r}\right]$
$\displaystyle=\left(\sum_{i=0}^{d}c_{i}^{2}\right)^{r}\operatorname*{\mathbb{E}}_{\sigma\sim\mathcal{D}}\left[\left(\sum_{i=0}^{d}y_{i}^{2}(\sigma)\right)^{r}\right]=\left\lVert
p\right\rVert^{2r}_{2}\operatorname*{\mathbb{E}}_{\sigma\sim\mathcal{D}}\left[\sum_{k_{0}+\dots+k_{d}=r}\binom{r}{k_{0},\dots,k_{d}}y^{2k_{0}}_{0}(\sigma)\cdots
y^{2k_{d}}_{d}(\sigma)\right]$
$\displaystyle\leqslant_{\text{\lx@cref{creftypecap~refnum}{lem:l1-to-
lt}}}\left\lVert
p\right\rVert^{2r}_{2}\operatorname*{\mathbb{E}}_{\sigma\sim\mathcal{D}}\left[d^{r-1}\sum_{i=0}^{d}y_{i}^{2r}\right]\leqslant\left\lVert
p\right\rVert_{2}^{2r}d^{r}o(n^{2r})=o(n^{2r})\left\lVert
p\right\rVert_{2}^{2r}\,.$ (6.8)

Since $d,r=O_{\varepsilon}(1)$, this yields that,

$\displaystyle\operatorname*{\mathbb{E}}\left[\left(\frac{1}{n}\operatorname{Tr}[p(X(\sigma))]-\operatorname*{\mathbb{E}}_{x\sim\mu_{sc}}[p(x)]\right)^{2r}\right]\leqslant\left\lVert
p\right\rVert^{2r}_{2}f(\varepsilon)o(n^{2r})\,,$

for some bounded function $f(\varepsilon)$, concluding the proof. ∎

Alternatively, a better dependency on the constants can be obtained using
Lemma C.12, however, it is irrelevant to the use case.

######  Corollary 6.7 (Large Hilbert-Schmidt Correlation between Bernstein
Polynomial and Hessian).

Given the Bernstein matrix polynomial $B^{L^{\alpha}_{\gamma}}$ as defined in
Proposition A.2, it is the case that for any high-entropy step distribution
$\mathcal{D}$,

$\operatorname*{\mathbb{E}}_{\sigma\sim\mathcal{D}}\left(\frac{1}{n}\left[\langle
B^{L^{\alpha}_{\gamma}}(X),X\rangle\right]-\psi\right)^{2p}\leqslant
o_{n}(1)\,,$ (6.9)

where $a\operatorname{\mathrm{Id}}_{n}\preceq X\preceq
b\operatorname{\mathrm{Id}}_{n}$, with $\psi\geqslant 1-\varepsilon$.

###### Proof.

To compute the correlation with the matrix Bernstein polynomial that
approximates the ramp function, we will use the definition of the Bernstein
polynomial followed by linearity (again) in conjunction with the bounds above.

$\displaystyle\frac{1}{n}\left[\left\langle
B^{L_{\alpha,\gamma}}_{d}(X),X\right\rangle\right]$
$\displaystyle=\frac{1}{n}\operatorname{Tr}\left(X\cdot
B^{L^{\alpha}_{\gamma}}_{d}(X)\right),$

By Lemma 6.6,

$\displaystyle\left(\frac{1}{\operatorname*{\mathbb{E}}_{x\sim\mu_{sc}}\left[B^{L^{\alpha}_{\gamma}}(x)\right]}\right)^{2p}\operatorname*{\mathbb{E}}_{\sigma\sim\mathcal{D}}\left(\frac{1}{n}\left[\langle
B^{L^{\alpha}_{\gamma}}(X),X\rangle\right]-\operatorname*{\mathbb{E}}_{x\sim\mu_{sc}}\left[x\cdot
B^{L^{\alpha}_{\gamma}}(x)\right]\right)^{2p}\leqslant\left(\frac{1}{\operatorname*{\mathbb{E}}_{x\sim\mu_{sc}}\left[B^{L^{\alpha}_{\gamma}}(x)\right]}\right)^{2p}o_{n}(1)\leqslant
o_{n}(1)\,,$

provided that
$\operatorname*{\mathbb{E}}_{x\sim\mu_{sc}}\left[B^{L^{\alpha}_{\gamma}}(x)\right]=C(\varepsilon)$
and $p=O_{\varepsilon}(1)$, where the former will follow from Lemma B.1 with
the appropriate choice of constants made later.
This reduces the goal to showing that,

$\frac{\operatorname*{\mathbb{E}}_{x\sim\mu_{sc}}\left[x\cdot
B^{L^{\alpha}_{\gamma}}(x)\right]}{\operatorname*{\mathbb{E}}_{x\sim\mu_{sc}}\left[B^{L^{\alpha}_{\gamma}}(x)\right]}\geqslant
1-\varepsilon\,,$ (6.10)

for an appropriate choice of $d(\varepsilon)$. Since this does not involve any
SoS program variable, it need not be a SoS proof.
By the definition in Appendix A and appropriately rescaling the interval
$\left[-2\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{1/2},2\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{1/2}\right]\to[-1,1]$
and denoting the re-scaled area of the semi-circle containing the valid
fraction of eigenvectors as

$\displaystyle\phi:=\frac{A(\varepsilon)}{2\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{1/2}}\,,$

we have the following,

$\displaystyle L^{\alpha}_{\gamma}(x)=\left\\{\begin{array}[]{l
l}0\,,&x\in\left[-1,\,1-\phi\right),\\\
\frac{1}{\phi^{2}}(x-(1-\phi))\,,&x\in\left[1-\phi,\,1-\phi+\phi^{2}\right],\\\
1\,,&x\in\left(1-\phi+\phi^{2},\,1\right],\end{array}\right.$ (6.14)

where $A(\varepsilon)$ is chosen later. The Weierstrass approximation (Theorem
3.29) gives,

$\displaystyle\|B^{L^{\alpha}_{\gamma}}(x)-L^{\alpha}_{\gamma}(x)\|\leqslant\varepsilon\,,\forall
x\in[-1,1]\,,$

provided $\mathsf{deg}(B^{L^{\alpha}_{\gamma}}(\cdot))\geqslant
O\left(\frac{1}{\varepsilon^{3}}\right)$.

##### Average correlation between $x$ and $B^{L^{\alpha}_{\gamma}}(x)$ under
$\mu_{sc}$

We now compute the scalar proxy for the correlation between the hessian and
the matrix Bernstein polynomial under the semi-circle measure. This will be
normalized appropriately using Lemma B.1.

$\displaystyle\operatorname*{\mathbb{E}}_{x\sim\mu_{sc}}\left[x\cdot
B^{L^{\alpha}_{\gamma}}(x)\right]=\int_{-1}^{1}x\cdot
B^{L^{\alpha}_{\gamma}}(x)\frac{2}{\pi}\sqrt{1-x^{2}}dx$
$\displaystyle\geqslant\int_{-1}^{1-\phi}x\left(L^{\alpha}_{\gamma}+\varepsilon\right)\frac{2}{\pi}\sqrt{1-x^{2}}dx+\int_{1-\phi}^{1}x\left(L^{\alpha}_{\gamma}(x)-\varepsilon\right)\frac{2}{\pi}\sqrt{1-x^{2}}dx$
$\displaystyle=\varepsilon\int_{-1}^{1-\phi}xd\mu_{sc}(x)-\varepsilon\int_{1-\phi}^{1}xd\mu_{sc}(x)+\int_{1-\phi}^{1-\phi+\phi^{2}}x\cdot\left(\frac{1}{\phi^{2}}(x-(1-\phi))\right)d\mu_{sc}(x)+\int_{1-\phi+\phi^{2}}^{1}xd\mu_{sc}(x)$
$\displaystyle=\underbrace{\varepsilon\left(\int_{-1}^{1-\phi}xd\mu_{sc}(x)-\int_{1-\phi}^{1}xd\mu_{sc}(x)\right)}_{\mathsf{I}(\varepsilon,\phi)}+\underbrace{\int_{1-\phi}^{1-\phi+\phi^{2}}x\cdot\left(\frac{1}{\phi^{2}}(x-(1-\phi))\right)d\mu_{sc}(x)}_{\mathsf{II}(\phi)}+\underbrace{\int_{1-\phi+\phi^{2}}^{1}xd\mu_{sc}(x)}_{\mathsf{III}(\phi)}$
$\displaystyle\geqslant\frac{-2\varepsilon}{3\pi}\left((1-(1-\phi)^{2})^{3/2}+(1-(1-\phi)^{2})^{3/2}\right)+\frac{2}{3\pi}\left(1-(1-\phi+\phi^{2})^{2}\right)^{3/2}$
$\displaystyle=\frac{2}{3\pi}\left[(1-(1-\phi+\phi^{2})^{2})^{3/2}-2\varepsilon\phi^{3/2}(2-\phi)^{3/2}\right]$

where the first inequality follows from using the closeness of the Bernstein
polynomial to the $L^{\alpha}_{\gamma}$ function, and the final inequality
follows from the observation that $\mathsf{II}(\phi)\geqslant 0$ and the use
of Proposition 3.41 to evaluate $\mathsf{I}(\varepsilon,\phi)$ and
$\mathsf{III}(\phi)$.

##### The final lower bound

Combining the upper bound on
$\operatorname*{\mathbb{E}}_{x\sim\mu_{sc}}\left[B^{L^{\alpha}_{\gamma}}(x)\right]$
via Lemma B.1 with the lower bound on the correlation term computed above
yields,

$\displaystyle\frac{\operatorname*{\mathbb{E}}_{x\sim\mu_{sc}}\left[x\cdot
B^{L^{\alpha}_{\gamma}}(x)\right]}{\operatorname*{\mathbb{E}}_{x\sim\mu_{sc}}\left[B^{L^{\alpha}_{\gamma}}(x)\right]}$
$\displaystyle\geqslant\frac{\frac{2}{3\pi}\left[(1-(1-\phi+\phi^{2})^{2})^{3/2}-2\varepsilon\phi^{3/2}(2-\phi)^{3/2}\right]}{1+\varepsilon-\mathsf{CDF}(1-\phi)}$
$\displaystyle=\frac{\frac{2}{3\pi}\left[(2\phi+2\phi^{3}-3\phi^{2}-\phi^{4})^{3/2}-2\varepsilon\phi^{3/2}(2-\phi)^{3/2}\right]}{\frac{4\sqrt{2}}{3\pi}\phi^{3/2}+o(\phi^{3/2})}$
$\displaystyle\geqslant\frac{\frac{2}{3\pi}2^{3/2}\phi^{3/2}\left((1-4\phi)^{3/2}-\frac{2}{2^{3/2}}\varepsilon(2-\phi)^{3/2}\right)}{\frac{4\sqrt{2}}{3\pi}\phi^{3/2}+o(\phi^{3/2})}$
$\displaystyle\geqslant\frac{\frac{4\sqrt{2}}{3\pi}\phi^{3/2}\left(\left(1-\frac{3}{2}4\phi+O(\phi^{2})\right)-2\varepsilon\right)}{\frac{4\sqrt{2}}{3\pi}\phi^{3/2}+o(\phi^{3/2})}$
$\displaystyle=\frac{1-20\phi-2\varepsilon}{1+o(\phi)}\,.$

For completeness, we also demonstrate an upper bound of (roughly) the same
quantitative strength on the normalized correlation between $x$ and
$B^{L^{\alpha}_{\gamma}}(x)$ under $\mu_{sc}$. To do this, we utilize a lower
bound on the normalization term given in Lemma B.1 and compute an upper bound
on $x\cdot B^{L^{\alpha}_{\gamma}}(x)$ using the uniform upper bound of
$\phi^{2}$ on the ramp part.

##### The final upper bound

A similar computation as above yields,

$\displaystyle\frac{\operatorname*{\mathbb{E}}_{x\sim\mu_{sc}}\left[x\cdot
B^{L^{\alpha}_{\gamma}}(x)\right]}{\operatorname*{\mathbb{E}}_{x\sim\mu_{sc}}\left[B^{L^{\alpha}_{\gamma}}(x)\right]}\leqslant\frac{\int_{-1}^{1}x\cdot\left(L^{\alpha}_{\gamma}(x)+\varepsilon\right)d\mu_{sc}(x)}{(1+2\varepsilon)\frac{4\sqrt{2}}{3\pi}\phi^{3/2}+O(\phi^{5/2})+o(\phi^{5/2})+\varepsilon}$
$\displaystyle=\frac{\int_{-1}^{1-\phi}x\cdot
0d\mu_{sc}(x)+\int_{1-\phi}^{1-\phi+\phi^{2}}\frac{x}{\phi^{2}}\cdot\left(x-(1-\phi)\right)d\mu_{sc}(x)+\int_{1-\phi+\phi^{2}}^{1}xd\mu_{sc}(x)+\varepsilon\operatorname*{\mathbb{E}}_{x\sim\mu_{sc}}\left[x\right]}{(1+2\varepsilon)\frac{4\sqrt{2}}{3\pi}\phi^{3/2}+O(\phi^{5/2})+o(\phi^{5/2})+\varepsilon}$
$\displaystyle\leqslant\frac{\int_{1-\phi}^{1-\phi+\phi^{2}}xd\mu_{sc}(x)+\int_{1-\phi+\phi^{2}}^{1}xd\mu_{sc}(x)}{(1+2\varepsilon)\frac{4\sqrt{2}}{3\pi}\phi^{3/2}+O(\phi^{5/2})+o(\phi^{5/2})+\varepsilon}=\frac{\int_{1-\phi}^{1}xd\mu_{sc}(x)}{(1+2\varepsilon)\frac{4\sqrt{2}}{3\pi}\phi^{3/2}+O(\phi^{5/2})+o(\phi^{5/2})+\varepsilon}$
$\displaystyle=\frac{\frac{2}{3\pi}\left(1-(1-\phi)^{2}\right)^{3/2}}{(1+2\varepsilon)\frac{4\sqrt{2}}{3\pi}\phi^{3/2}+O(\phi^{5/2})+o(\phi^{5/2})+\varepsilon}=\frac{\frac{2}{3\pi}\left(\phi(2-\phi)\right)^{3/2}}{(1+2\varepsilon)\frac{4\sqrt{2}}{3\pi}\phi^{3/2}+O(\phi^{5/2})+o(\phi^{5/2})+\varepsilon}$
$\displaystyle\leqslant\frac{\frac{2}{3\pi}2^{3/2}\phi^{3/2}}{(1+2\varepsilon)\frac{4\sqrt{2}}{3\pi}\phi^{3/2}+O(\phi^{5/2})+o(\phi^{5/2})+\varepsilon}=\frac{1}{(1+3\varepsilon)+O(\phi)+o(\phi)}\,.$

A choice of $A(\varepsilon)=O(\varepsilon)$ suffices for the desired
approximation ratios, yielding
$\phi=\frac{O(\varepsilon)}{2\nu^{\prime\prime}(\left\lVert\sigma\right\rVert^{2}_{2})^{1/2}}$.
∎

##  7 Feasibility

###  7.1 Feasibility via a near-optimal randomized algorithm

To show feasibility we will work with Algorithm 2, which is a variant of
Subag’s algorithm (Algorithm 1). Fix $0<\varepsilon<1$ and $A(\varepsilon)$.

$v_{1}\sim\mathcal{S}^{A(\varepsilon)n-1}(1)$.

$w_{1}=\frac{1}{\sqrt{k}}v_{1}$.

while $i\in\\{2,\dots,k\\}$ do

$v_{i}\sim\mathcal{S}^{\mathsf{eigen_{A(\varepsilon)n-1}}(\frac{1}{n}\nabla^{2}H(w_{i-1}))}(1)$

$w_{i}=w_{i-1}+\frac{v_{i}}{\sqrt{k}}$.

$i\leftarrow i+1$.

end while

Output $v_{k}$.

Algorithm 2 Randomized Hessian-Ascent [huang2022computational]

The only difference between Algorithm 1 and Algorithm 2 is that $x$ is chosen
uniformly at random from the unit sphere on the $\delta n$-dimensional
subspace corresponding to the top eigenvectors of the Hessian projected
orthogonal to $v_{i-1}$ as $\nabla^{2}_{\perp}H(v_{i-1})$. Critically,
Algorithm 2 continues to be nearly optimal for the same reason as the original
algorithm is [subag2021following, Lemma 3 & Theorem 4] in conjunction with the
fact that the first step can be arbitrary. A more detailed explanation for the
equivalence between these two algorithms can be found in
[huang2022computational, Section 3.2.7].

Consequently, showing feasibility amounts to showing that Algorithm 2
satisfies, at every iteration $i\in[k]$, the constraints of the SoS program at
the corresponding cumulant SoS variable. Demonstrating this is done by showing
that:

  1. 1.

The uniform distribution on the unit-sphere satisfies the conditional maximum-
entropy constraints imposed on valid pseudo-distributions at every iterate
$i\in[k]$.

  2. 2.

The cumulative norm increments have the same norms _and_ are incrementally
orthogonal. This is trivially true by construction and is omitted in the
formal proof.

  3. 3.

The conditional pseudo-expectation of the quadratic form of the hessian under
$v_{i}$ is at least $(1-\varepsilon)2\nu^{\prime\prime}(q_{i})$.

The statements above are proved formally in Lemma 7.1 below.

######  Lemma 7.1 (Feasibility of the SoS program).

Given $0<\varepsilon<1$ and $A(\varepsilon)$, the following facts hold:

  1. 1.

For $i=1,\dots,k$,

$\displaystyle\left\\{\widetilde{\kappa}_{{j}}\left[{v_{i}\,\mid\,v_{i-1}}\right]=o_{n}(1)\right\\}_{1,\dots,i-2}\,,j\neq
2\,,$
$\displaystyle\left\\{\left\lVert\operatorname*{\mathbb{E}}[v_{i}v_{i}^{\mathsf{T}}\,\mid\,v_{i-1}]\right\rVert_{{}_{\mathsf{op}}}\leqslant\frac{1+o_{n}(1)}{A(\varepsilon)n}\right\\}_{\mid
1,\dots,i-2}\,.$

  2. 2.

For $i=1,\dots,k$,

$\displaystyle\operatorname*{\mathbb{E}}_{x\sim\mathcal{S}^{A(\varepsilon)n-1}(\mathsf{eigen}\frac{1}{n}(\nabla^{2}_{\perp}H(w_{i-1})))}[x^{\mathsf{T}}\frac{1}{n}\nabla^{2}_{\perp}H(w_{i-1})x]\geqslant
2\nu^{\prime\prime}(q)^{1/2}-\varepsilon\,.$

###### Proof.

Note that every probability distribution is a pseudo-distribution, since
$\operatorname*{\mathbb{E}}_{x}[1]=1$ (by normalization of probability) and
$\operatorname*{\mathbb{E}}_{x}[p^{2}(x)]\geqslant 0$ for any polynomial $p$
since it can be viewed as an inner product between non-negative quantities.
Therefore, we set
$\operatorname*{\widetilde{\mathbb{E}}}=\operatorname*{\mathbb{E}}_{x}$ to
show that the maximum-entropy constraints bounds are satisfied.
[subag2021following, Lemma 3] immediately implies that the expected quadratic
form of the hessian is arbitrarily close to the largest eigenvalue under the
uniform distribution on the sphere of the top $A(\varepsilon)n$-dimensional
subspace of the Hessian orthogonal to the current iterate.

##### The $\left\lVert\cdot\right\rVert_{{}_{\mathsf{op}}}$ constraint

Proving this constraint amounts to showing that the operator norm of the
covariance of $v_{i}$ conditioned on $v_{1},\dots,v_{i-1}$ is bounded from
above by $\frac{1}{A(\varepsilon)n}$. Since
$v_{i}\sim\mathcal{S}^{A(\varepsilon)n-1}(\mathsf{eigen}(\frac{1}{n}\nabla^{2}_{\perp}H(w_{i-1})))$,
note that conditioned on $v_{1},\dots,v_{i-2}$,

$\displaystyle\left\lVert\operatorname*{\mathbb{E}}[v_{i}v_{i}^{\mathsf{T}}\,\mid\,v_{i-1}]\right\rVert_{{}_{\mathsf{op}}}$
$\displaystyle=\left\lVert\operatorname*{\mathbb{E}}_{v\sim\mathcal{S}^{A(\varepsilon)n-1}(\mathsf{eigen}(\frac{1}{n}\nabla^{2}_{\perp}H(w_{i-1})))}\left[vv^{\mathsf{T}}\right]\right\rVert_{{}_{\mathsf{op}}}$
$\displaystyle=_{\text{Rotational Invariance of Haar-
Measure}}\left\lVert\operatorname*{\mathbb{E}}_{v^{\prime}\sim\mathcal{S}^{n-1}(1^{A(\varepsilon)n}\oplus\,0^{(1-A(\varepsilon))n})}\left[v^{\prime}v^{\prime\mathsf{T}}\right]\right\rVert_{{}_{\mathsf{op}}}$
$\displaystyle=_{\text{\lx@cref{creftypecap~refnum}{lem:denom-conc} with
}q=1}\left\lVert\frac{1-o_{n}(1)}{A(\varepsilon)n}\operatorname*{\mathbb{E}}_{x\sim\mathcal{N}(0,\operatorname{\mathrm{Id}}_{A(\varepsilon)n}\oplus
0_{(1-A(\varepsilon))n})}\left[xx^{\mathsf{T}}\right]\right\rVert_{{}_{\mathsf{op}}}$
$\displaystyle=\frac{1-o_{n}(1)}{A(\varepsilon)n}\left\lVert\begin{bmatrix}1&\\\
&\ddots&\\\ &&1\\\ &&&0&\\\ &&&&\ddots&\\\
&&&&&0\end{bmatrix}\right\rVert_{{}_{\mathsf{op}}}\leqslant\frac{1+o_{n}(1)}{A(\varepsilon)n}\,,$

where in the above calculation
$v^{\prime}\sim\mathcal{S}^{n-1}(1^{A(\varepsilon)n}\oplus
0^{(1-A(\varepsilon))n})$ refers to a unit vector chosen uniformly at random
from the $A(\varepsilon)n$-dimensional sphere where the last
$(1-A(\varepsilon))n$-coordinates are $0$.

##### The $\kappa_{j}[v_{i}\mid v_{i-1}]$ constraint

To show this constraint is satisfied, we follow a similar strategy. Note that,
$\forall j\neq 2$, the following holds for every entry of the $j$-th
conditional cumulant tensor,

$\displaystyle\kappa_{j}[v_{i}\mid v_{i-1}]$ $\displaystyle=_{\text{by moment-
cumulant
formula}}\sum_{\pi\in\mathcal{P}(k_{1},\dots,k_{j})}(-1)^{|\pi|-1}(|\pi|-1)!\prod_{b\in\pi}\operatorname*{\mathbb{E}}\left[\prod_{a\in
b}v_{i_{a}}\right]$ $\displaystyle=_{\text{by
\lx@cref{creftypecap~refnum}{lem:denom-conc} and rotational
invariance}}\sum_{\pi\in\mathcal{P}(k_{1},\dots,k_{j})}(-1)^{|\pi|-1}(|\pi|-1)!\prod_{b\in\pi}\frac{1-o_{n}(1)}{\left(A(\varepsilon)n\right)^{|b|}}\operatorname*{\mathbb{E}}\left[\prod_{a\in
b}v^{\prime}_{i_{a}}\right]$
$\displaystyle=\sum_{\pi\in\mathcal{P}(k_{1},\dots,k_{j})}\frac{(-1)^{|\pi|-1}(|\pi|-1)!(1-o_{n}(1))}{(A(\varepsilon)n)^{j}}\prod_{b\in\pi}\operatorname*{\mathbb{E}}\left[\prod_{a\in
b}v^{\prime}_{i_{a}}\right]\,,$

where in the above calculation we again used 3.25. Now, since this is a sum of
the form,

$\sum_{\pi}\frac{(1-o_{n}(1))}{(\delta n)^{j}}a_{\pi}\,,$

we can conclude that,

$\left\lVert\kappa_{j}[v_{i}\mid
v_{i-1}]-\kappa_{j}\left[\mathcal{N}(0,\operatorname{\mathrm{Id}}_{A(\varepsilon)n\oplus
0_{(1-A(\varepsilon))n}})\right]\right\rVert_{F}^{2}\leqslant o_{n}(1)\,.$

##### The approximate-LDP constraint

We now show that the original algorithm’s updates satisfy the approximate-LDP
constraint in the program by demonstrating that the choice of polynomial
corresponding to a Bernstein polynomial of sufficiently large (but constant)
degree has the following properties:

  1. 1.

It permits a (weakly) orthogonal decomposition using the Hermite-Weiner
fourier basis, by virtue of being square integrable on $\mathcal{S}^{n-1}(1)$.

  2. 2.

By Proposition A.2, it is arbitrarily close in Frobenius norm to the projector
into the top $A(\varepsilon)n$-dimensional eigenspace of
$\nabla^{2}H(\sigma)$.

Note that any $f\in\ell_{2}(\Sigma)$ can be expressed as,

$f(x)=\sum_{\alpha\in\mathbb{N}^{n}}\left\langle
f,\mathsf{He}_{\alpha}(x\mid\Sigma)\right\rangle_{\ell^{2}(\Sigma)}\mathsf{He}_{\alpha}(x\mid\Sigma)\,.$

Note that for $x=\sum_{j=1}^{i-1}v_{j}$ with $\left\lVert
v_{j}\right\rVert^{2}=1$ and
$f=p_{d}:=B^{d}_{L^{\alpha}_{\gamma}}\left(c_{i}\operatorname{\mathrm{Id}}_{n}+\frac{1}{n}\nabla^{2}H(\sigma)\right)$
with the Fourier decomposition taken element-wise in the matrix polynomial, it
is the case that,

$\displaystyle\left\lVert f_{ij}\right\rVert^{2}_{2}:=\langle
f_{ij}(x),f_{ij}(x)\rangle_{\ell^{2}(\Sigma)}=\int_{x\sim\mathcal{N}(0,\Sigma)}\frac{1}{\delta
n}B^{d}_{L^{\alpha}_{\gamma}}\left(c_{i}(\operatorname{\mathrm{Id}}_{n})_{ij}+\frac{1}{n}\nabla^{2}H(x)_{ij}\right)dg(x)<\infty\,.$

Since the algorithm samples uniformly at random from the
$A(\varepsilon)n$-dimensional sphere corresponding to the top eigenspace of
$\nabla^{2}H(\sigma)$, we can then apply Lemma 7.2 to conclude that the
appropriately scaled conditional covariance of this algorithm is basically a
projector into this subspace times the identity matrix, and a triangle
inequality then immediately implies feasibility as demonstrated below:

$\displaystyle\left\lVert\operatorname*{\mathbb{E}}\left[v_{i}v_{i}^{\mathsf{T}}\mid
w_{i-1}\right]-p_{d}(w_{i-1})\right\rVert_{F}^{2}=\left\lVert\operatorname*{\mathbb{E}}\left[v_{i}v_{i}^{\mathsf{T}}\mid
w_{i-1}\right]-\sum_{\alpha,|\alpha|\leqslant
d}\langle\operatorname*{\mathbb{E}}\left[v_{i}v_{i}^{\mathsf{T}}\mid
w_{i-1}\right],\mathsf{He}_{d}\rangle_{\ell^{2}(\Pi_{\delta
n})}\right\rVert_{F}^{2}$
$\displaystyle=\left\lVert\operatorname*{\mathbb{E}}\left[v_{i}v_{i}^{\mathsf{T}}\mid
w_{i-1}\right]-\frac{1}{\delta
n}B^{d}(c_{i}\operatorname{\mathrm{Id}}_{n}+\frac{1}{n}\nabla^{2}H(\sigma))\right\rVert_{F}^{2}$
$\displaystyle\leqslant_{\text{triangle
inequality}}\left\lVert\operatorname*{\mathbb{E}}\left[v_{i}v_{i}^{\mathsf{T}}\mid
w_{i-1}\right]-\frac{1}{\delta n}\Pi_{\delta
n}\right\rVert^{2}_{F}+\frac{1}{\delta n}\left\lVert\Pi_{\delta
n}-\frac{1}{\delta
n}B^{d}(c_{i}\operatorname{\mathrm{Id}}_{n}+\frac{1}{n}\nabla^{2}H(\sigma))\right\rVert^{2}_{F}$
$\displaystyle=_{\text{by \lx@cref{creftypecap~refnum}{prop:matrix-bernstein-
linear-
ramp}}}\left\lVert\operatorname*{\mathbb{E}}\left[v_{i}v_{i}^{\mathsf{T}}\mid
w_{i-1}\right]-\Pi_{\delta n}\right\rVert^{2}_{F}+o_{n}(1)$
$\displaystyle=_{\text{by \lx@cref{creftypecap~refnum}{lem:denom-
conc}}}\left\lVert\frac{1-o_{n}(1)}{\delta
n}\operatorname*{\mathbb{E}}_{x\sim\mathcal{N}\left(0,\,\Pi_{\mathsf{eigen}\left(\nabla^{2}_{\perp}H(w_{i-1})\right)}\right)}[(xx^{\mathsf{T}})]-\frac{1}{\delta
n}\Pi_{\delta n}\right\rVert_{F}^{2}+o_{n}(1)$
$\displaystyle=_{\operatorname*{\mathbb{E}}[xx^{\mathsf{T}}]=\Pi_{\delta
n}}\left\lVert\frac{1-o(n(1))}{\delta n}\Pi_{\delta n}-\frac{1}{\delta
n}\Pi_{\delta n}\right\rVert_{F}^{2}+o_{n}(1)$
$\displaystyle=\left\lVert-\frac{o_{n}(1)}{\delta n}\Pi_{\delta
n}\right\rVert_{F}^{2}+o_{n}(1)=o_{n}(1)+o_{n}(1)=o_{n}(1)\,,$

where we used the fact that the covariance of the Gaussian is given by an
orthogonal transformation into the top $\delta n$-dimensional subspace of
$\nabla^{2}H(w_{i-1})$, which in this case is just $\Pi_{\delta
n}=\Pi_{\mathsf{eigen}\left(\nabla^{2}_{\perp}H(w_{i-1})\right)}$. For every
$q\in\mathbb{N}$,

$\displaystyle\operatorname*{\mathbb{E}}[(v_{i}v_{i}^{\mathsf{T}})^{\otimes
q}\mid w_{i-1}]$
$\displaystyle=_{\text{\lx@cref{creftypecap~refnum}{lem:denom-
conc}}}\frac{1-o_{n}(1)}{\left(\delta
n\right)^{q}}\operatorname*{\mathbb{E}}_{x\sim\mathcal{N}(0,\operatorname{\mathrm{Id}}_{A(\varepsilon)n}\oplus
0_{(1-A(\varepsilon))n})}\left[(xx^{\mathsf{T}})^{\otimes q}\right]$
$\displaystyle=_{\text{\lx@cref{creftypecap~refnum}{lem:isserlis}}}(q-1)!!\,\mathrm{sym}\left[\left(\operatorname*{\mathbb{E}}[xx^{\mathsf{T}}]\right)^{\otimes
q}\right]$
$\displaystyle=(q-1)!!\,\mathrm{sym}\left[\left((1-o_{n}(1))\operatorname*{\mathbb{E}}[v_{i}v_{i}^{\mathsf{T}}\mid
w_{i-1}]\right)^{\otimes q}\right]$
$\displaystyle=(1-o_{n}(1))(q-1)!!\,\mathrm{sym}\left[\left(\operatorname*{\mathbb{E}}[v_{i}v_{i}^{\mathsf{T}}\mid
w_{i-1}]\right)^{\otimes q}\right]\,.$

Rewrite the conditional covariance as
$\operatorname*{\mathbb{E}}[v_{i}v_{i}^{\mathsf{T}}\mid w_{i-1}]=p(w_{i-1})+B$
where $B$ is a matrix with $o_{n}(1)$ Frobenius norm (as shown above). This
yields that,

$\operatorname*{\mathbb{E}}[(v_{i}v_{i}^{\mathsf{T}})^{\otimes q}\mid
w_{i-1}]=(1-o_{n}(1))(q-1)!!\,\mathrm{sym}\left[\left(p(w_{i-1})+B\right)^{\otimes
q}\right]\,,$

and, therefore,

$\displaystyle\left\lVert\operatorname*{\mathbb{E}}[(v_{i}v_{i}^{\mathsf{T}})^{\otimes
q}\mid w_{i-1}]-p(w_{i-1})^{\otimes q}\right\rVert_{F}^{2}$
$\displaystyle\leqslant(q-1)!!\left\lVert\sum_{i\in[q],j=q-i}\mathrm{sym}\left(B^{\otimes
i}\otimes p(w_{i-1})^{\otimes j}\right)\right\rVert_{F}^{2}$
$\displaystyle\leqslant_{\text{$\mathrm{sym}$ is
contractive}}(q-1)!!\left\lVert\sum_{i\in[q],j=q-i}B^{\otimes i}\otimes
p(w_{i-1})^{\otimes j}\right\rVert_{F}^{2}$
$\displaystyle=_{\text{\lx@cref{creftypecap~refnum}{lem:l1-to-
lt}}}q^{2}(q-1)!!\sum_{i\in[q],j=q-i}\left\lVert B^{\otimes i}\otimes
p(w_{i-1})^{\otimes j}\right\rVert_{F}^{2}$
$\displaystyle\leqslant_{\text{\lx@cref{creftypecap~refnum}{lem:scalar-am-gm}
+ \lx@cref{creftypecap~refnum}{prop:matrix-bernstein-linear-
ramp}}}\frac{q^{3}(q-1)!!}{\delta^{q}}\frac{1}{n^{q}}\,.$

##### Soundness & Optimality of Taylor Expansion

The last step is to invoke soundness for all distributions (Section 4.1.10)
which implies that the corresponding pseudo-distribution comprising the low-
degree moments of the “Subag distribution” implicit in Algorithm 2 satisfies
the entropy constraints. Invoking [arous2020geometry, Corollary 59] in
conjunction with a Taylor series argument immediately implies that the hessian
has the largest contribution, and the final energy is lower-bounded by
$\sum_{i=1}^{k}\frac{1}{k}\nu^{\prime\prime}(q_{i})^{1/2}-\varepsilon\cdot
k(\varepsilon)-O\left(\frac{1}{k(\varepsilon)^{1/2}}\right)$, which is
arbitrarily close to $\int_{0}^{1}\nu^{\prime\prime}(q)^{1/2}dq$ with the
appropriate choice of $k(\varepsilon)$. The upper bound follows from Lemma
3.7, as no configuration can have value greater than the true maximum. ∎

###  7.2 Equivalence of rescaled spherical and Gaussian moments

We state and prove a small lemma that allows us to freely go between the low-
degree moments of the uniform measure on the real unit sphere and an
appropriate rescaling of the comparable Gaussian moments.

######  Lemma 7.2 (Equivalence of re-scaled gaussian and uniform moments).

Given $x\sim\mathcal{N}(0,\operatorname{\mathrm{Id}}_{n})$ and
$v\sim\mathsf{Unif}(\mathcal{S}^{n-1}(1))$, the following holds,

$\operatorname*{\mathbb{E}}_{v}[(vv^{\mathsf{T}})^{\otimes
q}]=\frac{(1-O_{n}(1/n))}{n^{q}}\operatorname*{\mathbb{E}}_{x}[(xx^{\mathsf{T}})^{\otimes
q}]\,.$ (7.1)

###### Proof.

The proof will use the rotational invariance of both distributions and an
explicit entry-by-entry comparison between the rescaled Gaussian moment matrix
and the real Haar distribution moment matrix.

$\displaystyle\operatorname*{\mathbb{E}}_{x}[(xx^{\mathsf{T}})^{\otimes
q}]=\sum_{a_{1},\dots,a_{2q}=1}^{n}\operatorname*{\mathbb{E}}_{x}[x_{a_{1}}\cdots
x_{a_{2q}}](e_{a_{1}}\otimes\cdots\otimes
e_{a_{q}})(e_{a_{q+1}}\otimes\cdots\otimes e_{a_{2q}})^{\mathsf{T}}\,.$

By exactly the same argument,

$\displaystyle\operatorname*{\mathbb{E}}_{v}[(vv^{\mathsf{T}})^{\otimes
q}]=\sum_{a_{1},\dots,a_{2q}=1}^{n}\operatorname*{\mathbb{E}}_{v}[v_{a_{1}}\cdots
v_{a_{2q}}](e_{a_{1}}\otimes\cdots\otimes
e_{a_{q}})(e_{a_{q+1}}\otimes\cdots\otimes e_{a_{2q}})^{\mathsf{T}}\,.$

It is elementary to see from Lemma 3.10 and the definition of the
$\Gamma$-function that the $2m_{i}$-th moment of $x_{i}\sim\mathcal{N}(0,1)$
is,

$\operatorname*{\mathbb{E}}_{x_{i}\sim\mathcal{N}(0,1)}[x^{2m_{i}}_{i}]=(2m_{i}-1)(2m_{i}-3)\cdots(3)(1)=\frac{(2m_{i})!}{2^{m_{i}}m_{i}!}=\frac{\Gamma(m_{i}+1/2)2^{m_{i}}}{\sqrt{\pi}}\,.$
(7.2)

Symmetry implies that the odd moments of the $1$-marginals of the uniform
distribution are $0$. By [folland2001integrate, Theorem, Pg. 2], the even
moments are

$\displaystyle\operatorname*{\mathbb{E}}_{v}[v_{a_{1}}\cdots v_{a_{2q}}]$
$\displaystyle=\operatorname*{\mathbb{E}}_{v}\left[\prod_{i=1}^{\ell}(v_{i}^{2})^{m_{i}}\right]$
$\displaystyle=\int_{v\in\mathbb{R}^{n}}\left(\prod_{i=1}^{\ell}(v_{i}^{2})^{m_{i}}\right)d\sigma=\frac{\int_{v\in\mathcal{S}^{n-1}(1)}v^{\alpha}d\sigma}{\int_{v\in\mathcal{S}^{n-1}(1)}d\sigma}$
$\displaystyle=\frac{2\Gamma(m_{1}+1/2)\Gamma(m_{2}+1/2)\cdots\Gamma(m_{\ell}+1/2)\left(\Gamma(1/2)\right)^{n-\ell}}{\Gamma(q+n/2)}\frac{\Gamma(n/2)}{2\Gamma(1/2)^{n}}$
$\displaystyle=\frac{\operatorname*{\mathbb{E}}_{x}[x_{a_{1}}\cdots
x_{a_{2q}}]}{2^{q}}\frac{\Gamma(n/2)}{\Gamma(q+n/2)}$
$\displaystyle=\frac{\operatorname*{\mathbb{E}}_{x}[x_{a_{1}}\cdots
x_{a_{2q}}]}{n(n+2)(n+4)\cdots(n+2q-2)}\,,$ (7.3)

where $2(m_{1}+\cdots+m_{\ell})=2q$ with $m_{i}$ denoting the number of times
an index $a_{i}\in[n]$ appears in $\\{a_{1},\dots,a_{2q}\\}$ and $\ell$ being
the number of distinct indices. Combining Section 7.2 and Equation 7.2
immediately yields that,

$\operatorname*{\mathbb{E}}_{v}[v_{a_{1}}\cdots
v_{a_{2q}}]=\frac{(1-O_{n}(1/n))}{n^{q}}\operatorname*{\mathbb{E}}_{x}[x_{a_{1}}\cdots
x_{a_{2}q}]\,,$ (7.4)

for every $a_{1},\dots,a_{2q}\in\\{1,\dots,n\\}$. ∎

##  8 Local iterative algorithms vs. HES SoS hierarchy

In this section, we provide three concrete examples of inputs from a
“modified” spherical spin glass ensemble, and then give proofs for the first
two (and a sketch for for the latter) of the sub-optimality of local iterative
algorithms on these (such as Hessian ascent [subag2021following] and its
randomized generalization (Algorithm 2)) along with arguments for why the HES
SoS relaxation provides certificates with better values. The last example, in
particular, motivates the conjecture that the HES SoS relaxation is capable of
“disentangling” arbitrary correlations in the input, provided they are
sufficiently smooth.

###  8.1 Invariance with respect to the spherical ideal

The first example is one that “fools” the randomized version (Algorithm 2) of
the Hessian ascent [subag2021following] algorithm by a scaling of the degree-2
component.

######  Example 8.1 (Homogeneous degree scaling).

Let the Hamiltonian $H$ be,

$H(\sigma)=\alpha_{2}\left\lVert\sigma\right\rVert^{100}_{2}H_{2}(\sigma)+H_{4}(\sigma)\,,$

for some choice of $\alpha>0$. Then, the algorithm given in Algorithm 2 is
provably suboptimal, but the HES SoS relaxation does better.

###### Proof.

A proof sketch for the claim above is given, whereby the geometry and scaling
of the top $\delta n$-dimensional eigenspace of the Hessian is calculated, and
then it is demonstrated that Algorithm 2 initially follows the _incorrect_
direction, thereby accumulating an incorrect amount of energy before it
changes direction. However, the HES relaxation is able to use the Schatten
$p$-norm of the collective Hessian to “distinguish” between the space where
the large (and dominant) amount of signal from the top part of the spectrum
resides, while the rank-$1$ terms contribute signal that is suppressed by the
moment scalings implied by the HES constraints.

##### Eigenspectrum of the Hessian

When $\sigma\in\mathcal{S}^{n-1}(1)$, note that $H(\sigma)$ is
indistinguishable from an instance of a spherical spin glass Hamiltonian. The
Hessian of $H(\sigma)$ is,

$\displaystyle\frac{1}{n}\nabla^{2}H(\sigma)=$
$\displaystyle\,\,\alpha_{2}\left\lVert\sigma\right\rVert^{100}_{2}\frac{1}{n}\nabla^{2}H_{2}(\sigma)+\frac{1}{n}\nabla^{2}H_{4}(\sigma)+100\alpha_{2}\left\lVert\sigma\right\rVert_{2}^{98}\frac{1}{n}H_{2}(\sigma)\operatorname{\mathrm{Id}}_{n}$
$\displaystyle+\underbrace{9800\alpha\left\lVert\sigma\right\rVert_{2}^{96}\frac{1}{n}H_{2}(\sigma)\sigma\sigma^{\mathsf{T}}+100\alpha_{2}\left\lVert\sigma\right\rVert_{2}^{98}\left(\sigma\left(\nabla\frac{1}{n}H_{2}(\sigma)\right)^{\mathsf{T}}+\frac{1}{n}\nabla
H_{2}(\sigma)\sigma^{\mathsf{T}}\right)}_{\text{rank-1 terms}}\,.$

By [subag2021following, Proof of Lemma-3], it is clear that,

$\alpha\left\lVert\sigma\right\rVert_{2}^{100}\frac{1}{n}\nabla^{2}H_{2}(\sigma)\sim\sqrt{2\alpha_{2}}\left\lVert\sigma\right\rVert_{2}^{50}\mathsf{GOE}(n)\,,\text{
and},\,\frac{1}{n}\nabla^{2}H_{4}(\sigma)\sim\sqrt{12}\mathsf{GOE}(n)\,.$

Furthermore,

$\displaystyle
H_{2}(\sigma)\leqslant\sup_{\sigma\in\mathcal{S}^{n-1}(\left\lVert\sigma\right\rVert^{2}_{2})}H_{2}(\sigma)=_{\text{\lx@cref{creftypecap~refnum}{prop:frsb-
equivalent}}}\sqrt{2}\int_{0}^{\left\lVert\sigma\right\rVert^{2}_{2}}dx=\sqrt{2}\left\lVert\sigma\right\rVert^{2}_{2}\,.$

This immediately implies that _all_ the (redundant) eigenvalues of
$100\alpha_{2}\left\lVert\sigma\right\rVert^{98}_{2}\frac{1}{n}H_{2}(\sigma)\operatorname{\mathrm{Id}}_{n}$
are no more than
$100\sqrt{2}\alpha_{2}\left\lVert\sigma\right\rVert^{100}_{2}$.

##### Suboptimality of the rank-$1$ components

Suboptimality of the rank-$1$ terms is a consequence of the analysis of Subag
[subag2018free, subag2021following] but can also be rigorously justified by
noticing that they have large alignment with the gradient which leads to
larger steps (under the HES constraints) and sub-optimal gain of energy after
a certain number of steps777This can be made rigorous via an analysis based on
a Taylor expansion argument which gives a precise iterate $i$ after which it
becomes suboptimal to follow the gradient term..

##### Threshold for suboptimality of local algorithm

We show the sub-optimality of the variant of Subag’s algorithm stated in
Algorithm 2. Note that, since $v\sim\mathcal{S}^{\mathsf{eigen}_{\delta
n}(\nabla^{2}H(\sigma))}(1)_{\perp\sigma}$, by construction,

$\displaystyle\left\langle
9800\alpha_{2}\left\lVert\sigma\right\rVert_{2}^{96}\frac{1}{n}H_{2}(\sigma)\sigma\sigma^{\mathsf{T}}+100\alpha_{2}\left\lVert\sigma\right\rVert_{2}^{98}\left(\sigma\left(\nabla\frac{1}{n}H_{2}(\sigma)\right)^{\mathsf{T}}+\frac{1}{n}\nabla
H_{2}(\sigma)\sigma^{\mathsf{T}}\right),vv^{\mathsf{T}}\right\rangle$
$\displaystyle=0+\left\langle
9800\alpha_{2}\left\lVert\sigma\right\rVert_{2}^{2}\sigma\left(\frac{1}{n}H_{2}(\sigma)\right)^{\mathsf{T}},vv^{\mathsf{T}}\right\rangle+0=0+0+9800\alpha_{2}\left\lVert\sigma\right\rVert^{2}_{2}\left\langle\frac{1}{n}H(\sigma),v\right\rangle\langle\sigma,v\rangle=0\,.$

Consequently, the algorithm chooses the next step uniformly at random from the
unit sphere in the space spanned by the top $\delta n$-dimensional eigenspace
of
$\alpha\left\lVert\sigma\right\rVert^{2}_{2}\frac{1}{n}\nabla^{2}_{\perp\sigma}H_{2}(\sigma)+\frac{1}{n}\nabla^{2}_{\perp\sigma}H_{4}(\sigma)+100\alpha\left\lVert\sigma\right\rVert^{98}_{2}\frac{1}{n}H_{2}(\sigma)\left(\operatorname{\mathrm{Id}}_{n}-\sigma\sigma^{\mathsf{T}}\right)$,
where the Hessians have been projected orthogonal to $\sigma$. Both the
Hessians are symmetric with independent Gaussian coefficients, and therefore,
by [subag2021following, Proof of Lemma 3],

$\displaystyle\alpha_{2}\left\lVert\sigma\right\rVert^{2}_{100}\frac{1}{n}\nabla^{2}_{\perp\sigma}H_{2}(\sigma)+\frac{1}{n}\nabla^{2}_{\perp\sigma}H_{4}(\sigma)+100\alpha\left\lVert\sigma\right\rVert^{98}_{2}\frac{1}{n}H_{2}(\sigma)\left(\operatorname{\mathrm{Id}}_{n}-\sigma\sigma^{\mathsf{T}}\right)$
$\displaystyle\sim\left(\sqrt{2\alpha_{2}}\left\lVert\sigma\right\rVert_{2}^{50}\right)\mathsf{GOE}(n-1)\boxplus\sqrt{12}\mathsf{GOE}(n-1)+100\alpha\left\lVert\sigma\right\rVert_{2}^{98}\frac{1}{n}H_{2}(\sigma)$
$\displaystyle\sim\left(\sqrt{2\alpha_{2}}\left\lVert\sigma\right\rVert_{2}^{50}\right)\mathsf{GOE}(n-1)\boxplus\sqrt{12}\mathsf{GOE}(n-1)+100\alpha\left\lVert\sigma\right\rVert_{2}^{98}\frac{1}{n}H_{2}(\sigma)$
$\displaystyle\sim_{\text{free
convolution~{}\cite[cite]{[\@@bibref{}{mingo2017free}{}{}]}}}\left(\sqrt{2\alpha_{2}\left\lVert\sigma\right\rVert_{2}^{100}+12}\right)\mathsf{GOE}(n-1)+100\alpha_{2}\left\lVert\sigma\right\rVert_{2}^{98}\frac{1}{n}H_{2}(\sigma)\,.$

Now, note that, for $\left\lVert\sigma\right\rVert^{2}_{2}$ sufficiently
small, as the top $\delta n$-dimensional eigenbases of
$\frac{1}{n}\nabla^{2}_{\perp\sigma}H_{2}(\sigma)$ and
$\frac{1}{n}\nabla^{2}_{\perp\sigma}H_{4}(\sigma)$ will be mostly orthogonal
(which can be verified using the fact that symmetric $\mathsf{GOE}$ matrices
are orthogonally invariant), one can see that the local algorithm will (with
high probability) move in the direction of
$\frac{1}{n}\nabla^{2}_{\perp\sigma}H_{4}(\sigma)$. However, as
$\left\lVert\sigma\right\rVert^{2}_{2}\to 1$, the dominating signal will come
(for sufficiently large $\alpha_{2}>0$) from
$\frac{1}{n}\nabla^{2}_{\perp\sigma}H(\sigma)$.

##### Robustness of the HES relaxation

The HES SOS relaxation will use the following low-degree extension,

$\tilde{H}(\sigma)=\alpha_{2}H_{2}(\sigma)+H_{4}(\sigma)\,,$

which is simply equivalent to the original spherical spin-glass Hamiltonian.
Consequently, Theorem 1.1 implies that the SoS HES program certifies the
correct bounds and outputs a solution with the same energy. ∎

_Note_ : While the above argument demonstrates that the HES relaxation is
automatically robust to the sort of perturbed Hamiltonian shown, a small
modification to the local algorithm to quotient the input Hamiltonian by the
ideal $\left\lVert\sigma\right\rVert_{2}^{2p}=\left(\frac{i}{k}\right)^{p}$
(which is just the $p$-th power of the length of the current step
$\sigma_{i}$) fixes this example.

###  8.2 Global information: The direct-sum case

We now use projectors $\\{\Pi,\Pi^{\perp}\\}$ to project into subspaces, and
then define different degree components to act in complementary subspaces.
This leads to examples of Hamiltonians which confuse the Hessian ascent
algorithm in a way that doesn’t seem to be fixable with _local_ information,
since one would require global knowledge of a path in the _future_ (on
average) to prevent the algorithm from getting “stuck” in the wrong subspace.

######  Example 8.2 (Disjoint-subspace projector).

Let the hamiltonian $H$ be,

$H(\sigma)=\alpha_{1}H_{2}(\Pi\sigma)+\alpha_{2}H_{4}(\Pi^{\perp}\sigma)\,.$

Then, there exists a choice of $\alpha_{1},\alpha_{2}\geqslant 0$ which causes
Algorithm 2 to be provably suboptimal, but the HES SoS relaxation does better.

###### Proof.

We will normalize contributions to the energy in the both sub-spaces, so that
the “history” of energy gains from the degree-$2$ component can _never_ fool
the HES SoS hierarchy to _only_ follow the degree-$2$ steps, but the Hessian
ascent [subag2021following] algorithm (and its randomized version) will always
follow the degree-$2$ component initially, thereby making it suboptimal.

##### Eigenspectrum of the Hessian

It is clear to see that the Hessian of the as-is Hamiltonian is,

$\frac{1}{n}\nabla^{2}H(\sigma)=\frac{1}{n}\left(\alpha_{2}\Pi\nabla^{2}H_{2}(\Pi\sigma)\Pi+\alpha_{4}\Pi^{\perp}H_{4}(\Pi^{\perp}\sigma)\Pi^{\perp}\right)\,,$

and that it has a spectral law given by,

$\frac{1}{n}\nabla^{2}H(\sigma)\sim\sqrt{2\alpha_{2}}\mathsf{GOE}(\operatorname{Tr}[\Pi])\oplus\left\lVert\Pi^{\perp}\sigma\right\rVert^{2}_{2}\sqrt{12\alpha_{4}}\mathsf{GOE}(\operatorname{Tr}[\Pi^{\perp}])\,.$

##### Normalization of contributions to $\Pi$ & $\Pi^{\perp}$

We now choose $\alpha_{2},\alpha_{4}\geqslant 0$, such that, the contribution
from the degree-2 part forces the Hessian ascent algorithm into the
$\Pi$-subspace, thereby, only following the degree-2 signal. Note that, near
the origin (at the $i$-th step), we want

$\sqrt{2\alpha_{2}}>>\left\lVert\Pi^{\perp}\sigma\right\rVert^{2}_{2}\sqrt{12\alpha_{4}}\,,$

which holds for a choice of
$\alpha_{4}<<\frac{1}{\left\lVert\Pi^{\perp}\sigma\right\rVert^{2}_{2}}\frac{\alpha_{2}}{6}$.
A concrete choice of
$\alpha_{4}=\frac{1}{\left\lVert\Pi^{\perp}\sigma\right\rVert^{2}_{2}}\frac{\alpha_{2}}{12}$
suffices to obstruct the local algorithm (since it will be mislead within the
first few steps).

##### Robustness of the SoS HES relaxation

The previous choice of normalization motivates the following extended
hamiltonian in the ball,

$\tilde{H}(\sigma)=\alpha_{2}\left\lVert\sigma\right\rVert^{2}_{2}H_{2}(\Pi\sigma)+\alpha_{4}H_{4}(\Pi^{\perp}\sigma)\,.$

After some algebra, it is easy to see that the Hessian, ignoring rank-$1$
terms, is given as,

$\frac{1}{n}\nabla^{2}\tilde{H}(\sigma)=\frac{1}{n}\left(\alpha_{2}\left\lVert\sigma\right\rVert^{2}_{2}\Pi\nabla^{2}\tilde{H}_{2}(\Pi\sigma)\Pi+\alpha_{4}\Pi^{\perp}\nabla^{2}\tilde{H}_{4}(\Pi^{\perp}\sigma)\Pi^{\perp}+\alpha_{2}\tilde{H}_{2}(\Pi\sigma)\operatorname{\mathrm{Id}}_{n}\right)\,.$

Assuming that $\sigma$ is the $i$-th step, the eigenspectrum of the above
Hessian (modulo the rank-$1$ contributions) is,

$\displaystyle\left\\{\sqrt{2\alpha_{2}\left\lVert\sigma\right\rVert^{2}_{2}}\mathsf{GOE}(\operatorname{Tr}[\Pi])\oplus\sqrt{12\alpha_{4}\left\lVert\Pi^{\perp}\sigma\right\rVert^{2}_{2}}\mathsf{GOE}(\operatorname{Tr}[\Pi^{\perp}])\right\\}+\alpha_{2}\left(\sum_{j=0}^{i}\frac{(j-i)}{n}\left\langle\nabla^{2}\tilde{H}_{2}(\Pi\sigma_{j-1}),(\Pi
v_{j})(\Pi v_{j})^{\mathsf{T}}\right\rangle\right)\,,$

where we intrinsically used the fact that the contributions of the higher-
order derivatives are suppressed under the HES distributions. Given the
normalizations above, this can be simplified further to purely depend on
$\alpha_{2}$ as,

$\left\\{\sqrt{2\alpha_{2}\left\lVert\sigma\right\rVert^{2}_{2}}\mathsf{GOE}(\operatorname{Tr}[\Pi])\oplus\sqrt{\alpha_{2}}\mathsf{GOE}(\operatorname{Tr}[\Pi^{\perp}])\right\\}+\alpha_{2}\left(\sum_{j=0}^{i}\frac{(j-i)}{n}\left\langle\nabla^{2}\tilde{H}_{2}(\Pi\sigma_{j-1}),(\Pi
v_{j})(\Pi v_{j})^{\mathsf{T}}\right\rangle\right)$

Since $\left\lVert\sigma\right\rVert^{2}_{2}<\frac{1}{2}$ for the first half
of the iterates and the energy from the
$\alpha_{2}\tilde{H}_{2}(\Pi\sigma)\operatorname{\mathrm{Id}}_{n}$ term
contributes equally to $\Pi$ and $\Pi^{\perp}$ (while being independent of
$v$), it is clear that the SoS HES certificates for the Hessian will be
dominated by the signals from the $\Pi^{\perp}$ space.

##### Suboptimality of local algorithm

Since the local algorithm can only see the as-is extended Hamiltonian, for the
first few steps it will _only_ follow the $\Pi$-subspace (because of the
reasoning above). Note that this implies that the algorithm will go into the
wrong subspace for a constant fraction of the time, and output a solution
strictly worse than the HES relaxation. ∎

Note: The only missing part of the argumentation presented is that the
certifiable threshold rank condition of Theorem 1.1 is proxied for using the
extremal point in the support of the bulk of the spectrum. Given that the
eigenspectrum involves a block-decomposed sum of (independent) $\mathsf{GOE}$
matrices and Lemma 6.1 already certifies the $p$-Schatten norm of matrices
that are Wigner-like, this is a reasonable proxy. To see why, note that,

$\displaystyle\left\lVert\frac{1}{n}\nabla^{2}\tilde{H}(\sigma)\right\rVert^{p}_{p}$
$\displaystyle=\|\frac{1}{n}\left\\{\sqrt{2\alpha_{2}\left\lVert\sigma\right\rVert^{2}_{2}}\mathsf{GOE}(\operatorname{Tr}[\Pi])\oplus\sqrt{\alpha_{2}}\mathsf{GOE}(\operatorname{Tr}[\Pi^{\perp}])\right\\}+$
$\displaystyle\alpha_{2}\left(\sum_{j=0}^{i}\frac{(j-i)}{n}\left\langle\nabla^{2}\tilde{H}_{2}(\Pi\sigma_{j-1}),(\Pi
v_{j})(\Pi
v_{j})^{\mathsf{T}}\right\rangle\right)\operatorname{\mathrm{Id}}_{n}\|^{p}_{p}$
$\displaystyle=\frac{1}{n^{p}}\sum_{k=0}^{p}\left(\alpha_{2}\tilde{H}_{2}(\sigma)\right)^{p-k}n\left\lVert\sqrt{2\alpha_{2}\left\lVert\sigma\right\rVert^{2}_{2}}\mathsf{GOE}(\operatorname{Tr}[\Pi])\oplus\sqrt{\alpha_{2}}\mathsf{GOE}(\operatorname{Tr}[\Pi^{\perp}])\right\rVert^{k}_{k}$
$\displaystyle=\frac{1}{n^{p-1}}\sum_{k=0}^{p}\left(\alpha_{2}\tilde{H}_{2}(\sigma)\right)^{p-k}\left(\left\lVert\sigma\right\rVert^{2k}_{2}\left\lVert\alpha_{2}\Pi\nabla^{2}\tilde{H}_{2}(\Pi\sigma)\Pi\right\rVert^{k}_{k}+\left\lVert\alpha_{4}\Pi^{\perp}\nabla^{2}\tilde{H}_{4}(\Pi^{\perp}\sigma)\Pi^{\perp}\right\rVert^{k}_{k}\right)\,,$

at which point one can invoke Lemma 6.1 to bound the Schatten norms of each of
the individual $\mathsf{GOE}$ components and then use the HES constraint on
the moment matrices of $\sigma$ to apply a $(p,q)$-Hölder’s inequality as in
Corollary 6.3. The fact that operator norm concentrates follows from standard
arguments [vershynin2018high, Corollary 4.4.8].

###  8.3 More involved SoS proofs

For this example, a SoS proof of stepwise boundedness will require upper
bounding the value gained over all of a set of possible choices that the
algorithm could’ve made in previous steps. Here, the Hamiltonian will involve
the direct sum of two usual spin glass Hamiltonians in different subspaces,
along with an interaction term between the two subspaces. These terms will be
set up so that the optimal choice will be to evenly balance $\sigma$ between
the two subspaces at every step, but showing an upper bound will require
considering every possible choice of how to split $\sigma$ between the two
subspaces at each step.

This will only be a very rough sketch of a possible argument, depending on the
following assumption:

######  Assumption 8.3.

Let $A\in\mathbb{R}^{(n/2)\times(n/2)}$ and
$C\in\mathbb{R}^{(n/2)\times(n/2)}$ be random dense symmetric-matrix-valued
polynomials of $\sigma$ and $\tau$ respectively, and
$B\in\mathbb{R}^{(n/2)\times(n/2)}$ be a random dense asymmetric-random-valued
matrix polynomial of both $\sigma$ and $\tau$. In this way, $A$, $B$, and $C$
will be GOE-distributed with respect to the randomness of the polynomial
coefficients and a HES distribution over $(\sigma,\tau)$. Let $a$, $b$, and
$c$, be the expected operator norms of $A$, $B$, and $C$, which may depend on
$\left\lVert\sigma\right\rVert_{2}^{2}$ and
$\left\lVert\tau\right\rVert_{2}^{2}$.

There exists some function $f(\mu^{2},a,b,c)$ that expresses the with-high-
probability maximum over $(u,v)$ satisfying $\left\lVert
u\right\rVert_{2}^{2}+\left\lVert v\right\rVert_{2}^{2}=1$ and $\left\lVert
u\right\rVert_{2}^{2}=\mu^{2}$ of

$\left(u\;\;v\right)M\left(\\!\\!\begin{array}[]{c}u\\\
v\end{array}\\!\\!\right):=\left(u\;\;v\right)\left(\begin{array}[]{cc}A&B\\\
B^{\mathsf{T}}&C\end{array}\right)\left(\\!\\!\begin{array}[]{c}u\\\
v\end{array}\\!\\!\right)$

There also exist $g_{X}(\mu^{2},a,b,c)$ for $X\in\\{A,B,C\\}$ that express the
with-high-probability values of $u^{\mathsf{T}}Au$, $2u^{\mathsf{T}}Bv$, and
$v^{\mathsf{T}}Cv$, respectively, for the aforementioned maximizing $(u,v)$.

The assumption is that there are SoS proofs that capture the fact that the
quantities described by $f(\mu^{2},a,b,c)$ and $g_{X}(\mu^{2},a,b,c)$ actually
concentrate strongly at the true values of those functions.

Furthermore, there is $b^{2}/(ac)$ large enough so that $f(\mu^{2},a,b,c)$ is
maximized at $\mu^{2}=1/2$.

One might hope for the existence of this SoS proof because there are
combinatorial arguments that compute the spectral moments of block-structured
Wigner-type random matrices via a trace power method [zhu2020graphon] which
has structure very similar to the proof of Lemma 6.1.

######  Example 8.4 (Shared-subspace projector).

Let the Hamiltonian $H$ be,

$H(\sigma)=\alpha_{4}H_{4}(\Pi\sigma)+\alpha_{6}H_{6}(\Pi\sigma,\Pi^{\perp}\sigma)+\alpha_{8}H_{8}(\Pi^{\perp}\sigma)\,,$

where $H_{6}(\rho,\tau)=\left\langle g_{6},\rho^{\otimes
2}\otimes\tau^{\otimes 4}\right\rangle$ for $g_{6}$ an order-6 tensor of
i.i.d. unit Gaussian entries. Then, there exists a choice of
$\alpha_{4},\alpha_{6},\alpha_{8}\geqslant 0$ which causes Algorithm 2 to be
provably suboptimal, but the HES SoS relaxation does better with 8.3.

The argument for the failure of the local algorithm will be very similar as
for the previous example, so we focus here on the success of the HES
hierarchy, especially on the stepwise upper bound on the extended Hamiltonian.
Let the extension be

$\tilde{H}(\sigma)=\alpha_{4}\left\lVert\sigma\right\rVert_{2}^{4}H_{4}(\Pi\sigma)+\alpha_{6}\left\lVert\sigma\right\rVert_{2}^{2}H_{6}(\Pi\sigma,\Pi^{\perp}\sigma)+\alpha_{8}H_{8}(\Pi^{\perp}\sigma)\,.$

Then its Hessian is

$\displaystyle\nabla^{2}\tilde{H}(\sigma)={}$
$\displaystyle\alpha_{4}\left(\left\lVert\sigma\right\rVert_{2}^{4}\nabla^{2}H_{4}(\Pi\sigma)+12\left\lVert\sigma\right\rVert_{2}^{2}H_{4}(\Pi\sigma)\operatorname{\mathrm{Id}}\right)$
$\displaystyle+2\alpha_{6}H_{6}(\Pi\sigma,\Pi^{\perp}\sigma)\operatorname{\mathrm{Id}}$
$\displaystyle+\alpha_{6}\left\lVert\sigma\right\rVert_{2}^{2}\left(2\nabla_{1}^{2}H_{6}(\Pi\sigma,\Pi^{\perp}\sigma)+12\nabla_{2}^{2}H_{6}(\Pi\sigma,\Pi^{\perp}\sigma)+16\nabla_{1}\nabla_{2}H_{6}(\Pi\sigma,\Pi^{\perp}\sigma)\right)$
$\displaystyle+\alpha_{8}\nabla^{2}H_{8}(\Pi^{\perp}\sigma)+\text{[rank-1
terms]}\,,$

where $\nabla_{1}$ refers to the gradient in the first argument of a multi-
argument function and similarly for $\nabla_{2}$, and the rank-1 terms make
negligible contribution over an HES distribution.

Now we fit this into 8.3, characterizing each of these terms by their norm and
whether they are part of $A$, $B$, $C$, or constant with respect to $v_{i}$.
We halve the contributions to $B$ due to the fact that $B$ needs to appear
twice.

Term | Norm | Type
---|---|---
$\alpha_{4}\left\lVert\sigma\right\rVert_{2}^{4}\nabla^{2}H_{4}(\Pi\sigma)$ | $\sqrt{6}\,\alpha_{4}\left\lVert\sigma\right\rVert_{2}^{4}\left\lVert\Pi\sigma\right\rVert_{2}^{2}$ | $A$
$12\alpha_{4}\left\lVert\sigma\right\rVert_{2}^{2}H_{4}(\Pi\sigma)\operatorname{\mathrm{Id}}$ | $12\alpha_{4}\left\lVert\sigma\right\rVert_{2}^{2}H_{4}(\Pi\sigma)$ | constant
$2\alpha_{6}H_{6}(\Pi\sigma,\Pi^{\perp}\sigma)\operatorname{\mathrm{Id}}$ | $2\alpha_{6}H_{6}(\Pi\sigma,\Pi^{\perp}\sigma)$ | constant
$2\alpha_{6}\left\lVert\sigma\right\rVert_{2}^{2}\nabla_{1}^{2}H_{6}(\Pi\sigma,\Pi^{\perp}\sigma)$ | $\sqrt{2}\alpha_{6}\left\lVert\sigma\right\rVert_{2}^{2}\left\lVert\Pi^{\perp}\sigma\right\rVert_{2}^{4}$ | $A$
$12\alpha_{6}\left\lVert\sigma\right\rVert_{2}^{2}\nabla_{2}^{2}H_{6}(\Pi\sigma,\Pi^{\perp}\sigma)$ | $\sqrt{12}\alpha_{6}\left\lVert\sigma\right\rVert_{2}^{2}\left\lVert\Pi\sigma\right\rVert_{2}^{2}\left\lVert\Pi^{\perp}\sigma\right\rVert_{2}^{2}$ | $C$
$16\alpha_{6}\left\lVert\sigma\right\rVert_{2}^{2}\nabla_{1}\nabla_{2}H_{6}(\Pi\sigma,\Pi^{\perp}\sigma)$ | $\sqrt{8}\alpha_{6}\left\lVert\sigma\right\rVert_{2}^{2}\left\lVert\Pi\sigma\right\rVert_{2}^{1}\left\lVert\Pi^{\perp}\sigma\right\rVert_{2}^{3}$ | $B$
$\alpha_{8}\nabla^{2}H_{8}(\Pi^{\perp}\sigma)$ | $\sqrt{28}\left\lVert\Pi^{\perp}\sigma\right\rVert_{2}^{6}$ | $C$

At this point, we invoke 8.3 and unfurl it over $i$ steps to get a formula for
the total expected value of the Hamiltonian
$H(\sigma_{i})\approx\sum_{j\in[i]}v_{j}^{\mathsf{T}}\nabla^{2}H(\sigma_{j-1})v_{j}$
at the end of $i$ steps, as a function of $\left\lVert\Pi
v_{j}\right\rVert_{2}^{2}$ for each $j\in[i]$.

Showing that this formula is maximized when $\left\lVert\Pi
v_{j}\right\rVert_{2}^{2}=1/2$ for all $j$ will involve a “payment plan”
argument where, while the value gain from the $H_{8}$ term is immediate, the
value gain from the $H_{4}$ and $H_{6}$ terms are spread out over all future
steps via the $\operatorname{\mathrm{Id}}$ terms in the Hessian. The task will
be to show that the extra value from the “payment plans” when taking
$\left\lVert\Pi v_{j}\right\rVert_{2}^{2}>1/2$ fails to offset the future
losses to the $H_{8}$ term, since the $H_{8}$ term scales with a relatively
high power of $\left\lVert\Pi^{\perp}\sigma\right\rVert_{2}^{2}$, and
similarly but conversely when taking $\left\lVert\Pi
v_{j}\right\rVert_{2}^{2}<1/2$.

Note: The parameters in the examples above (the subspaces projected to) are
carefully balanced (for instance, $\\{\Pi,\Pi^{\perp}\\}$ are chosen so that
$\operatorname{Tr}[\Pi]=\operatorname{Tr}[\Pi^{\perp}]$) to make the argument
easier, but we expect that the argument generalizes to unbalanced choices as
well, in part by using the $R$-transform from free probability theory
[speicher2019lecture] to more precisely calculate spectral distributions of
sums of matrices. Given this flexibility of choice one has to “fool” the local
algorithm in ways that the SoS HES relaxation seems to be robust to, it leads
to a natural conjecture - That the SoS HES relaxation can “peel” out
correlations between the entries (when drawn from some sub-gaussian
distribution), provided these correlations are arbitrarily well-approximated
by a low-degree polynomial. This conjecture is formally stated in 1.5.

##  9 Discussion

###  9.1 Hardness equivalence between LDPs and $O(1)$-degree SoS?

In the introduction, we described the desire in the literature to relate
various meta-algorithmic families, such as LDPs and SoS hierarchies. To
demonstrate two-way reductions between these families for the family of random
optimization problems of interest, it is not enough to merely show that the
(for instance) SoS HES hierarchy can match the various LDP algorithms in use
on the family of problems, one must also demonstrate that any SoS hierarchy
cannot do _better_ in terms of outputting better solutions (under the
utilization of some reasonable rounding schemes [barak2014sum, ss24rounding]).

Partly, the “hardness” half of this task would involve demonstrating that,
under the presence of various so-called overlap-gap properties
[chen2019suboptimality, gamarnik2020low, chou2022limitations, huang2021tight,
jones2022random] (OGPs), there exist _no_ $O(1)$-degree SoS relaxations that
perform above a certain algorithmic threshold [huang2021tight, Eq. (1.9) & Eq.
(2.5)][jones2022random, Corollary 6.9]. Even more ambitiously, one would
perhaps attempt to rule out the possibility of $\operatorname{poly}(n)$ sized
semi-definite programs (SDPs) in the presence of OGPs altogether.

However, as it stands, this seems to be a rather daunting technical challenge
and a good first step in this direction would be to demonstrate the overlap-
concentration of (certain families of) low-degree polynomial algorithms
[jones2022random, Conjecture 7.1] on mean-field spin glasses and sparse random
Max-CSPs.

###  9.2 SoS algorithms for the SK model & bounded-degree CSPs under fRSB?

There is a precise way to modify the parameterized HES distribution, so that
the SoS relaxation can be adapted to work when the solution domain is the
discrete hypercube. This requires certifying a “mirror” hessian bound and
defining the SoS program over a newly parameterized HES distribution
[auffinger2015parisi, montanari2021optimization, panchenko2013parisi]. Doing
this would require generalizing the graph matrix framework [ahn2016graph] to
the setting of random orthogonal matrices that are conjugated by diagonal
matrices with _non_ i.i.d. entries. This is a significant technical challenge,
but ideas from asymptotic freeness introduced in the free probability theory
literature in conjunction with appropriately modified cumulant constraints
allow the current SoS relaxation template to work for spin glasses on the
hypercube as well, and this is demonstrated in forthcoming work.

While an AMP algorithm is already known for the SK model
[montanari2021optimization], it would be interesting to know if a SoS
algorithm using some parameterized family of certificates can be made to work
in the setting of _all_ bounded degree CSPs. Adapting the current SoS program
for the sparse setting, however, would require surmounting the formidable
technical challenge of bounding the operator norms of _extremely_ sparse graph
matrices.

###  9.3 SoS certificate complexity of the Ghirlanda-Guerra identities?

Note that the success of the HES SoS relaxation given in our work crucially
relies on the information about the distribution of the eigenspectrum of the
Hessian and its scaling equivalence matching the optimal value given by the
Chen and Sen formula (Lemma 3.7). More precisely, the precise optimal value is
known to us by knowledge that is _external_ to the SoS solver. A natural
question then is whether one can certify key mathematical properties of many
spin glass models that are known to be the reason why the ground state
configurations reside in an ultrametric space. As mentioned before, due to a
landmark result of Panchenko [panchenko2010connection, panchenko2013parisi],
it is known that ultrametricity is implied by the so-called Ghirlanda-Guerra
(GG) identities [ghirlanda1998general]. As a result, a first natural question
along this direction is asking whether there even _exist_ SoS certificates to
check whether an instance of certain spin glass models satisfies some low-
degree weakening of the GG identities, written over _some_ parameterizable
step distribution.

######  Open Question 9.1 (SoS Certificate for “weakened” Ghirlanda-Guerra
Identities of Low-Degree Moments of the Gibbs Distribution).

Does there exist a SoS relaxation with degree parametrized by the support of
the overlap CDF $\mu$ (and $n$) as defined in the Parisi functional
$\mathcal{P}_{\beta}(\mu)$ that can certify whether an instance $H_{n}(.)$ of
a mean-field spin glass (approximately) satisfies the GG identities on the
_low degree_ moments of some parameterizable distribution (to proxy those of
the Gibbs distribution)? Are there other weakenings of the GG identities that
a SoS relaxation can certify?

Note that it is unclear how to even state this problem in the sum-of-squares
proof system, let alone find a certificate for it, as the Ghirlanda-Guerra
identities are a property satisfied by the Gibbs measure of (typical)
instances of spin glass models. However, if it is possible to understand the
low-degree cumulants of the Gibbs precisely, one may state a relaxed version
of this property over the low-degree moments of the Gibbs distribution. An
answer to the question above sheds light on the power of the SoS proof system
to “recognize” when solutions to a problem instance satisfy a certain
geometric structure (albeit under a different relaxation), which in this case
is an ultrametric structure. In a regime (if any) in which certificates of
low-complexity exist, one can conceivably show the existence of a SoS
algorithm that relies _even less_ on external knowledge about the problem
structure.

###  9.4 Consequences for quantum optimization

By a union of the results in [chou2022limitations, jones2022random,
alaoui2023local, huang2021tight] along with the conjecture that LDPs, AMP and
$O(1)$-degree SoS should be equivalent in performance on mean-field spin
glasses and sparse random Max-CSPs, a compelling picture is emerging that
there is _no_ scope for demonstrable quantum advantage (optimization) on these
family of problems [farhi2014quantum, farhi2019quantum] in the fRSB regime. In
the OGP regime, there also seems to be little scope for advantage (at least
for shallow quantum circuits [chou2022limitations, jones2022random,
anshu2023concentration]).

In fact, in the presence of fRSB, the outlined family of classical algorithms
exhibit (or are expected to exhibit) PTAS’ for these models. This provides
formal ground to “bypass” the claim that $QAOA_{p}$ does better than the
standard SDP relaxation [farhi2019quantum] on optimizing the SK model, since
this claim’s validity critically relies on the fact that no _current_ SDP
relaxation is known that exhibits a PTAS on the SK model. However, this claim
would cease to be true if the outlined conjecture of equivalence holds.

Based on the combination of known rigorous theoretical results and the
emerging average-case complexity landscape, it seems that the _only_ way to
have rigorous theoretical evidence for demonstrable quantum optimization
advantage on these family of problems is to demonstrate an efficient quantum
algorithm that can provably beat the “OGP” algorithmic threshold
[huang2021tight, Eq. (1.9) & Eq. (2.5)][jones2022random, Corollary 6.9].

###  9.5 Bridging the HES and Lasserre Hierarchies

We observe that Lemma 1.3 suggests that the HES approach might generally work
when there are no large spectral weights in the higher-degree components of
the input polynomial. The existence of these heavy spectral directions would
suggest that techniques based on reweighing or lifting into higher-dimensional
spaces would be non-trivial [hopkins2016fast, ma2016polynomial]. It would be
interesting to bridge the gap between these two regimes.

As well, it would be interesting to know under what conditions a low-bit-
complexity proof is guaranteed in the HES hierarchy, as this would ensure the
relaxation can be solved in polynomial-time in the general case
[raghavendra2017bit].

## Acknowledgements

We thank Alexandra Kolla, Boaz Barak & Peter J. Love for helpful comments and
feedback. We are grateful to Sam Hopkins for pointing out a simplification of
the notation for pseudo-distributions over distributions and for various
helpful discussions. We are also indebted to David Steurer for extended
insightful discussions while hosting us at ETH Zurich during Summer 2022. We
would like to thank Chris Jones for helpful comments on graph matrix norm
bounds and certifications. We would also like to thank Pravesh Kothari for
pointing us to an error in an earlier version of this argument and for
providing helpful feedback on a prior draft. We are grateful to the anonymous
reviewers for suggestions on how to improve presentation and include a
discussion on input-oblivious rounding, and we’d like to specially thank Kunal
Marwaha for giving detailed feedback and suggestions on the introduction.

Much of this work began after a long series of discussions with Antares Chen
in Fall 2021 and with Tommaso d’Orsi in Spring and Summer 2022, seeding the
development of these ideas, for which we are extremely grateful.

JS would like to thank Luca Trevisan for instigating his introduction to spin
glass theory and for discussions on SoS and local-improvement-based methods.
JSS would like to thank Peter J. Love for excellent mentorship and sharing
thoughts on “algorithmic equivalence” which (in part) motivated this paper.
JSS partly worked on this project as a visiting scholar at Bocconi University,
hosted by Luca Trevisan during Summer 2022.

Lastly, JS & JSS would like to extend many thanks to Prof. SueYeon Chung for
graciously hosting them at the Flatiron Institute in New York City over
various periods in Summer-2023 and Fall-2023 where the final stages of this
paper were completed.

JSS would like to dedicate this work and the companion work [ss24rounding] to
his late grandfather, Dr. Rashpal Malhotra.

## References

  * [ABS15] Sanjeev Arora, Boaz Barak, and David Steurer.  Subexponential algorithms for unique games and related problems.  Journal of the ACM (JACM), 62(5):1–25, 2015.
  * [AC15] Antonio Auffinger and Wei-Kuo Chen.  The Parisi formula has a unique minimizer.  Communications in Mathematical Physics, 335(3):1429–1444, 2015.
  * [ADPM12] Luigi Ambrosio, Giuseppe Da Prato, and Andrea Mennucci.  Introduction to measure theory and integration, volume 10.  Springer Science & Business Media, 2012.
  * [AGZ10] Greg W Anderson, Alice Guionnet, and Ofer Zeitouni.  An introduction to random matrices.  Cambridge university press, 2010.
  * [AM23] Anurag Anshu and Tony Metger.  Concentration bounds for quantum states and limitations on the QAOA from polynomial approximations.  Quantum, 7:999, 2023.
  * [AMP21] Kwangjun Ahn, Dhruv Medarametla, and Aaron Potechin.  Graph matrices: norm bounds and applications.  arXiv preprint arXiv:1604.03423v5, 2021.
  * [AMS21] Ahmed El Alaoui, Andrea Montanari, and Mark Sellke.  Optimization of mean-field spin glasses.  The Annals of Probability, 49(6):2922–2960, 2021.
  * [AMS22] Ahmed El Alaoui, Andrea Montanari, and Mark Sellke.  Sampling from the Sherrington-Kirkpatrick Gibbs measure via algorithmic stochastic localization.  In 2022 IEEE 63rd Annual Symposium on Foundations of Computer Science (FOCS), pages 323–334. IEEE, 2022.
  * [AMS23] Ahmed El Alaoui, Andrea Montanari, and Mark Sellke.  Local algorithms for maximum cut and minimum bisection on locally treelike regular graphs of large degree.  Random Structures & Algorithms, 63(3):689–715, 2023.
  * [ASZ20] Gérard Ben Arous, Eliran Subag, and Ofer Zeitouni.  Geometry and temperature chaos in mixed spherical spin glasses at low temperature: the perturbative regime.  Communications on Pure and Applied Mathematics, 73(8):1732–1828, 2020.
  * [AT09] Robert J Adler and Jonathan E Taylor.  Random fields and geometry.  Springer Science & Business Media, 2009.
  * [BAH+22] Afonso S Bandeira, Ahmed El Alaoui, Samuel B Hopkins, Tselil Schramm, Alexander S Wein, and Ilias Zadik.  The Franz-Parisi criterion and computational trade-offs in high dimensional statistics.  Advances in Neural Information Processing Systems, 35:33831–33844, 2022.
  * [Bar14] Boaz Barak.  Sum of squares upper bounds, lower bounds, and open questions.  Lecture notes, 2014.
  * [BBH+12] Boaz Barak, Fernando GSL Brandao, Aram W Harrow, Jonathan Kelner, David Steurer, and Yuan Zhou.  Hypercontractivity, sum-of-squares proofs, and their applications.  In Proceedings of the forty-fourth annual ACM symposium on Theory of computing, pages 307–326, 2012.
  * [BBK+21] Mitali Bafna, Boaz Barak, Pravesh K Kothari, Tselil Schramm, and David Steurer.  Playing unique games on certified small-set expanders.  In Proceedings of the 53rd Annual ACM SIGACT Symposium on Theory of Computing, pages 1629–1642, 2021.
  * [BDH+20] Ainesh Bakshi, Ilias Diakonikolas, Samuel B Hopkins, Daniel Kane, Sushrut Karmalkar, and Pravesh K Kothari.  Outlier-robust clustering of gaussians and other non-spherical mixtures.  In 2020 IEEE 61st Annual Symposium on Foundations of Computer Science (FOCS), pages 149–159. IEEE, 2020.
  * [Ber12] S. Bernstein.  D´emoistration du th’eor‘eme de Weierstrass fond´ee sur le calcul des probabilit´es.  Communications de la Soci´et´e math´ematique de Kharkow., 13(1):1–2, 1912.
  * [BGL17] Vijay Bhattiprolu, Venkatesan Guruswami, and Euiwoong Lee.  Sum-of-Squares Certificates for Maxima of Random Tensors on the Sphere.  In Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques (APPROX/RANDOM 2017). Schloss Dagstuhl-Leibniz-Zentrum fuer Informatik, 2017.
  * [BKS14] Boaz Barak, Jonathan A Kelner, and David Steurer.  Rounding sum-of-squares relaxations.  In Proceedings of the forty-sixth annual ACM symposium on Theory of computing, pages 31–40, 2014.
  * [BMO+15] Boaz Barak, Ankur Moitra, Ryan O’Donnell, Prasad Raghavendra, Oded Regev, David Steurer, Luca Trevisan, Aravindan Vijayaraghavan, David Witmer, and John Wright.  Beating the Random Assignment on Constraint Satisfaction Problems of Bounded Degree.  Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques, page 110, 2015.
  * [BMR21] Jess Banks, Sidhanth Mohanty, and Prasad Raghavendra.  Local statistics, semidefinite programming, and community detection.  In Proceedings of the 2021 ACM-SIAM Symposium on Discrete Algorithms (SODA), pages 1298–1316. SIAM, 2021.
  * [BR22] Andrei A Bulatov and Akbar Rafiey.  On the complexity of csp-based ideal membership problems.  In Proceedings of the 54th Annual ACM SIGACT Symposium on Theory of Computing, pages 436–449, 2022.
  * [BRS11] Boaz Barak, Prasad Raghavendra, and David Steurer.  Rounding semidefinite programming hierarchies via global correlation.  In 2011 ieee 52nd annual symposium on foundations of computer science, pages 472–481. IEEE, 2011.
  * [BS14] Boaz Barak and David Steurer.  Sum-of-squares proofs and the quest toward optimal algorithms.  In Proceedings of International Congress of Mathematicians (ICM), 2014.
  * [BS16] Boaz Barak and David Steurer.  Proofs, beliefs, and algorithms through the lens of sum-of-squares.  Course notes: http://www. sumofsquares. org/public/index. html, 1, 2016.
  * [CGPR19] Wei-Kuo Chen, David Gamarnik, Dmitry Panchenko, and Mustazee Rahman.  Suboptimality of local algorithms for a class of max-cut problems.  Annals of Probability, 47(3):1587–1618, 2019.
  * [Chi11] Theodore S Chihara.  An introduction to orthogonal polynomials.  Courier Corporation, 2011.
  * [CHM23] Antares Chen, Neng Huang, and Kunal Marwaha.  Local algorithms and the failure of log-depth quantum advantage on sparse random CSPs.  arXiv preprint arXiv:2310.01563, 2023.
  * [Cim12] J. Cimprič.  Real algebraic geometry for matrices over commutative rings.  Journal of Algebra, 359:89 – 103, 2012.
  * [CLSS22] Chi-Ning Chou, Peter J Love, Juspreet Singh Sandhu, and Jonathan Shi.  Limitations of local quantum algorithms on random max-k-xor and beyond.  In 49th International Colloquium on Automata, Languages, and Programming (ICALP 2022). Schloss Dagstuhl-Leibniz-Zentrum für Informatik, 2022.
  * [CMŚS07] Benoît Collins, James A Mingo, Piotr Śniady, and Roland Speicher.  Second order freeness and fluctuations of random matrices. III: Higher order freeness and free cumulants.  Documenta Mathematica, 12:1–70, 2007.
  * [Con18] Francesco Concetti.  The full replica symmetry breaking in the Ising spin glass on random regular graph.  Journal of Statistical Physics, 173(5):1459–1483, 2018.
  * [Con19] Francesco Concetti.  The full replica symmetry breaking solution in mean-field spin glass models.  arXiv preprint arXiv:1911.00557, 2019.
  * [COP19] Amin Coja-Oghlan and Will Perkins.  Spin systems on Bethe lattices.  Communications in Mathematical Physics, 372:441–523, 2019.
  * [CPRT22] Flavio Chierichetti, Alessandro Panconesi, Giuseppe Re, and Luca Trevisan.  Spectral Robustness for Correlation Clustering Reconstruction in Semi-Adversarial Models.  In International Conference on Artificial Intelligence and Statistics, pages 10852–10880. PMLR, 2022.
  * [CPS18] Wei-Kuo Chen, Dmitry Panchenko, and Eliran Subag.  Generalized TAP Free Energy.  Communications on Pure and Applied Mathematics, 2018.
  * [CPS21] Wei-Kuo Chen, Dmitry Panchenko, and Eliran Subag.  The generalized TAP free energy II.  Communications in Mathematical Physics, 381(1):257–291, 2021.
  * [CS92] Andrea Crisanti and H-J Sommers.  The sphericalp-spin interaction spin glass model: the statics.  Zeitschrift für Physik B Condensed Matter, 87(3):341–354, 1992.
  * [CS17] Wei-Kuo Chen and Arnab Sen.  Parisi formula, disorder chaos and fluctuation for the ground state energy in the spherical mixed p-spin models.  Communications in Mathematical Physics, 350(1):129–173, 2017.
  * [CST22] Antares Chen, Jonathan Shi, and Luca Trevisan.  Cut Sparsification of the Clique Beyond the Ramanujan Bound: A Separation of Cut Versus Spectral Sparsification.  In Proceedings of the 2022 Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), pages 3693–3731. SIAM, 2022.
  * [DB59] Louis De Branges.  The stone-weierstrass theorem.  Proceedings of the American Mathematical Society, 10(5):822–824, 1959.
  * [DGJ+10] Ilias Diakonikolas, Parikshit Gopalan, Ragesh Jaiswal, Rocco A Servedio, and Emanuele Viola.  Bounded independence fools halfspaces.  SIAM Journal on Computing, 39(8):3441–3462, 2010.
  * [DKK+22] Ilias Diakonikolas, Daniel M Kane, Sushrut Karmalkar, Ankit Pensia, and Thanasis Pittas.  Robust sparse mean estimation via sum of squares.  In Conference on Learning Theory, pages 4703–4763. PMLR, 2022.
  * [FGG14] Edward Farhi, Jeffrey Goldstone, and Sam Gutmann.  A quantum approximate optimization algorithm.  arXiv preprint arXiv:1411.4028, 2014.
  * [FGGZ19] Edward Farhi, Jeffrey Goldstone, Sam Gutmann, and Leo Zhou.  The quantum approximate optimization algorithm and the sherrington-kirkpatrick model at infinite size.  arXiv preprint arXiv:1910.08187, 2019.
  * [Fol01] Gerald B Folland.  How to integrate a polynomial over a sphere.  The American Mathematical Monthly, 108(5):446–448, 2001.
  * [For23] Peter J Forrester.  A review of exact results for fluctuation formulas in random matrix theory.  Probability Surveys, 20:170–225, 2023.
  * [Gam21] David Gamarnik.  The overlap gap property: A topological barrier to optimizing over random structures.  Proceedings of the National Academy of Sciences, 118(41):e2108492118, 2021.
  * [GG98] Stefano Ghirlanda and Francesco Guerra.  General properties of overlap probability distributions in disordered spin systems. Towards Parisi ultrametricity.  Journal of Physics A: Mathematical and General, 31(46):9149, 1998.
  * [GJ21] David Gamarnik and Aukosh Jagannath.  The overlap gap property and approximate message passing algorithms for $p$-spin models.  The Annals of Probability, 49(1):180–205, 2021.
  * [GJJ+20] Mrinalkanti Ghosh, Fernando Granha Jeronimo, Chris Jones, Aaron Potechin, and Goutham Rajendran.  Sum-of-squares lower bounds for sherrington-kirkpatrick via planted affine planes.  In 2020 IEEE 61st Annual Symposium on Foundations of Computer Science (FOCS), pages 954–965. IEEE, 2020.
  * [GJW20] David Gamarnik, Aukosh Jagannath, and Alexander S Wein.  Low-degree hardness of random optimization problems.  In 2020 IEEE 61st Annual Symposium on Foundations of Computer Science (FOCS), pages 131–140. IEEE, 2020.
  * [GLR18] Rong Ge, Holden Lee, and Andrej Risteski.  Simulated tempering langevin monte carlo ii: An improved proof using soft markov chain decomposition.  arXiv preprint arXiv:1812.00793, 2018.
  * [Gri01] Dima Grigoriev.  Linear lower bound on degrees of Positivstellensatz calculus proofs for the parity.  Theoretical Computer Science, 259(1-2):613–622, 2001.
  * [GS14] David Gamarnik and Madhu Sudan.  Limits of local algorithms over sparse random graphs.  In Proceedings of the 5th conference on Innovations in theoretical computer science, pages 369–376, 2014.
  * [Har13] Aram W Harrow.  The church of the symmetric subspace.  arXiv preprint arXiv:1308.6595, 2013.
  * [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 2017 IEEE 58th Annual Symposium on Foundations of Computer Science (FOCS), pages 720–731. IEEE, 2017.
  * [HL18] Samuel B Hopkins and Jerry Li.  Mixture models, robustness, and sum of squares proofs.  In Proceedings of the 50th Annual ACM SIGACT Symposium on Theory of Computing, pages 1021–1034, 2018.
  * [Hop18a] Sam Hopkins.  Clustering and Sum of Squares Proofs: Six Blog Posts on Unsupervised Learning.  https://www.samuelbhopkins.com/clustering.pdf, 2018.  2022-11-07.
  * [Hop18b] Samuel Hopkins.  Statistical inference and the sum of squares method.  PhD thesis, Cornell University, 2018.
  * [HS22] Brice Huang and Mark Sellke.  Tight lipschitz hardness for optimizing mean field spin glasses.  In 2022 IEEE 63rd Annual Symposium on Foundations of Computer Science (FOCS), pages 312–322. IEEE, 2022.
  * [HSSS16] Samuel B Hopkins, Tselil Schramm, Jonathan Shi, and David Steurer.  Fast spectral algorithms from sum-of-squares proofs: tensor decomposition and planted sparse vectors.  In Proceedings of the forty-eighth annual ACM symposium on Theory of Computing, pages 178–191, 2016.
  * [Hua22] Brice Huang.  Computational Hardness in Random Optimization Problems from the Overlap Gap Property.  PhD thesis, Massachusetts Institute of Technology, 2022.
  * [IS23] Misha Ivkov and Tselil Schramm.  Semidefinite programs simulate approximate message passing robustly.  arXiv preprint arXiv:2311.09017, 2023.
  * [JKO98] Richard Jordan, David Kinderlehrer, and Felix Otto.  The Variational Formulation of the Fokker–Planck Equation.  SIAM Journal on Mathematical Analysis, 29(1):1–17, 1998.
  * [JKR19] Vishesh Jain, Frederic Koehler, and Andrej Risteski.  Mean-field approximation, convex hierarchies, and the optimality of correlation rounding: a unified perspective.  In Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing, pages 1226–1236, 2019.
  * [JMSS23] Chris Jones, Kunal Marwaha, Juspreet Singh Sandhu, and Jonathan Shi.  Random Max-CSPs Inherit Algorithmic Hardness from Spin Glasses.  In 14th Innovations in Theoretical Computer Science Conference (ITCS 2023), volume 251 of Leibniz International Proceedings in Informatics (LIPIcs), pages 77:1–77:26. Schloss Dagstuhl – Leibniz-Zentrum für Informatik, 2023.
  * [JPR+22] Chris Jones, Aaron Potechin, Goutham Rajendran, Madhur Tulsiani, and Jeff Xu.  Sum-of-squares lower bounds for sparse independent set.  In 2021 IEEE 62nd Annual Symposium on Foundations of Computer Science (FOCS), pages 406–416. IEEE, 2022.
  * [KB21] Dmitriy Kunisky and Afonso S Bandeira.  A tight degree 4 sum-of-squares lower bound for the Sherrington–Kirkpatrick Hamiltonian.  Mathematical Programming, 190(1):721–759, 2021.
  * [KKK19] Sushrut Karmalkar, Adam Klivans, and Pravesh Kothari.  List-decodable linear regression.  Advances in neural information processing systems, 32, 2019.
  * [KSS18] Pravesh K Kothari, Jacob Steinhardt, and David Steurer.  Robust moment estimation and improved clustering via sum of squares.  In Proceedings of the 50th Annual ACM SIGACT Symposium on Theory of Computing, pages 1035–1046, 2018.
  * [KV05] Subhash Khot and Nisheeth K Vishnoi.  On the unique games conjecture.  In FOCS, volume 5, page 3, 2005.
  * [Lat06] Rafał Latała.  Estimates of moments and tails of Gaussian chaoses.  The Annals of Probability, 2006.
  * [LM21] Allen Liu and Ankur Moitra.  Settling the robust learnability of mixtures of gaussians.  In Proceedings of the 53rd Annual ACM SIGACT Symposium on Theory of Computing, pages 518–531, 2021.
  * [LM23] Euiwoong Lee and Pasin Manurangsi.  Hardness of Approximating Bounded-Degree Max 2-CSP and Independent Set on k-Claw-Free Graphs.  arXiv preprint arXiv:2309.04099, 2023.
  * [LRS15] James R Lee, Prasad Raghavendra, and David Steurer.  Lower bounds on the size of semidefinite programming relaxations.  In Proceedings of the forty-seventh annual ACM symposium on Theory of computing, pages 567–576, 2015.
  * [LS17] Zhixian Lei and Yueqi Sheng.  Sum of square proof for Brascamp-Lieb type inequality.  arXiv preprint arXiv:1710.01458, 2017.
  * [Mon19] A. Montanari.  Optimization of the Sherrington-Kirkpatrick Hamiltonian.  In 2019 IEEE 60th Annual Symposium on Foundations of Computer Science (FOCS), pages 1417–1433, Los Alamitos, CA, USA, nov 2019. IEEE Computer Society.
  * [MP01] Marc Mézard and Giorgio Parisi.  The Bethe lattice spin glass revisited.  The European Physical Journal B-Condensed Matter and Complex Systems, 20:217–233, 2001.
  * [MR20] Ankur Moitra and Andrej Risteski.  Fast convergence for Langevin diffusion with manifold structure.  arXiv preprint arXiv:2002.05576, 2020.
  * [MS17] James A Mingo and Roland Speicher.  Free probability and random matrices, volume 35.  Springer, 2017.
  * [MSS16] Tengyu Ma, Jonathan Shi, and David Steurer.  Polynomial-time tensor decompositions with sum-of-squares.  In 2016 IEEE 57th Annual Symposium on Foundations of Computer Science (FOCS), pages 438–446. IEEE, 2016.
  * [Pan10] Dmitry Panchenko.  A connection between the Ghirlanda–Guerra identities and ultrametricity.  The Annals of Probability, 38(1):327–347, 2010.
  * [Pan13a] Dmitry Panchenko.  The Parisi ultrametricity conjecture.  Annals of Mathematics, pages 383–393, 2013.
  * [Pan13b] Dmitry Panchenko.  The Ruelle Probability Cascades.  In The Sherrington-Kirkpatrick Model, pages 33–77. Springer, 2013.
  * [Pan15] Dmitry Panchenko.  Hierarchical exchangeability of pure states in mean field spin glass models.  Probability Theory and Related Fields, 161(3-4):619–650, 2015.
  * [Pan16] Dmitry Panchenko.  Structure of finite-RSB asymptotic Gibbs measures in the diluted spin glass models.  Journal of statistical physics, 162(1):1–42, 2016.
  * [Par79] Giorgio Parisi.  Infinite number of order parameters for spin-glasses.  Physical Review Letters, 43(23):1754, 1979.
  * [Par80] Giorgio Parisi.  A sequence of approximated solutions to the SK model for spin glasses.  Journal of Physics A: Mathematical and General, 13(4):L115, 1980.
  * [Par04] Pablo A Parrilo.  Sum of squares programs and polynomial inequalities.  In SIAG/OPT Views-and-News: A Forum for the SIAM Activity Group on Optimization, volume 15, pages 7–15, 2004.
  * [Rag08] Prasad Raghavendra.  Optimal algorithms and inapproximability results for every CSP?  In Proceedings of the fortieth annual ACM symposium on Theory of computing, pages 245–254, 2008.
  * [Rah17] Sharif Rahman.  Wiener–Hermite polynomial expansion for multivariate Gaussian probability measures.  Journal of Mathematical Analysis and Applications, 454(1):303–334, 2017.
  * [Ris16] Andrej Risteski.  How to calculate partition functions using convex programming hierarchies: provable bounds for variational methods.  In Conference on Learning Theory, pages 1402–1416. PMLR, 2016.
  * [RT23] Goutham Rajendran and Madhur Tulsiani.  Concentration of polynomial random matrices via Efron-Stein inequalities.  In Proceedings of the 2023 Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), pages 3614–3653. SIAM, 2023.
  * [Rue87] David Ruelle.  A mathematical reformulation of Derrida’s REM and GREM.  Communications in Mathematical Physics, 108(2):225–239, 1987.
  * [RV+17] Mustazee Rahman, Balint Virag, et al.  Local algorithms for independent sets are half-optimal.  The Annals of Probability, 45(3):1543–1577, 2017.
  * [RW17] Prasad Raghavendra and Benjamin Weitz.  On the Bit Complexity of Sum-of-Squares Proofs.  In 44th International Colloquium on Automata, Languages, and Programming (ICALP 2017). Schloss Dagstuhl-Leibniz-Zentrum fuer Informatik, 2017.
  * [RY20] Prasad Raghavendra and Morris Yau.  List decodable learning via sum of squares.  In Proceedings of the Fourteenth Annual ACM-SIAM Symposium on Discrete Algorithms, pages 161–180. SIAM, 2020.
  * [Sch] Tselil Schramm.  private communication.
  * [Sch20] Tselil Schramm.  Reconciling Statistical Queries and the Low-Degree Likelihood Ratio.  https://www.youtube.com/watch?v=TBKlAcaqmwo, 2020.  Talk given for the The 1W-MINDS Seminar, 2020-09-14.
  * [SK75] David Sherrington and Scott Kirkpatrick.  Solvable model of a spin-glass.  Physical review letters, 35(26):1792, 1975.
  * [Spe14] Joel Spencer.  Asymptopia, volume 71.  American Mathematical Soc., 2014.
  * [Spe19] Roland Speicher.  Lecture Notes on” Free Probability Theory”.  arXiv preprint arXiv:1908.08125, 2019.
  * [SS24] Juspreet Singh Sandhu and Jonathan Shi.  A Sum-of-Squares Hierarchy in the Absence of Pointwise Proofs II: Rounding High-Entropy Steps.  To Appear, 2024.
  * [Sub18] Eliran Subag.  Free energy landscapes in spherical spin glasses.  arXiv preprint arXiv:1804.10576, 2018.
  * [Sub21] Eliran Subag.  Following the Ground States of Full-RSB Spherical Spin Glasses.  Communications on Pure and Applied Mathematics, 74(5):1021–1044, 2021.
  * [Tal06] Michel Talagrand.  Free energy of the spherical mean field model.  Probability theory and related fields, 134(3):339–382, 2006.
  * [Tal10a] Michel Talagrand.  Construction of pure states in mean field models for spin glasses.  Probability theory and related fields, 148:601–643, 2010.
  * [Tal10b] Michel Talagrand.  Mean field models for spin glasses: Volume I: Basic examples, volume 54.  Springer Science & Business Media, 2010.
  * [Ter21] György Terdik.  Multivariate statistical methods: going beyond the linear.  Springer Nature, 2021.
  * [Ver18] Roman Vershynin.  High-dimensional probability: An introduction with applications in data science, volume 47.  Cambridge university press, 2018.
  * [Wei02] Eric W Weisstein.  Lambert W-function.  https://mathworld. wolfram. com/, 2002.
  * [Wei22] Alexander S Wein.  Optimal low-degree hardness of maximum independent set.  Mathematical Statistics and Learning, 4(3):221–251, 2022.
  * [Zhu20] Yizhe Zhu.  A graphon approach to limiting spectral distributions of Wigner-type matrices.  Random Structures & Algorithms, 56(1):251–279, 2020.

##  Appendix A Univariate Matrix Bernstein Polynomials

In this section, we define the _matrix_ bernstein basis polynomials over
symmetric matrices with eigenvalues in a compact interval
$I\subset\mathbb{R}$. The matrix Bernstein polynomial $B^{f}_{d}(M)$ of degree
$d$ approximating a continuous function $f:I\to\mathbb{R}$ will have the same
desirable properties as in the scalar regime when the underlying inputs are
symmetric matrices, including a low-degree SoS proof of positivity over the
Loewner order (Appendix B).

######  Definition A.1 (Matrix Bernstein Basis Polynomials).

Given a degree $d\geqslant 0$, the $d+1$ matrix Bernstein basis polynomials on
the interval $I=[a,b]\subset\mathbb{R}$ are defined as,

$B_{d,i}=\binom{d}{i}\frac{(M-a\operatorname{\mathrm{Id}})^{i}(b\operatorname{\mathrm{Id}}-M)^{d-i}}{(b-a)^{d}}\,,$
(A.1)

where $M$ is a symmetric matrix such that $a\operatorname{\mathrm{Id}}\preceq
M\preceq b\operatorname{\mathrm{Id}}$.

Given a step-function $S_{\alpha}:I=[a,b]\to\mathbb{R}$ with the step at
$\alpha\in(a,b)$, we define the linear ramp function
$L_{\alpha,\gamma}:[a,b]\to\mathbb{R}$ as a continuous approximation of the
step-function as follows,

$L_{\alpha,\gamma}(x)=\begin{cases}0&x\in[a,\alpha-\gamma]\,,\\\
x&x\in(\alpha-\gamma,\alpha+\gamma)\,,\\\
1&x\in[\alpha+\gamma,b]\,.\end{cases}$ (A.2)

Typically, $\gamma$ will be chosen to be significantly smaller than $\alpha$.
For the Hessian $\nabla^{2}_{\perp}H(\cdot)$, the linear ramp function
approximates its top $\delta n$-dimensional eigenspace with the following
choice: $a=-2\nu^{\prime\prime}(\cdot)^{1/2}$,
$b=2\nu^{\prime\prime}(\cdot)^{1/2}$,
$\alpha=2\nu^{\prime\prime}(\cdot)^{1/2}-2\varepsilon/3$ and
$\gamma=\varepsilon/3$. The choice of $\varepsilon$ and $\delta$ correspond to
[subag2021following, Lemma 3].

######  Proposition A.2 (Matrix Bernstein Approximation of
$L_{\alpha,\gamma}$).

Let $I=[a,b]\subset\mathbb{R}$. Then, given a linear-ramp function
$L_{\alpha,\gamma}:I\to\mathbb{R}$, the Bernstein polynomial
$B^{L_{\alpha,\gamma}}_{d}(\cdot)$ defined over symmetric matrices $M$ with
$a\operatorname{\mathrm{Id}}\preceq M\preceq b\operatorname{\mathrm{Id}}$ that
approximates it is given as,

$B^{L_{\alpha,\gamma}}_{d}(M):=\sum_{j=0}^{d}L_{\alpha,\gamma}\left(a+\frac{j(b-a)}{d}\right)\binom{d}{i}\frac{(M-a\operatorname{\mathrm{Id}})^{j}(b\operatorname{\mathrm{Id}}-M)^{d-j}}{(b-a)^{d}}\,.$
(A.3)

Furthermore, $B^{L_{\alpha,\gamma}}_{d}$ has the following properties:

1\. Symmetric Positive Semi-Definiteness:

$B^{L_{\alpha,\gamma}}_{d}(M)\succeq 0$.

2\. Operator Norm Bound:

$B^{L_{\alpha,\gamma}}_{d}(M)\preceq\operatorname{\mathrm{Id}}$.

3\. Close in Frobenius Norm:

Let $\Pi_{\alpha+\gamma,b}$ be the orthogonal projector to the subspace of $M$
with eigenvalues in $[\alpha+\gamma,b]$ and
$p=\frac{B^{L_{\alpha},\gamma}_{d}(M)}{\delta^{\prime}n}$. Then, for any
$d\geqslant\frac{8\nu^{\prime\prime}(q)^{1/2}}{\varepsilon^{3}}$,

$\left\lVert
p-\frac{1}{\delta^{\prime}n}\Pi_{\alpha+\gamma,b}\right\rVert_{F}\leqslant\sqrt{\frac{2\delta^{\prime\prime\prime}}{\delta^{\prime
2}}}\frac{1}{\sqrt{n}}\,,$

where,

$\delta^{\prime\prime\prime}=\max\left\\{\varepsilon^{2},\frac{3}{4}(\delta-\delta^{\prime\prime})\right\\}\,.$

###### Proof.

We prove each item in sequence.

1\. Symmetric Positive Semi-Definiteness:

This follows immediately from the following observations:

  * •

The value of the scaled linear ramp function
$\binom{d}{i}\frac{L_{\alpha,\gamma}(a+j(b-a)/d)}{(b-a)^{d}}\geqslant 0$ for
every $j\in[d]$ since the function is non-negative and $b>a$.

  * •

By invoking Lemma B.2 with
$p_{1}(M-a\operatorname{\mathrm{Id}})=(M-a\operatorname{\mathrm{Id}})^{j}$ and
$p_{2}(b\operatorname{\mathrm{Id}}-M)=(b\operatorname{\mathrm{Id}}-M)^{d-j}$,
we have that $p_{1}(\cdot)p_{2}(\cdot)$ is a symmetric and positive definite
matrix since $a\operatorname{\mathrm{Id}}\preceq M\preceq
b\operatorname{\mathrm{Id}}$.

  * •

A positive linear sum of symmetric and positive semi-definite matrices is a
symmetric positive semi-definite matrix.

2\. Operator Norm Bound:

Diagonalizing the symmetric matrix immediately implies that,

$B^{L_{\alpha,\gamma}}_{d}(M)=O\begin{bmatrix}B^{L_{\alpha,\gamma}}_{d}(\lambda_{1})&&\\\
&\ddots&\\\
&&B^{L_{\alpha,\gamma}}_{d}(\lambda_{n})\end{bmatrix}O^{\mathsf{T}}\,,$ (A.4)

where the Bernstein polynomials applied to the eigenvalue are as defined in
Equation 3.21. Then, by Definition 3.28,

$\displaystyle B^{L_{\alpha,\gamma}}_{d}(\lambda_{i})$
$\displaystyle=\sum_{j=0}^{d}L_{\alpha,\gamma}(a+j(b-a)/d)\binom{d}{j}\left(\frac{\lambda_{i}-a}{b-a}\right)^{j}\left(\frac{b-\lambda_{i}}{b-a}\right)^{d-j}\,,$
$\displaystyle\leqslant\sum_{j=0}^{d}\binom{d}{j}\left(\frac{\lambda_{i}-a}{b-a}\right)^{j}\left(\frac{b-\lambda_{i}}{b-a}\right)^{d-j}$
$\displaystyle=\sum_{j=0}^{d}\binom{d}{j}y^{j}(1-y)^{d-j}=(y+(1-y))^{d}$
$\displaystyle=1\,,$

where we used the facts that $L_{\alpha,\gamma}(a+j(b-a)/d)\leqslant 1$ for
every $\lambda_{i}\in[a,b]$ and the substitution
$y=\frac{\lambda_{i}-a}{b-a}\in[0,1]$. The last equality follows by the
binomial formula. This immediately implies that,

$B^{L_{\alpha,\gamma}}_{d}(M)\leqslant
O\operatorname{\mathrm{Id}}O^{\mathsf{T}}=\operatorname{\mathrm{Id}}\,.$

3\. Close in Frobenius Norm:

We apply [subag2021following, Lemma 3] repetitively to obtain
$0<\delta^{\prime\prime}<\delta^{\prime}<\delta<1$ corresponding to the
choices of $0<\varepsilon/3<2\varepsilon/3<\varepsilon<1$ respectively. Then,
applying Theorem 3.29 for $f=L_{\alpha,\gamma}$ with the choice of parameters
$a,b,\gamma$ stated above, we have,

$\displaystyle\left\lVert
p-\frac{1}{\delta^{\prime}n}\Pi_{\alpha+\gamma,b}\right\rVert_{F}$
$\displaystyle=\frac{1}{\delta^{\prime}n}\sqrt{\sum_{i=1}^{\delta^{\prime\prime}n}\varepsilon^{2}+\sum_{i=\delta^{\prime\prime}n+1}^{\delta^{\prime}n}(1/2+\varepsilon)^{2}+\sum_{i=\delta^{\prime}n+1}^{\delta
n}(1/2+\varepsilon)^{2}+\sum_{i=\delta n}^{n}\varepsilon^{2}}$
$\displaystyle=\frac{1}{\delta^{\prime}n}\left(\sqrt{(\delta^{\prime\prime}n)\varepsilon^{2}+(\delta^{\prime}-\delta^{\prime\prime})n\left(\frac{1}{2}+\varepsilon\right)^{2}+(\delta-\delta^{\prime})n\left(\frac{1}{2}+\varepsilon\right)^{2}+(1-\delta)n\varepsilon^{2}}\right)$
$\displaystyle=\frac{1}{\delta^{\prime}n}\left(\sqrt{(1-\delta-\delta^{\prime\prime})\varepsilon^{2}n+(\delta-\delta^{\prime\prime})\left(\frac{1}{2}+\varepsilon\right)^{2}n}\right)$
$\displaystyle=\frac{1}{\delta^{\prime}n}\left(\sqrt{(1-2\delta^{\prime\prime})\varepsilon^{2}n+\left(\frac{\delta-\delta^{\prime\prime}}{4}\right)n+\left(\frac{\delta-\delta^{\prime\prime}}{2}\right)\varepsilon
n}\right)$
$\displaystyle\leqslant\frac{1}{\delta^{\prime}n}\left(\sqrt{(1-2\delta^{\prime\prime})\varepsilon^{2}n+\frac{3n}{4}(\delta-\delta^{\prime\prime})}\right)$
$\displaystyle\leqslant\frac{1}{\delta^{\prime}n}\left(\sqrt{\varepsilon^{2}n+\frac{3n}{4}(\delta-\delta^{\prime\prime})}\right)\leqslant\sqrt{\frac{2\delta^{\prime\prime\prime}}{\delta^{\prime
2}}}\frac{1}{\sqrt{n}}\,,$

where we used the definition of the Frobenius norm, in conjunction with a
rewrite of the projector diagonalized into the eigenbasis of
$B^{L_{\alpha,\gamma}}_{d}(M)$ and the fact that $\varepsilon<1$.

∎

##  Appendix B Analytic Sparsity Proofs

###  B.1 Normalizing the Bernstein polynomial under the semi-circle measure

######  Lemma B.1 (Average Area of Bernstein Polynomial under $\mu_{sc}$).

For the semi-circle measure $d\mu_{sc}$ supported on $[-1,1]$ and
$B^{L^{\alpha}_{\gamma}}(x)$ as defined in Equation A.3 with a choice of
$L^{\alpha}_{\gamma}(x)$ given in Equation 6.14, the following holds,

$\operatorname*{\mathbb{E}}_{x\sim\mu_{sc}}\left[B^{L^{\alpha}_{\gamma}}(x)\right]\leqslant
1+\varepsilon-\mathsf{CDF}(1-\phi)\,,$

and,

$\operatorname*{\mathbb{E}}_{x\sim\mu_{sc}}\left[B^{L^{\alpha}_{\gamma}}(x)\right]\geqslant(1+2\varepsilon)\frac{4\sqrt{2}}{3\pi}\phi^{3/2}+O(\phi^{5/2})+o(\phi^{5/2})+\varepsilon\,.$

###### Proof.

We compute the area of the semi-circle where the Bernstein polynomial “lights
up” to normalize the correlation appropriately. To avoid unnecessary
complications, we will use the uniform lower bound of $0$ and uniform upper
bound of $\phi^{2}$ on the terms, leading us to a lower and upper bound of the
correct _leading_ order.

##### The upper bound

We being by computing the upper bound by using the uniform upper bound of
$\phi^{2}$ on the ramp part,

$\displaystyle\operatorname*{\mathbb{E}}_{x\sim\mu_{sc}}\left[B^{L^{\alpha}_{\gamma}}(x)\right]=\int_{-1}^{1}B^{L^{\alpha}_{\gamma}}(x)\frac{2}{\pi}\sqrt{1-x^{2}}dx\leqslant\int_{-1}^{1}\left(L^{\alpha}_{\gamma}(x)+\varepsilon\right)\frac{2}{\pi}\sqrt{1-x^{2}}dx$
$\displaystyle=\int_{-1}^{1}L^{\alpha}_{\gamma}(x)d\mu_{sc}(x)+\varepsilon\int_{-1}^{1}d\mu_{sc}(x)$
$\displaystyle=\int_{-1}^{1-\phi}0\cdot
d\mu_{sc}(x)+\frac{1}{\phi^{2}}\int_{1-\phi}^{1-\phi+\phi^{2}}\left(x-(1-\phi)\right)d\mu_{sc}(x)+\int_{1-\phi+\phi^{2}}^{1}d\mu_{sc}(x)+\varepsilon$
$\displaystyle\leqslant
0+\frac{1}{\phi^{2}}\phi^{2}\int_{1-\phi}^{1-\phi+\phi^{2}}d\mu_{sc}(x)+\int_{1-\phi+\phi^{2}}^{1}d\mu_{sc}(x)+\varepsilon$
$\displaystyle=\int_{1-\phi}^{1}d\mu_{sc}(x)+\varepsilon=(1+\varepsilon)-\mathsf{CDF}(1-\phi)\,.$

To evaluate the above, we take the Taylor expansion of $\mathsf{CDF}(1-\phi)$
around $1$ (using the fact that $\phi(\varepsilon)<<1$) and have,

$\displaystyle 1-\mathsf{CDF}(1-\phi)$
$\displaystyle=1-\left(\frac{1}{2}+\frac{(1-\phi)\sqrt{1-(1-\phi)^{2}}}{\pi}+\frac{\arcsin{1-\phi}}{\pi}\right)$
$\displaystyle\overset{\text{Taylor
expansion}}{=}1-\left(1-\frac{4\sqrt{2}}{3\pi}\phi^{3/2}+O(\phi^{5/2})+o(\phi^{5/2})\right)$
$\displaystyle=\frac{4\sqrt{2}}{3\pi}\phi^{3/2}+O(\phi^{5/2})+o(\phi^{5/2})\,,$

where $o(\phi^{5/2})$ denotes terms of higher power than $5/2$ and utilizes
the fact that $\phi(\varepsilon)<<1$.

##### The lower bound

The lower bound follows a similar computation, using the uniform lower bound
of $0$ on the ramp part,

$\displaystyle\operatorname*{\mathbb{E}}_{x\sim\mu_{sc}}\left[B^{L^{\alpha}_{\gamma}}(x)\right]=\int_{-1}^{1}B^{L^{\alpha}_{\gamma}}(x)\frac{2}{\pi}\sqrt{1-x^{2}}dx$
$\displaystyle\geqslant\int_{-1}^{1-\phi}\left(L^{\alpha}_{\gamma}+\varepsilon\right)\frac{2}{\pi}\sqrt{1-x^{2}}dx+\int_{1-\phi}^{1}\left(L^{\alpha}_{\gamma}-\varepsilon\right)\frac{2}{\pi}dx$
$\displaystyle=\int_{1-\phi}^{1-\phi+\phi^{2}}\frac{1}{\phi^{2}}\left(x-(1-\phi)\right)d\mu_{sc}(x)+\int_{1-\phi+\phi^{2}}^{1}d\mu_{sc}+\varepsilon\left(\int_{-1}^{1-\phi}d\mu_{sc}(x)-\int_{1-\phi}^{1}d\mu_{sc}(x)\right)$
$\displaystyle\geqslant
1-\mathsf{CDF}(1-\phi+\phi^{2})+\varepsilon\left(\mathsf{CDF}(1-\phi)-0-1+\mathsf{CDF}(1-\phi)\right)$
$\displaystyle=1-\mathsf{CDF}(1-\phi+\phi^{2})+2\varepsilon(\mathsf{CDF}(1-\phi)-\frac{1}{2})$
$\displaystyle=1-\left(\frac{1}{2}+\frac{(1-\phi(1-\phi))\sqrt{1-(1-\phi(1-\phi))^{2}}}{\pi}+\frac{\arcsin{1-\phi(1-\phi)}}{\pi}\right)$
$\displaystyle\qquad\qquad\qquad\qquad+2\varepsilon\left(\frac{1}{2}\frac{4\sqrt{2}}{3\pi}\phi^{3/2}+O(\phi^{5/2})+o(\phi^{5/2})\right)$
$\displaystyle=\frac{4\sqrt{2}}{3\pi}\phi^{3/2}(1-\phi)^{3/2}+O(\phi^{5/2}(1-\phi)^{5/2})+o(\phi^{5/2}(1-\phi^{5/2})+2\varepsilon\frac{4\sqrt{2}}{3\pi}\phi^{3/2}+\varepsilon$
$\displaystyle\geqslant(1+2\varepsilon)\frac{4\sqrt{2}}{3\pi}\phi^{3/2}+O(\phi^{5/2})+o(\phi^{5/2})+\varepsilon\,,$

where fourth equality uses a similar taylor expansion as the one used in
computing the upper bound. ∎

###  B.2 Auxiliary lemmata

We state and prove some helper lemmata, used to reason about matrix Bernstein
polynomials and low-degree SoS proofs.

######  Lemma B.2 (Products of univariate non-negative symmetric matrix
polynomials).

Given univariate matrix polynomials of symmetric matrices
$p_{1},p_{2}:M_{n}(\mathbb{R})\to\mathbb{R}$ with
$\mathsf{deg}(p_{1}),\mathsf{deg}(p_{2})\leqslant d$ and $0\preceq p_{1}(M)$,
$0\preceq p_{2}(M)$ for every $M\in M_{n}(\mathbb{R})$, the matrix polynomial
$p_{1}(M)p_{2}(M)$ is symmetric and non-negative in the Loewner order with
$\mathsf{deg}(p_{1}p_{2})\leqslant 2d$.

###### Proof.

By Proposition B.3, $p_{1}(M)p_{2}(M)$ has non-negative non-zero eigenvalues.
It remains to check that $p_{1}(M)p_{2}(M)$ is symmetric, as Proposition B.3
then also implies that $p_{1}(M)p_{2}(M)$ is a symmetric PSD matrix. This
follows by expressing the matrix polynomials explicitly over the monomial
basis with their real coefficients $\\{a_{i}\\}_{i=0}^{d}$ and
$\\{b_{i}\\}_{i=0}^{d}$.

$\displaystyle p_{1}(M)p_{2}(M)$
$\displaystyle=\left(\sum_{i=0}^{d}a_{i}M^{i}\right)\left(\sum_{j=0}^{d}b_{j}M^{j}\right)$
$\displaystyle=\sum_{i,j=0}^{d}a_{i}b_{j}M^{i}M^{j}=\sum_{i=0}^{d}\sum_{j=0}^{d}a_{i}b_{j}M^{i+j}$
$\displaystyle=\sum_{i=0}^{d}\sum_{j=0}^{d}a_{i}b_{j}(M^{\mathsf{T}})^{i}(M^{\mathsf{T}})^{j}=\sum_{i=0}^{d}\sum_{j=0}^{d}a_{i}b_{j}(M^{i})^{\mathsf{T}}(M^{j})^{\mathsf{T}}=\sum_{i=0}^{d}\sum_{j=0}^{d}a_{i}b_{j}(M^{j}M^{i})^{\mathsf{T}}$
$\displaystyle=\sum_{i=0}^{d}\sum_{j=0}^{d}a_{i}b_{j}(M^{i}M^{j})^{\mathsf{T}}=\left(\sum_{i=0}^{d}\sum_{j=0}^{d}a_{i}b_{j}M^{i}M^{j}\right)^{\mathsf{T}}=\left(p_{1}(M)p_{2}(M)\right)^{\mathsf{T}}\,.$

Note that $\mathsf{deg}(p_{1}p_{2})\leqslant 2d$ also follows from the
reasoning above. ∎

Stated below is an elementary proposition which shows that the products of
matrices that commute preserves the Lowener order.

######  Proposition B.3 (Non-Zero Eigenspectrum Equivalence of conjugated PSD
matrices).

Let $A$ and $B$ by symmetric PSD matrices in $M_{n}(\mathbb{R})$. Then,

$\displaystyle\lambda_{i}(AB)=\lambda_{i}(A^{1/2}BA^{1/2})\,,$ (B.1)

for all non-zero eigenvalues $\lambda_{i}(\cdot)$. Furthermore, if
$AB=(AB)^{\mathsf{T}}$, then $AB$ is a symmetric PSD matrix.

###### Proof.

Since $A$ is PSD,

$\displaystyle AB$ $\displaystyle=A^{1/2}A^{1/2}B\,.$

For every non-zero eigenvector $u_{i}$ of $A^{1/2}BA^{1/2}$ with eigenvalue
$\lambda_{i}$, the vector $v_{i}=A^{1/2}u_{i}$ is an eigenvector with
eigenvalue $u_{i}$ as,

$\displaystyle AB(A^{1/2}u_{i})=ABv_{i}$
$\displaystyle=A^{1/2}(A^{1/2}B)v_{i}=A^{1/2}(A^{1/2}BA^{1/2})u_{i}$
$\displaystyle=A^{1/2}\lambda_{i}u_{i}=\lambda_{i}(A^{1/2}u_{i})\,.$

For every $w\in\mathsf{ker}(A^{1/2})$ that is also a non-zero eigenvector of
$AB$, we can choose $w^{\prime}=A^{1/2}Bw$ as an eigenvector of
$A^{1/2}BA^{1/2}$ with the same eigenvalue. If there is some vector
$w\in\mathsf{ker}(A^{1/2})\cap\mathsf{ker}(B)$, then it is a zero eigenvector
of both $AB$ and $A^{1/2}A^{1/2}B$.
It remains to show that $A^{1/2}BA^{1/2}$ is PSD. This follows by the
elementary observation that, for every $v\in\mathbb{R}^{n}$,

$\displaystyle v^{\mathsf{T}}(A^{1/2}BA^{1/2})v=\left\lVert
B^{1/2}w\right\rVert^{2}_{2}\geqslant 0\,,$

where $w=A^{1/2}v$. This implies that the non-zero eigenvalues of $AB$ are
non-negative.
If, furthermore,

$\displaystyle AB=(AB)^{\mathsf{T}}=B^{\mathsf{T}}A^{\mathsf{T}}=BA\,,$

the spectral theorem for symmetric matrices immediately implies that $AB$ is a
symmetric PSD matrix. ∎

######  Lemma B.4 (Asymptotic evaluation of Catalan Numbers).

For all $q>0$, the Catalan number $C_{q}$ can be equivalently expressed as,

$C_{q}=\left(1+O\left(\frac{1}{q}\right)\right)\left(\frac{4^{q}}{\sqrt{\pi}q^{3/2}}\right)\,.$
(B.2)

###### Proof.

The proof follows by a simple application of the Stirling bounds for the
factorial function on the definition of the Catalan numbers.

$\displaystyle C_{q}$
$\displaystyle:=\binom{2q}{q}\frac{1}{q+1}=\frac{(2q)!}{(q!)^{2}(q+1)}$