$z\in{\cal S}^{N_{1}-1}$ implies the trace of the leading $N_{1}\times N_{1}$
block of $X$ is equal to 1. The construction above leads to a semidefinite
relaxation with size

$\displaystyle n=(N+1)N_{1},\quad
m=1+nN_{2}+\mathfrak{t}(N_{1}-1)\times\mathfrak{t}(N).$ (63)

Due to the limitation of interior point methods, Cifuentes Cifuentes19arXiv-
stls  was only able to numerically verify the tightness of the relaxation for
very small problems (_e.g.,_ $~{}N_{1}\leq N_{2}\leq 10,N<20$).

We aim to compute globally optimal solutions of (STLS) with much larger
dimensions. We perform experiments on random instances of (STLS) where the
affine map ${\cal L}$ is structured to be a square Hankel matrix such that
$N=N_{1}+N_{2}-1,N_{1}=N_{2}$. We set $N_{1}\in\\{10,20,30,40\\}$, and at each
level we randomly generate three problem instances by drawing $\theta\sim{\cal
N}(0,I_{N})$ from the standard Gaussian distribution. We then solve the sparse
semidefinite relaxation using SDPT3, MOSEK, CDCS, SketchyCGAL, SDPNAL+, and
STRIDE. For STRIDE, since the nonlinear programing (STLS) does not admit any
manifold structure, we use fmincon with an interior point method as the nlp
method. To generate hypotheses for nlp from the moment matrix $X$, we follow

$\displaystyle X=\sum_{i=1}^{n}\lambda_{i}v_{i}v_{i}^{\mathsf{T}},\quad
z_{i}=\frac{v_{i}[z]}{\left\lVert v_{i}[z]\right\rVert},\quad
u^{i}_{j}=z_{i}^{\mathsf{T}}v_{i}[u_{j}z],\;j=1,\dots,N,$ (64)

where we first perform a spectral decomposition, then round $z_{i}$ by
projecting the entries of $v_{i}$ corresponding to block $z$ onto the unit
sphere, and round $u^{i}_{j}$ by computing the inner product between $z_{i}$
and the entries of $v_{i}$ corresponding to block $u_{j}z$ (again, the
rationale for this rounding comes from the lifting (62)). We generate $r=5$
hypotheses by rounding $5$ eigenvectors. In order to set $M_{b}$ for computing
$\eta_{s}$, we make the assumption that the search variable $u$ of (STLS)
contains random vectors that follow ${\cal N}(0,I_{N})$, and its squared norm
follows a chi-square distribution of degree $N$. As a result, we choose the
quantile corresponding to a $99.9\%$ probability, denoted as $\overline{M}$,
as the bound on $\left\lVert u\right\rVert^{2}$, such that
$M_{b}=\overline{M}+1$ can upper bound the trace of $X$.

Table 4 gives the numerical results of different solvers. Notice that at
$N_{1}=40$, we have increased the maximum runtime of CDCS and SDPNAL+ to
$20000$ seconds to make a fair comparison with STRIDE. Generally, Table 4 has
similar results as Table 3. (i) IPMs can solve small problems to high
accuracy, but cannot handle large problems due to memory issues. (ii) First-
order solvers (CDCS and SketchyCGAL) are memory efficient but exhibit slow
convergence and cannot attain high accuracy to certify global optimality and
tightness of the relaxation. SDPNAL+ can solve problems with $N_{1}=10$ and
$20$ to high accuracy, but failed for $N_{1}=30$ and $40$. The degraded
performance of CDCS, SketchyCGAL, and SDPNAL+ compared to Tables 1-2 again
suggests the difficulty in solving sparse relaxations compared to dense
relaxations. (iii) STRIDE computed solutions of high accuracy for all test
instances and it is about $5$ times faster than SDPNAL+.

Table 4: Results on solving sparse second-order relaxation of random (STLS)
instances. “$**$” indicates solver out of memory.

Dimension Run Metric SDPT3 Toh99OMS-sdpt3  MOSEK mosek  CDCS Zheng20MP-CDCS
SketchyCGAL Yurtsever21SIMDS-scalableSDP  SDPNAL+ Yang15MPC-sdpnalplus  STRIDE
$\begin{subarray}{c}N_{1}:\ 10\\\ \\\ n:\ 200\\\ \\\ m:\ 10,551\end{subarray}$
$\\#1$ $\eta_{p}$ $4.8\mathrm{e}{-12}$ $1.6\mathrm{e}{-11}$
$3.7\mathrm{e}{-9}$ $0.3602$ $1.0\mathrm{e}{-10}$ $6.2\mathrm{e}{-16}$
$\eta_{d}$ $1.2\mathrm{e}{-10}$ $2.2\mathrm{e}{-15}$ $1.5\mathrm{e}{-5}$
$4.3\mathrm{e}{-14}$ $1.2\mathrm{e}{-11}$ $6.3\mathrm{e}{-15}$ $\eta_{g}$
$5.5\mathrm{e}{-6}$ $1.1\mathrm{e}{-12}$ $9.4\mathrm{e}{-9}$
$9.1\mathrm{e}{-4}$ $5.3\mathrm{e}{-9}$ $4.8\mathrm{e}{-12}$ $\eta_{s}$
$5.7\mathrm{e}{-8}$ $8.7\mathrm{e}{-14}$ $0.0078$ $0.8808$ $1.2\mathrm{e}{-8}$
$1.1\mathrm{e}{-11}$ time $94.0$ $48.2$ $138.6$ $339.7$ $13.3$ $14.1$ $\\#2$
$\eta_{p}$ $1.2\mathrm{e}{-10}$ $1.8\mathrm{e}{-10}$ $4.5\mathrm{e}{-6}$
$0.3787$ $1.7\mathrm{e}{-11}$ $1.5\mathrm{e}{-15}$ $\eta_{d}$
$1.7\mathrm{e}{-10}$ $4.1\mathrm{e}{-11}$ $3.8\mathrm{e}{-4}$
$2.6\mathrm{e}{-5}$ $1.0\mathrm{e}{-10}$ $2.5\mathrm{e}{-14}$ $\eta_{g}$
$5.5\mathrm{e}{-5}$ $5.2\mathrm{e}{-9}$ $1.9\mathrm{e}{-4}$ $0.2894$
$6.3\mathrm{e}{-10}$ $2.1\mathrm{e}{-11}$ $\eta_{s}$ $6.4\mathrm{e}{-6}$
$2.4\mathrm{e}{-9}$ $0.1179$ $0.8544$ $1.5\mathrm{e}{-8}$ $2.5\mathrm{e}{-11}$
time $92.2$ $43.9$ $140.9$ $344.2$ $48.6$ $18.4$ $\\#3$ $\eta_{p}$
$1.6\mathrm{e}{-10}$ $9.0\mathrm{e}{-12}$ $8.3\mathrm{e}{-7}$ $0.3943$
$7.5\mathrm{e}{-11}$ $8.2\mathrm{e}{-16}$ $\eta_{d}$ $2.3\mathrm{e}{-11}$
$9.7\mathrm{e}{-9}$ $1.4\mathrm{e}{-4}$ $1.6\mathrm{e}{-5}$
$8.9\mathrm{e}{-11}$ $1.0\mathrm{e}{-13}$ $\eta_{g}$ $1.9\mathrm{e}{-6}$
$1.7\mathrm{e}{-10}$ $5.3\mathrm{e}{-5}$ $0.0164$ $7.3\mathrm{e}{-11}$
$8.0\mathrm{e}{-11}$ $\eta_{s}$ $2.4\mathrm{e}{-8}$ $1.5\mathrm{e}{-11}$
$0.1406$ $0.8993$ $3.9\mathrm{e}{-13}$ $1.9\mathrm{e}{-10}$ time $99.2$ $44.4$
$139.8$ $338.1$ $38.9$ $17.4$ $\begin{subarray}{c}N_{1}:\ 20\\\ \\\ n:\ 800\\\
\\\ m:\ 164,201\end{subarray}$ $\\#1$ $\eta_{p}$ ** ** $2.7\mathrm{e}{-6}$
$0.9608$ $3.2\mathrm{e}{-13}$ $1.9\mathrm{e}{-15}$ $\eta_{d}$
$3.4\mathrm{e}{-4}$ $0.0325$ $7.8\mathrm{e}{-10}$ $3.8\mathrm{e}{-14}$
$\eta_{g}$ $2.2\mathrm{e}{-4}$ $0.1207$ $2.8\mathrm{e}{-7}$
$1.8\mathrm{e}{-10}$ $\eta_{s}$ $0.3945$ $0.9863$ $1.2\mathrm{e}{-8}$
$2.8\mathrm{e}{-10}$ time $4206$ $2171$ $1189$ $292.5$ $\\#2$ $\eta_{p}$ ** **
$2.5\mathrm{e}{-6}$ $0.9326$ $5.6\mathrm{e}{-13}$ $1.5\mathrm{e}{-15}$
$\eta_{d}$ $2.7\mathrm{e}{-4}$ $0.0109$ $1.7\mathrm{e}{-10}$
$2.4\mathrm{e}{-15}$ $\eta_{g}$ $1.6\mathrm{e}{-4}$ $0.4961$
$2.0\mathrm{e}{-10}$ $2.8\mathrm{e}{-10}$ $\eta_{s}$ $0.4805$ $0.9714$
$9.1\mathrm{e}{-8}$ $7.2\mathrm{e}{-10}$ time $4187$ $2163$ $870.6$ $257.7$
$\\#3$ $\eta_{p}$ ** ** $4.3\mathrm{e}{-6}$ $0.9450$ $5.3\mathrm{e}{-11}$
$1.6\mathrm{e}{-15}$ $\eta_{d}$ $3.4\mathrm{e}{-4}$ $0.0200$
$7.9\mathrm{e}{-11}$ $2.2\mathrm{e}{-15}$ $\eta_{g}$ $3.1\mathrm{e}{-4}$
$0.1536$ $1.3\mathrm{e}{-8}$ $1.2\mathrm{e}{-10}$ $\eta_{s}$ $0.4493$ $0.9798$
$4.1\mathrm{e}{-9}$ $2.5\mathrm{e}{-10}$ time $4224$ $2178$ $1169$ $324.4$
$\begin{subarray}{c}N_{1}:\ 30\\\ \\\ n:\ 1800\\\ \\\ m:\
823,951\end{subarray}$ $\\#1$ $\eta_{p}$ ** ** $1.4\mathrm{e}{-5}$ $1.3862$
$1.7\mathrm{e}{-12}$ $5.4\mathrm{e}{-15}$ $\eta_{d}$ $1.0\mathrm{e}{-3}$
$1.3752$ $2.4\mathrm{e}{-5}$ $4.9\mathrm{e}{-14}$ $\eta_{g}$
$1.7\mathrm{e}{-4}$ $0.5175$ $0.1060$ $1.6\mathrm{e}{-9}$ $\eta_{s}$ $0.87$
$0.9997$ $0.2083$ $6.8\mathrm{e}{-10}$ time $12234$ $6047$ $10000$ $1679$
$\\#2$ $\eta_{p}$ ** ** $1.6\mathrm{e}{-5}$ $1.3826$ $2.0\mathrm{e}{-14}$
$1.6\mathrm{e}{-8}$ $\eta_{d}$ $0.0011$ $1.3544$ $1.2\mathrm{e}{-5}$
$2.3\mathrm{e}{-11}$ $\eta_{g}$ $3.3\mathrm{e}{-4}$ $0.4171$ $0.1077$
$2.8\mathrm{e}{-7}$ $\eta_{s}$ $0.8248$ $0.9996$ $0.0454$ $2.7\mathrm{e}{-7}$
time $12324$ $5994$ $10001$ $1813$ $\\#3$ $\eta_{p}$ ** ** $1.4\mathrm{e}{-5}$
$1.3819$ $2.2\mathrm{e}{-12}$ $4.2\mathrm{e}{-15}$ $\eta_{d}$ $0.0010$
$1.3376$ $5.8\mathrm{e}{-6}$ $2.3\mathrm{e}{-15}$ $\eta_{g}$
$3.7\mathrm{e}{-4}$ $0.3730$ $0.0569$ $6.9\mathrm{e}{-10}$ $\eta_{s}$ $0.8191$
$0.9996$ $0.0640$ $4.6\mathrm{e}{-10}$ time $12360$ $5976$ $10000$ $2050$
$\begin{subarray}{c}N_{1}:\ 40\\\ \\\ n:\ 3200\\\ \\\ m:\
2,592,801\end{subarray}$ $\\#1$ $\eta_{p}$ ** ** ** $1.6887$
$9.7\mathrm{e}{-14}$ $5.6\mathrm{e}{-15}$ $\eta_{d}$ $2.8283$
$2.3\mathrm{e}{-5}$ $2.2\mathrm{e}{-12}$ $\eta_{g}$ $0.4648$ $0.0439$
$4.3\mathrm{e}{-10}$ $\eta_{s}$ $1.0$ $0.2345$ $1.3\mathrm{e}{-7}$ time
$17222$ $20003$ $14541$ $\\#2$ $\eta_{p}$ ** ** ** $1.6865$
$2.3\mathrm{e}{-5}$ $4.9\mathrm{e}{-15}$ $\eta_{d}$ $2.8528$
$2.6\mathrm{e}{-5}$ $1.0\mathrm{e}{-13}$ $\eta_{g}$ $0.4391$ $0.0750$
$2.1\mathrm{e}{-9}$ $\eta_{s}$ $1.0$ $0.3296$ $3.5\mathrm{e}{-9}$ time $17270$
$20008$ $9994$ $\\#3$ $\eta_{p}$ ** ** ** $1.6875$ $3.5\mathrm{e}{-4}$
$7.9\mathrm{e}{-15}$ $\eta_{d}$ $2.9918$ $1.4\mathrm{e}{-4}$
$5.3\mathrm{e}{-12}$ $\eta_{g}$ $0.5329$ $0.1743$ $2.0\mathrm{e}{-8}$
$\eta_{s}$ $1.0$ $0.6361$ $2.8\mathrm{e}{-7}$ time $17227$ $20345$ $10120$

##  6 Concluding Remarks

In this paper, we have designed an efficient algorithmic framework, STRIDE, to
solve large-scale tight semidefinite programming (SDP) relaxations of
polynomial optimization problems (POPs). STRIDE employs an inexact projected
gradient method (iPGM) as its backbone, but leverages fast nonlinear
programming (NLP) methods to seek rapid primal acceleration. By conducting a
novel convergence analysis for the iPGM and taking safeguarded steps, we have
proved the global convergence of STRIDE. For solving the projection subproblem
in every iPGM step, we have designed a modified limited-memory BFGS method and
proved its global convergence. In addition to the contributions in algorithmic
design and convergence analysis, we have studied several important classes of
POPs and their dense or sparse (tight) relaxations. In our extensive numerical
experiments, STRIDE globally solved all the POPs and the corresponding SDP
relaxations to very high accuracy, and it offered state-of-the-art efficiency
and robustness. We hope the practical performance of STRIDE can encourage more
theoretical study towards investigating when and why the NLP iterates can
produce effective acceleration to the iPGM backbone.

###### Acknowledgements.

The authors would like to thank Diego Cifuentes for useful discussions about
the nearest structured rank deficient matrix problem, and Jie Wang for
clarifications about sparse semidefinite relaxations. H. Yang and L. Carlone
were partially funded by ARL DCIST CRA W911NF-17-2-0181 and NSF CAREER award
“Certifiable Perception for Autonomous Cyber-Physical Systems”.

##  Appendix A Proof of Theorem 2.1

Let us first prove the following useful lemma for later convenience.

###### Lemma 1

Under the inexactness conditions in (17), for any $X\in\mathbb{S}_{+}^{n}$, it
holds that

$\left\langle C,X-X^{k}\right\rangle\geq\frac{1}{\sigma_{k}}\left\lVert
X^{k}-X^{k-1}\right\rVert^{2}+\frac{1}{\sigma_{k}}\left\langle
X^{k-1}-X,X^{k}-X^{k-1}\right\rangle+\left\langle\mathcal{A}^{*}y^{k},X-X^{k}\right\rangle.$

###### Proof

From the last two conditions in (17),
$\mathcal{A}^{*}y^{k}+S^{k}-C=\frac{1}{\sigma_{k}}\left(X^{k}-X^{k-1}\right)$
and $\left\langle X^{k},S^{k}\right\rangle=0$, we have that

$\displaystyle\left\langle C,X-X^{k}\right\rangle=$
$\displaystyle\;\left\langle\mathcal{A}^{*}y^{k}+S^{k}-\frac{1}{\sigma_{k}}\left(X^{k}-X^{k-1}\right),X-X^{k}\right\rangle$
$\displaystyle=$ $\displaystyle\;\frac{1}{\sigma_{k}}\left\langle
X^{k}-X^{k-1},X^{k}-X\right\rangle+\left\langle\mathcal{A}^{*}y^{k},X-X^{k}\right\rangle+\underbrace{\left\langle
S^{k},X\right\rangle}_{\geq 0}-\underbrace{\left\langle
S^{k},X^{k}\right\rangle}_{=0}$ $\displaystyle\geq$
$\displaystyle\;\frac{1}{\sigma_{k}}\left\lVert
X^{k}-X^{k-1}\right\rVert^{2}+\frac{1}{\sigma_{k}}\left\langle
X^{k-1}-X,X^{k}-X^{k-1}\right\rangle+\left\langle\mathcal{A}^{*}y^{k},X-X^{k}\right\rangle,$

where the last inequality is due to the following fact:

$\displaystyle\frac{1}{\sigma_{k}}\left\langle
X^{k}-X^{k-1},X^{k}-X\right\rangle=$
$\displaystyle\displaystyle\frac{1}{\sigma_{k}}\left\langle
X^{k}-X^{k-1},X^{k}-X^{k-1}+X^{k-1}-X\right\rangle$ $\displaystyle=$
$\displaystyle\displaystyle\frac{1}{\sigma_{k}}\left\lVert
X^{k}-X^{k-1}\right\rVert^{2}+\frac{1}{\sigma_{k}}\left\langle
X^{k-1}-X,X^{k}-X^{k-1}\right\rangle.$

Thus, the proof is completed. ∎

Then, we shall conduct our proof of Theorem 2.1 as follows.

###### Proof

For notational simplicity, denote $r^{k}=\mathcal{A}(X^{k})-b$ and
$e^{k}=X^{k}-X^{\star}$ for all $k\geq 0$. Then, $\mathcal{A}(e^{k})=r^{k}$
since $\mathcal{A}(X^{\star})=b$. Moreover, we have $\left\lVert
r^{k}\right\rVert\leq\varepsilon_{k}$ from the first condition in (17).
Therefore, the convergence of $\left\\{\left\lVert r^{k}\right\rVert\right\\}$
is obvious, since $\left\lVert r^{k}\right\rVert\leq\varepsilon_{k}$ and
$\\{k\varepsilon_{k}\\}$ is summable. Recall that $X^{\star}$ is an optimal
solution to problem (P). Let $X=X^{k-1}$ in Lemma 1, we get

$\left\langle C,X^{k-1}-X^{k}\right\rangle\geq\frac{1}{\sigma_{k}}\left\lVert
e^{k}-e^{k-1}\right\rVert^{2}+\left\langle\mathcal{A}^{*}y^{k},e^{k-1}-e^{k}\right\rangle$
(65)

On the other hand, let $X=X^{\star}$ in Lemma 1, we obtain

$\left\langle
C,X^{\star}-X^{k}\right\rangle\geq\frac{1}{\sigma_{k}}\left\lVert
e^{k}-e^{k-1}\right\rVert^{2}+\frac{1}{\sigma_{k}}\left\langle
e^{k-1},e^{k}-e^{k-1}\right\rangle-\left\langle
A^{*}y^{k},e^{k}\right\rangle.$ (66)

Multiplying $k-1$ to (65) and add the resulting inequality to (66) yields

$\displaystyle\;\left\langle
C,(k-1)\left(X^{k-1}-X^{\star}\right)-k\left(X^{k}-X^{\star}\right)\right\rangle$
(67) $\displaystyle\geq$ $\displaystyle\;\frac{k}{\sigma_{k}}\left\lVert
e^{k}-e^{k-1}\right\rVert^{2}+\frac{1}{\sigma_{k}}\left\langle
e^{k-1},e^{k}-e^{k-1}\right\rangle+\left\langle\mathcal{A}^{*}y^{k},(k-1)e^{k-1}-ke^{k}\right\rangle$
$\displaystyle=$ $\displaystyle\;\frac{1}{\sigma_{k}}\left(k\left\lVert
e^{k}\right\rVert^{2}+(k-1)\left\lVert
e^{k-1}\right\rVert^{2}-(2k-1)\left\langle
e^{k},e^{k-1}\right\rangle\right)+\left\langle
y^{k},(k-1)r^{k-1}-kr^{k}\right\rangle$ $\displaystyle\geq$
$\displaystyle\;\frac{1}{\sigma_{k}}\left(k\left\lVert
e^{k}\right\rVert^{2}+(k-1)\left\lVert
e^{k-1}\right\rVert^{2}-\frac{2k-1}{2}\left(\left\lVert
e^{k}\right\rVert^{2}+\left\lVert
e^{k-1}\right\rVert^{2}\right)\right)+\left\langle
y^{k},(k-1)r^{k-1}-kr^{k}\right\rangle$ $\displaystyle=$
$\displaystyle\;\frac{1}{2\sigma_{k}}\left\lVert
e^{k}\right\rVert^{2}-\frac{1}{2\sigma_{k}}\left\lVert
e^{k-1}\right\rVert^{2}+\left\langle y^{k},(k-1)r^{k-1}-kr^{k}\right\rangle$
$\displaystyle\geq$ $\displaystyle\;\frac{1}{2\sigma_{k}}\left\lVert
e^{k}\right\rVert^{2}-\frac{1}{2\sigma_{k-1}}\left\lVert
e^{k-1}\right\rVert^{2}+\left\langle y^{k},(k-1)r^{k-1}-kr^{k}\right\rangle,$

where in the last inequality, we use the fact that
$\sigma_{k}\geq\sigma_{k-1}>0$ for all $k\geq 1$. By summing the inequality in
(67) for $k=1,\dots,k$, we get (recall that $\left\lVert y^{k}\right\rVert\leq
M$ for all $k\geq 0$)

$\displaystyle-k\left\langle C,X^{k}-X^{\star}\right\rangle\geq$
$\displaystyle\;\frac{1}{2\sigma_{k}}\left\lVert
e^{k}\right\rVert^{2}-\frac{1}{2\sigma_{0}}\left\lVert
e_{0}\right\rVert^{2}+\sum_{i=1}^{k}\left\langle
y^{i},(i-1)r^{i-1}-ir^{i}\right\rangle$ $\displaystyle\geq$
$\displaystyle\;\frac{1}{2\sigma_{k}}\left\lVert
e^{k}\right\rVert^{2}-\frac{1}{2\sigma_{0}}\left\lVert
e_{0}\right\rVert^{2}-2M\sum_{i=1}^{k}i\varepsilon_{i}.$

The above inequality leads to the following upper bound on $\left\langle
C,X^{k}-X^{\star}\right\rangle$:

$\left\langle
C,X^{k}-X^{\star}\right\rangle\leq\frac{1}{k}\left(\frac{1}{2\sigma_{0}}\left\lVert
e^{0}\right\rVert^{2}+2M\sum_{i=1}^{k}i\varepsilon_{i}\right).$ (68)

Next, we shall provide a lower bound for $\left\langle
C,X^{k}-X^{\star}\right\rangle$ for $k\geq 1$. To this end, let us define the
set $\mathcal{F}_{k}$, which is an enlargement of the primal feasible set
$\mathcal{F}_{\mathrm{P}}$ with respect to $\varepsilon_{k}\geq 0$, as
follows:

$\mathcal{F}_{k}:=\left\\{X\in\mathbb{S}^{n}\;\middle|\;\left\lVert\mathcal{A}(X)-b\right\rVert\leq\varepsilon_{k},\;X\in\mathbb{S}_{+}^{n}\right\\},\quad\forall\,k\geq
1.$

Moreover, denote

$X^{\star}_{k}:=\operatorname*{arg\,min}\left\\{\left\langle
C,X\right\rangle\;\middle|\;X\in\mathcal{F}_{k}\right\\},\quad\forall\,k\geq
1.$

Since $(y^{\star},S^{\star})\in\mathbb{R}^{m}\times\mathbb{S}_{+}^{n}$ is a
dual optimal solution, by (Jiang2012, , Lemma 3.4), we have that

$0\leq\left\langle C,X^{\star}-X^{\star}_{k}\right\rangle\leq\left\lVert
y^{\star}\right\rVert\varepsilon_{k},\quad\forall\,k\geq 1.$

As a consequence, because $X^{k}\in{\cal F}_{k}$ and $\left\langle
C,X^{k}\right\rangle\geq\left\langle C,X^{\star}_{k}\right\rangle$, it holds
that

$\left\langle C,X^{k}-X^{\star}\right\rangle\geq\left\langle
C,X^{\star}_{k}-X^{\star}\right\rangle\geq-\left\lVert
y^{\star}\right\rVert\varepsilon_{k}=-\frac{\left\lVert y^{\star}\right\rVert
k\varepsilon_{k}}{k}.$ (69)

Combining (68) and (69), we see that

$-\frac{\left\lVert y^{\star}\right\rVert k\varepsilon_{k}}{k}\leq\left\langle
C,X^{k}-X^{\star}\right\rangle\leq\frac{1}{k}\left(\frac{1}{2\sigma_{0}}\left\lVert
e^{0}\right\rVert^{2}+2M\sum_{i=1}^{k}i\varepsilon_{i}\right).$

Finally, let us estimate the dual infeasibility. To this end, we get from (65)
that

$\frac{1}{\sigma_{k}}\left\lVert X^{k}-X^{k-1}\right\rVert^{2}\leq\left\langle
C,X^{k-1}-X^{\star}\right\rangle-\left\langle
C,X^{k}-X^{\star}\right\rangle+\left\lVert y^{k}\right\rVert\left(\left\lVert
r^{k-1}\right\rVert+\left\lVert r^{k}\right\rVert\right)$

which implies that

$\displaystyle\;\left\lVert\mathcal{A}^{*}y^{k}+S^{k}-C\right\rVert=\frac{1}{\sigma_{k}}\left\lVert
X^{k}-X^{k-1}\right\rVert$ $\displaystyle\leq$
$\displaystyle\;\frac{1}{\sqrt{\sigma_{k}}}\sqrt{\frac{1}{k-1}\left(\frac{1}{2\sigma_{0}}\left\lVert
e^{0}\right\rVert^{2}+2M\sum_{i=1}^{k-1}i\varepsilon_{i}\right)+\left\lVert
y^{\star}\right\rVert\varepsilon_{k}+M\left(\varepsilon_{k-1}+\varepsilon_{k}\right)}$
$\displaystyle\leq$
$\displaystyle\;O\left(\frac{1}{\sqrt{k\sigma_{k}}}\right),$

and establishes the convergence of the dual infeasibility. Therefore, the
proof is completed. ∎

##  Appendix B Proof for Theorem 4.3

From Step 1 of the modified L-BFGS algorithm, we have
$d^{k}=(-\beta_{k}I+Q_{k})\nabla\phi(\xi^{k})$, with $Q_{k}=0$ if $\left\lVert
d^{k}\right\rVert\geq K$. This leads to
$\left\langle\nabla\phi(\xi^{k}),d^{k}\right\rangle=\left\langle\nabla\phi(\xi^{k}),-(\beta_{k}I+Q_{k})\nabla\phi(\xi^{k})\right\rangle$,
which implies that

$\lambda_{\min}(\beta_{k}I+Q_{k})\leq\frac{\left\langle-\nabla\phi(\xi^{k}),d^{k}\right\rangle}{\left\lVert\nabla\phi(\xi^{k})\right\rVert^{2}},$
(70)

where $\lambda_{\min}(\cdot)$ denotes the minimum eigenvalue. This implies
that $d^{k}$ is a descent direction and the line search scheme is well-
defined.

Let the sequence $\\{\xi^{k}\\}$ be generated by the algorithm. We next show
the second part of the theorem, _i.e.,_ $\\{\xi^{k}\\}$ is a bounded sequence.
Since $\mathcal{A}\mathcal{A}^{*}$ is nonsigular and the Slater condition
holds, by (Rockafellar1974, , Theorem 17’ & 18’), the level set
$\mathcal{L}_{\phi}:=\left\\{\xi\in\mathbb{R}^{m}\,:\,\phi(\xi)\leq\phi(\xi^{0})\right\\}$
is compact. Thus, the boundedness of the sequence is proven.

Let us now prove the third part of the theorem, which states that every
accumulation point of $\\{\xi^{k}\\}$ is an optimal solution solution of
problem (39). To this end, without loss of generality, we assume that
$\\{\xi^{k}\\}$ is an infinite sequence (the case for a finite sequence is
trivial and we omit it here) and $\xi^{k}\rightarrow\hat{\xi}$ as
$k\rightarrow\infty$, and we shall prove that $\nabla\phi(\hat{\xi})=0$. We
plan to prove the latter result by contraction. Let us assume that
$\nabla\phi(\hat{\xi})\neq 0$. Under this assumption, we will deduce a
contradiction in three steps.

Step A.

We first claim that when $\nabla\phi(\hat{\xi})\neq 0$, $\left\lVert
d^{k}\right\rVert$ is bounded from below and above, _i.e.,_ there exist two
positive constants $M_{1}$ and $M_{2}$ such that $0<M_{1}\leq\left\lVert
d^{k}\right\rVert\leq M_{2}$ for all $k$.

On the one hand, by Step 1 of the L-BFGS Algorithm 5, we know that either
$d^{k}=-(\beta_{k}I+Q_{k})\nabla\phi(\xi^{k})$, in which case $\left\lVert
d^{k}\right\rVert\leq K$ (otherwise the algorithm rejects it), or
$d^{k}=-\beta_{k}\nabla\phi(\xi^{k})$, in which case $\left\lVert
d^{k}\right\rVert\leq\tau_{1}M_{0}^{1+\tau_{2}}$ (due to
$\beta_{k}=\tau_{1}\left\lVert\nabla\phi(\xi^{k})\right\rVert^{\tau_{2}}$).
Therefore, we have $\left\lVert
d^{k}\right\rVert\leq\max\left\\{K,\tau_{1}M_{0}^{1+\tau_{2}}\right\\}:=M_{2}<\infty$,
and $\left\lVert d^{k}\right\rVert$ is upper bounded.

On the other hand, if for some index set $J$, $\left\lVert
d^{k_{j}}\right\rVert\rightarrow 0$ for $j\in J$, then we can deduce that

$\left\lVert\nabla\phi(\xi^{k_{j}})\right\rVert=\left\lVert(\beta_{k_{j}}I+Q_{k_{j}})^{-1}d^{k_{j}}\right\rVert\leq\frac{1}{\beta_{k_{j}}}\left\lVert
d^{k_{j}}\right\rVert=\frac{1}{\tau_{1}\left\lVert\nabla\phi(\xi^{k_{j}}\right\rVert}\left\lVert
d^{k_{j}}\right\rVert,\quad j\in J.$

The above inequality implies that
$\left\lVert\nabla\phi(\hat{\xi})\right\rVert=\lim_{j\in
J}\left\lVert\nabla\phi(\xi^{k_{j}})\right\rVert=0$, a contradiction. Thus,
$\left\lVert d^{k}\right\rVert$ must also be lower bounded, by some constant
$M_{1}>0$.

Step B.

We next claim that $\\{\rho^{m_{k}}\\}$ is bounded away from zero.

Suppose that, for some index set $J$, $\rho^{m_{k_{j}}}\rightarrow 0$, $j\in
J$ (_i.e.,_ $m_{k_{j}}\rightarrow\infty$ due to $\rho\in(0,1)$). Recall that
$m_{k_{j}}$ is the smallest nonnegative integer to satisfy the inequality in
the line search scheme (so that $m_{k_{j}}-1$ does not satisfy the
inequality). Thus, it holds that

$\phi(\xi^{k_{j}}+\rho^{m_{k_{j}}-1}d^{k_{j}})>\phi(\xi^{k_{j}})+\mu\rho^{m_{k_{j}}-1}\left\langle\nabla\phi(\xi^{k_{j}}),d^{k_{j}}\right\rangle.$
(71)

Meanwhile, we can expand $\phi(\xi^{k_{j}}+\rho^{m_{k_{j}}-1}d^{k_{j}})$ as

$\phi(\xi^{k_{j}}+\rho^{m_{k_{j}}-1}d^{k_{j}})=\phi(\xi^{k_{j}})+\int_{0}^{1}\left\langle\nabla\phi(\xi^{k_{j}}+t\rho^{m_{k_{j}}-1}d^{k_{j}}),\rho^{m_{k_{j}}-1}d^{k_{j}}\right\rangle
dt.$ (72)

By substituting (72) into (71), we obtain the following inequality:

$\int_{0}^{1}\left\langle\nabla\phi(\xi^{k_{j}}+t\rho^{m_{k_{j}}-1}d^{k_{j}}),d^{k_{j}}\right\rangle
dt>\mu\left\langle\nabla\phi(\xi^{k_{j}}),d^{k_{j}}\right\rangle,\quad j\in
J.$ (73)

Since the sequence $\\{d^{k_{j}}\\}$ is bounded from both below and above, by
taking a subset of $J$ if necessary, we assume that there exists $\hat{d}$
such that $d^{k_{j}}\rightarrow\hat{d}$, $j\in J$. Taking limit with respect
to $j\in J$ in (73) yields that
$\left\langle\nabla\phi(\hat{\xi}),\hat{d}\right\rangle\geq\mu\left\langle\nabla\phi(\hat{\xi}),\hat{d}\right\rangle$,
which further implies that
$\left\langle\nabla\phi(\hat{\xi}),\hat{d}\right\rangle\geq 0$ by the fact
that $\mu\in(0,1/2)$. Also, since
$\left\langle\nabla\phi(\hat{\xi}),\hat{d}\right\rangle\leq 0$, we derive
$\left\langle\nabla\phi(\hat{\xi}),\hat{d}\right\rangle=0$. However, denoting
$\hat{\beta}$ as the limit of $\beta_{k_{j}}$ as $k_{j}\rightarrow\infty$,
$j\in J$, and noting that $Q_{k}\succeq 0$ for all $k$, the inequality in (70)
shows that

$0=-\left\langle\nabla\phi(\hat{\xi}),\hat{d}\right\rangle\geq\hat{\beta}\left\lVert\nabla\phi(\hat{\xi})\right\rVert^{2}\geq\tau_{1}\left\lVert\nabla\phi(\hat{\xi})\right\rVert^{2+\tau_{2}},$

which implies that $\nabla\phi(\hat{\xi})=0$. This is also a contradiction.
Therefore, $\\{\rho^{m_{k}}\\}$ is bounded away from zero.

Step C.

Now, using the fact that $\phi$ is bounded from below on the bounded sequence
$\\{\xi^{k}\\}$, we have that $\\{\phi(\xi^{k+1})-\phi(\xi^{k})\\}\rightarrow
0$ (since the sequence $\\{\phi(\xi^{k})\\}$ is nonincreasing and
$\phi(\cdot)$ is continuous). This implies, by the line search rule again,
that
$\rho^{m_{k}}\left\langle\nabla\phi(\xi^{k}),d^{k}\right\rangle\rightarrow 0$,
as $k\rightarrow\infty$. Therefore, it holds that
$\left\langle\nabla\phi(\hat{\xi}),\hat{d}\right\rangle=0$, since
$\\{\rho^{m_{k}}\\}$ is bounded away from zero, which again implies that
$\nabla\phi(\hat{\xi})=0$, and this contradicts our initial hypothesis that
$\nabla\phi(\hat{\xi})\neq 0$.

At this point, we have shown that $\nabla\phi(\hat{\xi})=0$, as desired. At
last, by the convexity of $\phi(\cdot)$, vanishing of the gradient at
$\hat{\xi}$ implies that $\hat{\xi}$ is an optimal solution of problem (39).
Therefore, the proof is completed. ∎

## References

  * (1) Abbe, E., Bandeira, A.S., Hall, G.: Exact recovery in the stochastic block model.  IEEE Transactions on Information Theory 62(1), 471–487 (2015)
  * (2) Absil, P.A., Mahony, R., Sepulchre, R.: Optimization algorithms on matrix manifolds.  Princeton University Press (2009)
  * (3) Aholt, C., Agarwal, S., Thomas, R.: A QCQP approach to triangulation.  In: European Conference on Computer Vision, pp. 654–667. Springer (2012)
  * (4) Alizadeh, F., Haeberly, J.P.A., Overton, M.L.: Complementarity and nondegeneracy in semidefinite programming.  Mathematical programming 77(1), 111–128 (1997)
  * (5) Alizadeh, F., Haeberly, J.P.A., Overton, M.L.: Primal-dual interior-point methods for semidefinite programming: convergence rates, stability and numerical results.  SIAM Journal on Optimization 8(3), 746–768 (1998)
  * (6) Antonante, P., Tzoumas, V., Yang, H., Carlone, L.: Outlier-robust estimation: Hardness, minimally-tuned algorithms, and applications.  IEEE Trans. Robotics (2021)
  * (7) ApS, M.: The MOSEK optimization toolbox for MATLAB manual. Version 9.0. (2019).  URL http://docs.mosek.com/9.0/toolbox/index.html
  * (8) Barak, B., Brandao, F.G., Harrow, A.W., Kelner, J., Steurer, D., Zhou, Y.: Hypercontractivity, sum-of-squares proofs, and their applications.  In: Proceedings of the forty-fourth annual ACM symposium on Theory of computing, pp. 307–326 (2012)
  * (9) Bay, H., Tuytelaars, T., Van Gool, L.: Surf: Speeded up robust features.  In: European conference on computer vision, pp. 404–417. Springer (2006)
  * (10) Beck, A., Teboulle, M.: A fast iterative shrinkage-thresholding algorithm for linear inverse problems.  SIAM journal on imaging sciences 2(1), 183–202 (2009)
  * (11) Bertsekas, D.: Nonlinear Programming.  Athena Scientific (1999)
  * (12) Blekherman, G., Parrilo, P.A., Thomas, R.R.: Semidefinite optimization and convex algebraic geometry.  SIAM (2012)
  * (13) Boumal, N., Mishra, B., Absil, P.A., Sepulchre, R.: Manopt, a Matlab toolbox for optimization on manifolds.  Journal of Machine Learning Research 15(42), 1455–1459 (2014).  URL https://www.manopt.org
  * (14) Boumal, N., Voroninski, V., Bandeira, A.S.: The non-convex Burer-Monteiro approach works on smooth semidefinite programs.  In: Conference on Neural Information Processing Systems (NeurIPS) (2016)
  * (15) Briales, J., Gonzalez-Jimenez, J.: Convex global 3d registration with lagrangian duality.  In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 4960–4969 (2017)
  * (16) Burer, S., Monteiro, R.D.: A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization.  Mathematical Programming 95(2), 329–357 (2003)
  * (17) Burer, S., Ye, Y.: Exact semidefinite formulations for a class of (random and non-random) nonconvex quadratic programs.  Mathematical Programming pp. 1–17 (2019)
  * (18) Candes, E.J., Eldar, Y.C., Strohmer, T., Voroninski, V.: Phase retrieval via matrix completion.  SIAM review 57(2), 225–251 (2015)
  * (19) Chen, L., Sun, D., Toh, K.: An efficient inexact symmetric gauss-seidel based majorized ADMM for high-dimensional convex composite conic programming.  Mathematical Programming 161, 237–270 (2017)
  * (20) Cifuentes, D.: A convex relaxation to compute the nearest structured rank deficient matrix.  SIAM Journal on Matrix Analysis and Applications 42(2), 708–729 (2021)
  * (21) Cifuentes, D., Agarwal, S., Parrilo, P.A., Thomas, R.R.: On the local stability of semidefinite relaxations.  Mathematical Programming (2021)
  * (22) Combettes, P.L., Pesquet, J.C.: Proximal splitting methods in signal processing.  In: Fixed-point algorithms for inverse problems in science and engineering, pp. 185–212. Springer (2011)
  * (23) Dai, Y.H., Yuan, Y.: A nonlinear conjugate gradient method with a strong global convergence property.  SIAM Journal on optimization 10(1), 177–182 (1999)
  * (24) De Klerk, E.: The complexity of optimizing over a simplex, hypercube or sphere: a short survey.  Central European Journal of Operations Research 16(2), 111–125 (2008)
  * (25) Doherty, A.C., Parrilo, P.A., Spedalieri, F.M.: Complete family of separability criteria.  Physical Review A 69(2), 022308 (2004)
  * (26) Fang, K., Fawzi, H.: The sum-of-squares hierarchy on the sphere and applications in quantum information theory.  Mathematical Programming pp. 1–30 (2020)
  * (27) Fukuda, M., Kojima, M., Murota, K., Nakata, K.: Exploiting sparsity in semidefinite programming via matrix completion i: General framework.  SIAM Journal on optimization 11(3), 647–674 (2001)
  * (28) Gamarnik, D., Kızıldağ, E.C.: Algorithmic obstructions in the random number partitioning problem.  arXiv preprint arXiv:2103.01369 (2021)
  * (29) Goemans, M.X., Williamson, D.P.: Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming.  Journal of the ACM (JACM) 42(6), 1115–1145 (1995)
  * (30) Helmberg, C., Rendl, F., Weismantel, R.: A semidefinite programming approach to the quadratic knapsack problem.  Journal of combinatorial optimization 4(2), 197–215 (2000)
  * (31) Henrion, D., Lasserre, J.B.: GloptiPoly: Global optimization over polynomials with Matlab and SeDuMi.  ACM Transactions on Mathematical Software (TOMS) 29(2), 165–194 (2003)
  * (32) Inc., W.R.: Mathematica, Version 12.2.  URL https://www.wolfram.com/mathematica.  Champaign, IL, 2020
  * (33) Jiang, K., Sun, D., Toh, K.C.: An inexact accelerated proximal gradient method for large scale linearly constrained convex sdp.  SIAM Journal on Optimization 22(3), 1042–1064 (2012)
  * (34) Josz, C., Henrion, D.: Strong duality in Lasserre’s hierarchy for polynomial optimization.  Optimization Letters 10(1), 3–10 (2016)
  * (35) Kaltofen, E., Yang, Z., Zhi, L.: Approximate greatest common divisors of several polynomials with linearly constrained coefficients and singular polynomials.  In: Proceedings of the 2006 international symposium on Symbolic and algebraic computation, pp. 169–176 (2006)
  * (36) Lasserre, J.B.: An explicit exact SDP relaxation for nonlinear 0-1 programs.  In: International Conference on Integer Programming and Combinatorial Optimization, pp. 293–303. Springer (2001)
  * (37) Lasserre, J.B.: Global optimization with polynomials and the problem of moments.  SIAM Journal on Optimization 11(3), 796–817 (2001)
  * (38) Lasserre, J.B.: Moments, positive polynomials and their applications, vol. 1.  World Scientific (2009)
  * (39) Laurent, M.: Semidefinite representations for finite varieties.  Mathematical programming 109(1), 1–26 (2007)
  * (40) Li, X., Sun, D., Toh, K.C.: A block symmetric Gauss–Seidel decomposition theorem for convex composite quadratic programming and its applications.  Mathematical Programming 175(1), 395–418 (2019)
  * (41) Ling, C., Nie, J., Qi, L., Ye, Y.: Biquadratic optimization over unit spheres and semidefinite programming relaxations.  SIAM Journal on Optimization 20(3), 1286–1310 (2010)
  * (42) Luo, Z.Q., Ma, W.K., So, A.M.C., Ye, Y., Zhang, S.: Semidefinite relaxation of quadratic optimization problems.  IEEE Signal Processing Magazine 27(3), 20–34 (2010)
  * (43) Mai, N.H.A., Magron, V., Lasserre, J.B.: A hierarchy of spectral relaxations for polynomial optimization.  arXiv preprint arXiv:2007.09027 (2020)
  * (44) Malick, J., Povh, J., Rendl, F., Wiegele, A.: Regularization methods for semidefinite programming.  SIAM Journal on Optimization 20(1), 336–356 (2009)
  * (45) Malick, J., Sendov, H.S.: Clarke generalized Jacobian of the projection onto the cone of positive semidefinite matrices.  Set-Valued Analysis 14(3), 273–293 (2006)
  * (46) Markovsky, I.: Structured low-rank approximation and its applications.  Automatica 44(4), 891–909 (2008)
  * (47) Markovsky, I., Usevich, K.: Software for weighted structured low-rank approximation.  Journal of Computational and Applied Mathematics 256, 278–292 (2014)
  * (48) Markovsky, I., Van Huffel, S.: Overview of total least-squares methods.  Signal processing 87(10), 2283–2302 (2007)
  * (49) Meneghetti, G., Danelljan, M., Felsberg, M., Nordberg, K.: Image alignment for panorama stitching in sparsely structured environments.  In: Scandinavian Conference on Image Analysis, pp. 428–439. Springer (2015)
  * (50) Mertens, S.: Number partitioning.  Computational Complexity and Statistical Physics p. 125 (2006)
  * (51) Nesterov, Y.: Lectures on convex optimization, vol. 137.  Springer (2018)
  * (52) Nie, J.: Polynomial optimization with real varieties.  SIAM Journal On Optimization 23(3), 1634–1646 (2013)
  * (53) Nie, J.: Optimality conditions and finite convergence of Lasserre’s hierarchy.  Mathematical programming 146(1), 97–121 (2014)
  * (54) Nocedal, J., Wright, S.: Numerical Optimization.  Springer Science & Business Media (2006)
  * (55) Parrilo, P.A.: Semidefinite programming relaxations for semialgebraic problems.  Mathematical programming 96(2), 293–320 (2003)
  * (56) Rockafellar, R.T.: Conjugate duality and optimization.  SIAM (1974)
  * (57) Rosen, D.M.: Scalable low-rank semidefinite programming for certifiably correct machine perception.  In: Intl. Workshop on the Algorithmic Foundations of Robotics (WAFR), vol. 3 (2020)
  * (58) Rosen, D.M., Carlone, L., Bandeira, A.S., Leonard, J.J.: SE-Sync: A certifiably correct algorithm for synchronization over the special Euclidean group.  The International Journal of Robotics Research 38(2-3), 95–125 (2019)
  * (59) Rosen, J.B., Park, H., Glick, J.: Total least norm formulation and solution for structured problems.  SIAM Journal on matrix analysis and applications 17(1), 110–126 (1996)
  * (60) Rusu, R., Blodow, N., Beetz, M.: Fast point feature histograms (fpfh) for 3d registration.  In: IEEE Intl. Conf. on Robotics and Automation (ICRA), pp. 3212–3217. Citeseer (2009)
  * (61) Shi, J., Yang, H., Carlone, L.: Optimal pose and shape estimation for category-level 3d object perception.  In: Robotics: Science and Systems (RSS) (2021)
  * (62) Shi, J., Yang, H., Carlone, L.: ROBIN: a graph-theoretic approach to reject outliers in robust estimation using invariants.  In: IEEE Intl. Conf. on Robotics and Automation (ICRA) (2021)
  * (63) Shor, N.Z.: Dual quadratic estimates in polynomial and boolean programming.  Annals of Operations Research 25(1), 163–168 (1990)
  * (64) Sun, D., Toh, K.C., Yang, L.: A convergent 3-block semiproximal alternating direction method of multipliers for conic programming with 4-type constraints.  SIAM journal on Optimization 25(2), 882–915 (2015)
  * (65) Toh, K.C.: Solving large scale semidefinite programs via an iterative solver on the augmented systems.  SIAM Journal on Optimization 14(3), 670–698 (2004)
  * (66) Toh, K.C., Todd, M.J., Tütüncü, R.H.: SDPT3—a MATLAB software package for semidefinite programming, version 1.3.  Optimization methods and software 11(1-4), 545–581 (1999)
  * (67) Waki, H., Kim, S., Kojima, M., Muramatsu, M.: Sums of squares and semidefinite programming relaxations for polynomial optimization problems with structured sparsity.  SIAM J. Optimization 17, 218–242 (2006)
  * (68) Wang, A.L., Kılınç-Karzan, F.: On the tightness of SDP relaxations of QCQPs.  Mathematical Programming pp. 1–41 (2021)
  * (69) Wang, J., Magron, V., Lasserre, J.B.: Chordal-TSSOS: a moment-SOS hierarchy that exploits term sparsity with chordal extension.  SIAM Journal on Optimization 31(1), 114–141 (2021)
  * (70) Wang, J., Magron, V., Lasserre, J.B.: TSSOS: A Moment-SOS hierarchy that exploits term sparsity.  SIAM Journal on Optimization 31(1), 30–58 (2021)
  * (71) Wu, G., Sun, J., Chen, J.: Optimal linear quadratic regulator of switched systems.  IEEE transactions on automatic control 64(7), 2898–2904 (2018)
  * (72) Yang, H., Antonante, P., Tzoumas, V., Carlone, L.: Graduated non-convexity for robust spatial perception: From non-minimal solvers to global outlier rejection.  IEEE Robotics and Automation Letters (RA-L) 5(2), 1127–1134 (2020)
  * (73) Yang, H., Carlone, L.: A polynomial-time solution for robust registration with extreme outlier rates.  In: Robotics: Science and Systems (RSS) (2019)
  * (74) Yang, H., Carlone, L.: A quaternion-based certifiably optimal solution to the Wahba problem with outliers.  In: Proceedings of the IEEE/CVF International Conference on Computer Vision, pp. 1665–1674 (2019)
  * (75) Yang, H., Carlone, L.: In perfect shape: Certifiably optimal 3d shape reconstruction from 2d landmarks.  In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 621–630 (2020)
  * (76) Yang, H., Carlone, L.: One ring to rule them all: Certifiably robust geometric perception with outliers.  In: Conference on Neural Information Processing Systems (NeurIPS) (2020)
  * (77) Yang, H., Shi, J., Carlone, L.: Teaser: Fast and certifiable point cloud registration.  IEEE Transactions on Robotics (2020)
  * (78) Yang, L., Sun, D., Toh, K.C.: SDPNAL+: a majorized semismooth Newton-CG augmented Lagrangian method for semidefinite programming with nonnegative constraints.  Mathematical Programming Computation 7(3), 331–366 (2015)
  * (79) Yurtsever, A., Tropp, J.A., Fercoq, O., Udell, M., Cevher, V.: Scalable semidefinite programming.  SIAM Journal on Mathematics of Data Science 3(1), 171–200 (2021)
  * (80) Zeng, A., Song, S., Nießner, M., Fisher, M., Xiao, J., Funkhouser, T.: 3dmatch: Learning the matching of local 3d geometry in range scans.  In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, vol. 1, p. 4 (2017)
  * (81) Zhang, R.Y., Lavaei, J.: Sparse semidefinite programs with guaranteed near-linear time complexity via dualized clique tree conversion.  Mathematical programming 188(1), 351–393 (2021)
  * (82) Zhao, X.Y., Sun, D., Toh, K.C.: A Newton-CG augmented Lagrangian method for semidefinite programming.  SIAM Journal on Optimization 20(4), 1737–1765 (2010)
  * (83) Zheng, Y., Fantuzzi, G., Papachristodoulou, A., Goulart, P., Wynn, A.: Chordal decomposition in operator-splitting methods for sparse semidefinite programs.  Mathematical Programming 180(1), 489–532 (2020)


# Yield-stress transition in suspensions of deformable droplets

Giuseppe Negro1, Livio Nicola Carenza2, Giuseppe Gonnella1, Fraser Mackay3,
Alexander Morozov3, Davide Marenduzzo3 1Dipartimento di Fisica, Universitá
degli Studi di Bari and INFN, Sezione di Bari, via Amendola 173, Bari,
I-70126, Italy,
2 Instituut-Lorentz, Universiteit Leiden, P.O. Box 9506, 2300 RA Leiden,
Netherlands,
3 SUPA, School of Physics and Astronomy, University of Edinburgh, Peter
Guthrie Tait Road, Edinburgh, EH9 3FD, UK

###### Abstract

Yield-stress materials, which require a sufficiently large forcing to flow,
are currently ill-understood theoretically. To gain insight into their
yielding transition, here we study numerically the rheology of a suspension of
deformable droplets under pressure-driven flow. We show that the suspension
displays yield-stress behaviour, with the droplets remaining motionless when
the applied body-force is below a critical value. In the non-flowing phase,
droplets jam to form an amorphous structure, whereas they order in the flowing
phase. Yielding is linked to a percolation transition in the contacts of
droplet-droplet overlaps, and requires suitable wetting boundary conditions
and strict conservation of the droplet area to exist. Close to the yielding
transition, we find strong oscillations in the droplet motion which closely
resemble those found experimentally in confined colloidal glasses under flow.
We show that even when droplets are static the underlying solvent moves by
permeation, so that the viscosity of the composite system is never truly
infinite, and, as we discuss, its precise value ceases to be a bulk material
property of the system.

Yield-stress fluids are materials which flow only when subject to a
sufficiently large stress, or external forcing Jones2002 ; Bonn2017 . The
critical stress above which there is flow is known as the yield stress.
Examples of yield-stress fluids abound in everyday materials and include
toothpaste, whipping or shaving cream, mayonnaise and cement. The defining
property of an ideal yield-stress fluid is that the apparent viscosity should
be infinite below yielding, so that the yield stress should mark a transition
between a solid-like and a fluid-like regime. Nevertheless, in practice it is
often arduous to distinguish this behaviour from that of a strongly shear-
thinning fluid for which the viscosity drops by orders of magnitude at the
yielding point, such that the material always flows albeit very slowly under
any external forcing, however small Bonn2017 .

Phenomenological theories for yield-stress fluids typically assume a non-
Newtonian and non-linear relation between the shear stress $\sigma$ and the
shear rate (or velocity gradient) $\dot{\gamma}$. A popular model is the
Herschel-Bulkley one Herschel1926 , which is based on the generic equation
$\sigma=\sigma_{y}+\eta_{\rm\infty}{\dot{\gamma}}^{n}$, with $\sigma_{y}$ the
yield stress, $\eta_{\infty}$ a material parameter and $n$ a generic exponent
found by fitting experimental data, and smaller than $1$ for shear-thinning
fluids. Phenomenological models like this are extremely useful to analyse and
compare experiments, but – by their nature – they do not address the
fundamental physical mechanisms underlying the existence of a yield stress.

Yield-stress fluids can be characterised according the softness of their
constituents Bonn2017 ; Zac2 , and range from bubble foams Heller1987 ;
Zinchenko2017 ; Cloitre2003  to suspensions of nearly-hard colloidal spheres
Pusey1986 ; Pham2006  (e.g., spherical particles stabilised sterically with a
thin polymer layer). In all cases, at large enough particle concentrations –
such that the system is the jammed, or glassy, phase respectively – these
materials are experimentally known to undergo a yielding transition. They also
display soft glassy rheology, as described by the Herschel-Bulkley model
Sollich1997 ; Holmes2003 ; Paredes2013 . In colloidal fluids, rheological
experiments further show that the effective viscosity of the system becomes
very large, and possibly diverges Bonn2017 ; Pham2006 . A confounding factor
hampering a conclusive demonstration of ideal yield-stress behaviour in
experiments is that the solid-like phase in a colloidal glass is amorphous,
and the fundamental physics of the amorphous state is not fully understood
Bonn2017 . Additionally, as we show to be relevant here, colloidal glasses or
foams are composite materials, so that the behaviour of the dispersed
particles and the underlying solvent may differ, thereby complicating the
picture.

Here we consider a generic universal model system for a yield-stress fluid: a
suspension of soft deformable droplets embedded in a Newtonian fluid
Foglino2017 ; Foglino2018 ; Tiribocchi2020 ; Zac1 . By changing the particle
softness – in particular interpenetrability and surface tension – our model
encompasses foams, stabilised oil-in-water emulsions and colloidal suspensions
as special cases. We show that these deformable suspensions display the
hallmark of yield-stress fluids, as the droplets are immobile even when
subjected to a (small) pressure difference, or forcing. In this immobile
phase, the droplets are arranged in an amorphous pattern, and the network of
droplet-droplet contacts, or overlaps, percolates. These overlaps provide the
soft analogue of frictional contacts, which are known to play a crucial role
in colloidal rheology Wyart2014 . Upon yielding, contact percolation is lost,
while droplets order as they flow. Importantly, we find that even in the phase
where droplets are static the solvent flows by permeation, meaning that the
viscosity of the overall system is never truly infinite. Close to the yielding
point, sustained velocity oscillations occur, similarly to what found
experimentally in flowing colloidal fluids close to the glass transition
Isa2009 . Our results allow to gain more insight into the microscopic physical
mechanisms underlying the yielding transition. In our system, the latter is
controlled by an inverse Bingham number, measuring the ratio between viscous
and interfacial forces. Notably, this is the same number which controls
discontinuous shear thinning at larger forces Foglino2017 ; Foglino2018 ;
Fei2020 . Finally, our simulations suggest that yielding disappears when the
droplets area is not strictly conserved, suggesting that systems featuring
evaporation-condensation phenomena can evade yield-stress behaviour and flow
under any forcing.

Figure 1: Yielding transition in the conserved model. (a)-(b) Color map of
$\phi=\sum_{i}\phi_{i}$ for $f<f_{c}$ ($f=2.0\times 10^{-6}$ in panel (a)) and
$f>f_{c}$ ($f=4.0\times 10^{-6}$ in panel (b)), for the conserved model. Black
and red regions correspond to $\phi=0$ and $\phi=2$ respectively. (c) Free
energy of overlaps (see text) as a function of body-force $f$. The insets of
panel (c) show clusters of contacting droplets, resulting from a density-based
spatial clustering analysis on the free energy of overlaps SI . Different
colors correspond to different clusters. Left and right inset correspond to
the configuration shown in panels (a) and (b) respectively. Movies of the
dynamics corresponding to (a) and (b) can be seen in Suppl. Movies 1 and 2
respectively SI .

To study the rheology of our soft droplet suspension, we work in $2D$ (Fig. 1)
and consider two models: the first strictly conserves the area of each
droplet, the second allows it to fluctuate around a target value, for instance
due to evaporation or condensation phenomena. We refer to these as the
conserved and non-conserved model, respectively. In both cases, the $N$
droplets in the system are non-coalescing, and we ensure this by describing
them in terms of $N$ distinct phase fields, $\phi_{i}$, with $i=1,\ldots,N$.
The hydrodynamics of the suspension can then be studied by following the
coupled evolution of the phase-field variables and of the velocity field
$\mathbf{v}$ of the underlying solvent. The use of phase-field means that
lubrication forces, which are notoriously challenging to accurately account
for in simulations Nguyen2002 ; Wyart2014  are altogether absent. At the same
time, frictional forces, which are known to play a crucial role in colloidal
rheology Wyart2014 , are present, in the form of overlaps between different
phase-fields.

The thermodynamics of the conserved model is governed by a free energy
$\mathcal{F}$ whose density is

$\sum_{i}^{N}\frac{\alpha}{4}\phi_{i}^{2}(\phi_{i}-\phi_{0})^{2}+\frac{K}{2}\sum_{i}^{N}(\nabla\phi_{i})^{2}+\sum_{i,j,i<j}\epsilon\phi_{i}^{2}\phi_{j}^{2}.$
(1)

Here, the first two terms favour the formation of circular droplets with
$\phi_{i}\simeq\phi_{0}$ in their interior, and $\phi_{i}\simeq 0$ outside.
The material constants $\alpha$ and $K$ determine the surface tension
$\gamma=\sqrt{8K\alpha/9}$ and the interfacial thickness
$\xi=\sqrt{2K/\alpha}$ of the droplets Pagonabarraga2002 . The term
proportional to $\epsilon>0$ describes soft repulsion pushing droplets apart
when overlapping. The phase-field variables evolve according to a set of
coupled Cahn-Hilliard equations,

$\frac{\partial\phi_{i}}{\partial
t}+\mathbf{v}\cdot\mathbf{\nabla}\phi_{i}=M\nabla^{2}\mu_{i}$ (2)

where $M$ is the mobility, $\mu_{i}=\delta{\cal F}/\delta\phi_{i}$ the
chemical potential of the $i$-th droplet. The flow obeys the Navier-Stokes
equation

$\rho\Big{(}\frac{\partial}{\partial
t}+{\mathbf{v}}\cdot{\nabla}\Big{)}{\mathbf{v}}=-{\nabla p}+\mathbf{f}^{\rm
th}+\eta_{0}\nabla^{2}{\mathbf{v}}+f\mathbf{e_{y}},$ (3)

where $\rho$ indicates the total fluid density, $p$ denotes the hydrodynamic
pressure and $\eta_{0}$ the solvent viscosity SI . The term
$\mathbf{f}^{th}=-\sum_{i}\phi_{i}{\nabla}\mu_{i}$ stands for the internal
thermodynamic force field due to the presence of non-trivial compositional
order parameters, while $f$ is the magnitude of the body-force, which we take
along the horizontal direction (Fig. 1a).

In our second model for the non-conserved concentration field, the free energy
${\cal F}$ is supplemented by an additional term,

$\mathcal{F}_{\rm constraint}=\lambda\left(1-\frac{1}{\pi
R^{2}\phi_{0}^{2}}\int dydz\,\phi_{i}^{2}\right)^{2},$ (4)

with $\lambda>0$ a constant which quantifies droplet compressibility, and
provides a soft constraint for the droplet area. The phase-fields evolve
according to a relaxational and overdamped dynamics defined by

$\frac{\partial\phi_{i}}{\partial
t}+\mathbf{v}\cdot\nabla\phi_{i}=-\frac{1}{\Gamma}\frac{\delta\mathcal{F^{\prime}}}{\delta\phi_{i}},$
(5)

where $\Gamma$ is a friction-like parameter and
$\mathcal{F^{\prime}}=\mathcal{F}+\mathcal{F}_{constraint}$. The equation for
$\mathbf{v}$ is still given by Eq. (3).

The dynamics are integrated with a parallel hybrid lattice Boltzmann approach
Tiribocchi2009 ; Marenduzzo2007 ; Carenza2019  where Eq. (3) is solved by a
lattice Boltzmann algorithm, and Eqs. (2),(5) are solved by finite difference
methods. We consider flow in a channel with no-slip boundary conditions at the
top and bottom walls, driven by a fixed pressure difference along the $y$
direction – leading to Poiseuille, parabolic, flow for a Newtonian fluid. At
the wall, neutral wetting boundary conditions are imposed for each droplet.
For more details, and a full list of parameters used, see SI .

We first study the rheological response of a droplet suspension (with packing
fraction $\varphi\simeq 0.5$) in the conserved model. A key result is that
there exists a critical body-force $f_{c}$ separating two fundamentally
different behaviours. For small forcing (Figs. 1a and 2a) the suspended
droplets are jammed and settle into a stationary non-flowing configuration
(Suppl. Movie 1) where they are immobile for the whole duration of the
numerical experiment ($\sim\mathcal{O}(10^{8})$ iterations). The snapshot
shown in Fig. 1a shows a typical droplet configuration for this regime. For
larger $f$, there is a subtle morphological rearrangement of the droplets
(Fig. 1b), which is accompanied by a yielding transition, as droplets now
steadily move (Fig. 2a). The snapshot shown in Fig. 1b shows a typical late
time configuration, which is travelling from left to right at a fixed velocity
(Fig. 2a and Suppl. Movie 2). An inspection of the configurations shows that
while the non-flowing state is amorphous (Fig. 1a), in the flowing state
droplets order (Fig. 1b; see SI  for a quantification of flow-induced
ordering). This morphological adjustment is accompanied by a fundamental
change in the patterns of contacts, or overlaps, between droplets. As shown in
the left inset of Fig. 1c, such overlaps create a percolating network in the
non-flowing state, whereas after yielding contacts no longer percolate along
the flow gradient direction (right inset of Fig. 1c and Fig S1). The change in
droplet contacts can be quantified by plotting the overlap free energy
($\mathcal{F}_{\textrm{overlaps}}=\epsilon\int
dydz\sum_{i,j}\phi_{i}^{2}\phi_{j}^{2}$) as a function of body-force (Fig.
1c): this quantity drops sharply at the yielding transition, corresponding to
the loss of contacts between droplets near the wall (right inset). As
discussed in more detail below, another key feature is that droplets need to
deform at least transiently when the system yields SI .

To quantify the yielding transition we compute two quantities: _(i)_ the mean
velocity of the droplets’ center of mass $\langle v_{y}\rangle_{d},$ (Figs.
2a,c) and _(ii)_ the throughput flow $Q=\int dydz\ \textrm{v}_{y}$ (Figs.
2b,d). While $\langle v_{y}\rangle_{d}$ quantifies the motility of the
suspended particles, $Q$ can be used to compute the effective viscosity of the
suspension, $\eta_{\rm eff}$. The latter quantity can be estimated as
$\eta_{\rm eff}=\eta_{0}\frac{Q_{0}}{Q}$, where
$Q_{0}=\frac{fL^{3}}{12\eta_{0}}$ is the throughput flow of a Newtonian fluid
with viscosity $\eta_{0}$ subject to a body-force $f$, and $L$ the channel
width. The yield-stress behavior is apparent from the plot of $\langle
v_{y}\rangle_{d}$ in Fig. 2a. Close to criticality, the mean droplet speed
behaves as $\sim(f-f_{c})^{\beta}$, with $\beta\simeq 0.54$ (dotted curve in
the inset of Fig. 2a). The phenomenology resembles that of the Prandtl-
Thomlinson model (where $\beta=1/2$) which describes a particle in a dashboard
potential, and provides a simple microscopic model for dry friction Popov2017
.

Importantly, even in the non-flowing phase in which droplets are at rest, the
underlying solvent flows (Fig. 2b): indeed $Q$ is non-zero for all values of
$f$. In more detail we find that there is a well-defined linear regime at
small forcing, which corresponds to a high but finite effective viscosity
(inset of Fig. 2b). In stark contrast, yield-stress fluid under shear exhibit
wall slip and an infinite effective viscosity Bonn2017 . This shows the exact
value of $\eta_{\rm eff}$ depends on the geometry of the system, and hence can
no longer be viewed one of its bulk material property. The flow at $0<f<f_{c}$
is purely permeative, as the solvent flows through an immobile network of
jammed droplets. The distinct behaviour of the droplet and solvent components
in the suspension is instructive, and shows that the composite material
behaves in a more complex way than what would be predicted for an ideal
single-phase yield-stress fluid. In our conserved model, yielding can
therefore be viewed as a continuous transition between a permeation regime
with jammed amorphous droplets where solvent flows mainly around them, and a
flowing ordered phase. In the latter phase, the flow is plug-like SI ;
Foglino2017 , as found experimentally for colloidal suspensions in a pressure-
driven flow Isa2007 .

It is interesting to contrast the behaviour we have just discussed with that
of the non-conserved model, where evaporation and condensation effects are
included. Surprisingly, replacing strict area conservation with a soft
constraint leads to a complete loss of yield-stress behaviour (Figs. 2c,d). In
the non-conserved model, droplets flow at any value of the forcing, however
small, so that it is not possible to define a yield stress. While a yielding-
like behaviour can still be observed as a smooth crossover, there is no longer
a singularity in the droplet velocity curve (Fig. 2c). An analysis of the area
of each droplet show that the droplet motion is accompanied by area
oscillations whose magnitude in controlled by $\lambda$, signalling that
motion occurs via evaporation-condensation (inset of Fig. 2c). The behaviour
of the throughput flow mirrors that of the droplet velocity in this non-
conserved model (Fig. 2d)

Figure 2: Flow behaviour in the conserved and non-conserved models. (a)-(b)
Average droplet velocity (a) and throughput flow (b) for the conserved model.
The inset of panel (a) shows the mean droplet speed close to criticality and
the result of the fit (dashed line) with the function $\langle
v\rangle\propto(f-f_{c})^{\beta}$, with $\beta\simeq 0.54$. The inset of panel
(b) shows the effective viscosity $\eta_{\textrm{eff}}$ as a function of the
body-force $f$. (c)-(d) Average droplet velocity (c) and throughput flow (d)
for the non-conserved model. The inset of panel (c) shows the area of three
nearby droplets versus time for $f=1.0\times 10^{-6}$.

More insight into the fundamental difference between the conserved and non-
conserved models can be gained by analysing the behaviour of a single droplet
at a solid wall under an external forcing, and with neutral wetting boundary
conditions (Fig. S5 and Suppl. Movie 3). While in the conserved model the
droplet sticks to the wall and requires a finite forcing to start moving, in
the non-conserved model evaporation and condensation provide another pathway
for contact line motion Bonn2009 , and the droplet drifts along the wall for
any value of the forcing. Therefore, besides the presence of a percolating
network of droplet overlaps (Fig. 1), the existence of a well-defined yielding
transition also requires a suitable behaviour of droplets close to the wall.

To understand the microscopic mechanism underlying yielding in the conserved
model more deeply, we now consider their dynamics close to $f_{c}$. Just after
yielding, we find a “stick-slip” behaviour where the emulsion alternates
between plug-like motion, where droplets flow, with stationary spells, where
they are almost jammed (Suppl. Movie 4). The throughput solvent flow and the
average droplet velocity both show irregular oscillations over time (red and
orange curves in Fig. 3a). The average variance (or amplitude) of the
stochastic oscillations increases with the forcing, and approaches zero at
$f_{c}$ (Fig. 3b). This behaviour is reminiscent of that found in velocity
oscillations of colloidal glasses close to the yielding transition Isa2009 .
There are indeed some key qualitative analogies between the two cases. In both
systems, the non-flowing and flowing states subtly differ in the typical
particle configuration. In our non-flowing emulsions, overlaps between
droplets abound and create a nearly percolating chain through the system, just
like frictional contacts for hard-sphere colloids. Instead, in the flowing
states there are gaps between most particles Isa2009 . Analysing the dynamics
in more details reveals an important distinction, though. In our system,
instantaneous yielding events – i.e., transitions from jammed to flowing
states – are typically accompanied by a sudden change in behaviour in the
deformation free energy
($\mathcal{F}_{\textrm{def}}=\frac{K}{2}\int\sum_{i}\left(\nabla\phi_{i}\right)^{2}dydz$,
see SI ). This suggest that yielding in our deformable suspensions requires a
transient change in droplet shape, which is instead essentially fixed for
colloids.

Figure 3: Oscillations near the yielding transition. (a) Throughput flow
versus time, for the conserved model, for different values of $f$ near
$f_{c}=3.15\times 10^{-6}$. (b) Plot of the variance of the oscillations as a
function of $f$.

To verify that our qualitative mechanism for yielding through interfacial
deformations is correct, we independently varied the parameters in Eq. 1 to
see how they affect the value of the critical forcing. We found that $f_{c}$
scales linearly with surface tension, $\gamma$ (Fig. 4a), and interfacial
width, $\xi$ (see Fig. S6). The only other parameters appreciably affecting
$f_{c}$ are the system size $L$ and the droplet radius $R$: increasing either
of these lengthscales leads to a decrease in $f_{c}$ (Fig. 4b). Our data
therefore suggest that a key dimensionless parameter may be the capillary
number $Ca=fLR^{2}/(\gamma\xi)$, which was also empirically found to determine
the physics of discontinuous shear thinning Foglino2017 ; Foglino2018 ;
Fei2020 . This can be viewed as an inverse Bingham number
$\sigma_{v}/\sigma_{y}$, with $\sigma_{v}\sim fL$ the viscous stress and
$\sigma_{y}\sim\gamma\xi/R^{2}$ an effective yield stress. The form of this
dimensionless control parameter suggests that in order for the suspension to
yield the external forcing has to overcome free energy barriers associated
with changes in particle shape, whose cost increases with $\gamma$ and $\xi$.

In summary, we studied the rheology of a soft droplet suspension under
pressure-driven flow. We found that the droplets only start moving when the
forcing they are subjected to exceeds a critical threshold, as in an ideal
yield-stress fluid. However, unlike one such material, even when droplets are
jammed the solvent flows through them via permeation, as in sheared
cholesteric cholestericpermeation  and smectic liquid crystals
elasticpermeation , leading to an effective viscosity which depends on the
system geometry and ceases to be a bulk property of the material. Yielding is
accompanied by a morphological transition. The jammed phase is amorphous and
the network of droplet-droplet contacts, or overlaps, percolates in the
direction perpendicular to the wall, conferring rigidity to the system. In the
flowing phase, droplets order and contact percolation is lost. Within this
picture, overlaps play a qualitatively similar role to frictional contacts in
hard colloids Wyart2014 . In our case, though, the transition between the
jammed and flowing phase requires a transient change in droplet shape. More
quantitatively, yielding occurs for a sufficiently large value of an inverse
Bingham number, controlling the balance between viscous and interfacial
stresses. The mechanism is therefore similar to that determining discontinuous
shear thinning at larger forcing Foglino2017 ; Foglino2018 ; Fei2020 : the
main difference is that at the yielding transition interfacial deformations
are spatially localised and transient in time, whereas at the discontinuous
shear thinning transitions they affect large portions of the system and occur
at all times.

Figure 4: Yielding phase diagram. (a) Phase diagram as a function of body-
force $f$ and surface tension $\gamma$. Orange squares: flowing systems;
purple circles: non-moving states. (b) Phase diagram as a function of $f$ and
system size $L$ plane. Red squares: flowing systems; blue circles: non-moving
states.

Strikingly, we predict the yield-stress behaviour can be completely eliminated
in our model by allowing droplet areas to fluctuate, for instance due to
evaporation/condensation phenomena. We hope that our results will stimulate
experiments to directly test our predictions, such as the importance of
permeation and the scaling of $f_{c}$. To assess the universality of our
results, one could investigate the yielding transition in other materials,
such as biological tissues Bi2016 ; Chiang2016 ; Loewe2020 , red blood cell
suspensions Lazaro2014 , and liquid crystalline emulsions lcemulsions .

##  I Acknowledgments

###### Acknowledgements.

The work has been performed under the Project HPC-EUROPA3
(INFRAIA-2016-1-730897), with the support of the EC Research Innovation Action
under the H2020 Programme. Part of this work was carried out on the Dutch
national e-infrastructure with the support of SURF through the Grant 2021.028
for computational time (L.N.C and G.N.). We acknowledge funding from MIUR
Project No. PRIN 2020/PFCXPE.

## References

  * (1) R. A. L. Jones, Soft Condensed Matter (Oxford University Press, New York, 2002).
  * (2) D. Bonn, M. M. Denn, L. Berthier, T. Divous and S. Manneville, Yield stress materials in soft condensed matter, Rev. Mod. Phys. 89, 035002 (2017).
  * (3) W. Herschel and R. Bulkley, Konsistenzmessungen von Gummi-Benzollösungen, Kolloid Z. 39, 291 (1926).
  * (4) J. Ruiz-Franco, F. Camerin, N. Gnan, and E. Zaccarelli, Phys. Rev. Materials 4, 045601 (2020).
  * (5) J. P. Heller and M. S. Kuntamukkula, Critical review of the foam rheology literature, Ind. Eng. Chem. Res. 26, 318-325 (1987).
  * (6) A. Z. Zinchenko and R. H. Davis, General rheology of highly concentrated emulsions with insoluble surfactants, J. Fluid Mech. 816, 661 (2017).
  * (7) M. Cloitre, R. Borrega, F. Monti and L. Leibler, Glassy dynamics and flow properties of soft colloidal pastes, Phys. Rev. Lett. 90, 068303 (2003).
  * (8) P. Pusey and W. van Megen, Phase behaviour of concentrated suspensions of nearly hard colloidal spheres, Nature 320, 340 (1986).
  * (9) K. N. Pham, G. Petekidis, D. Vlassopoulos, S. U. Egelhaaf, P. N. Pusey and W. C. K. Poon, Yielding of colloidal glasses, Europhys. Lett. 75, 624 (2006).
  * (10) P. Sollich, F. Lequeux, P. Hebraud and M. E. Cates, Rheology of soft glassy materials, Phys. Rev. Lett. 78, 2020 (1997).
  * (11) C. B. Holmes, M. Fuchs and M. E. Cates, Jamming transitions in a schematic model of suspension rheology, Europhys. Lett. 63, 240 (2003).
  * (12) J. Paredes, M. A. J. Michels and D. Bonn, Rheology across the zero-temperature jamming transition, Phys. Rev. Lett. 111, 105701 (2013).
  * (13) M. Foglino, A. N. Morozov, O. Henrich and D. Marenduzzo, Flow of deformable droplets: discontinuous shear thinning and velocity oscillations, Phys. Rev. Lett. 119, 208002 (2017).
  * (14) M. Foglino, A. N. Morozov and D. Marenduzzo, Rheology and microrheology of deformable droplet suspensions, Soft Matter 14, 9361 (2018).
  * (15) A. Tiribocchi, A. Montessori, F. Bonaccorso, M. Lauricella and S. Succi, Concentrated phase emulsion with multicore morphology under shear: a numerical study, Phys. Rev. Fluids 5, 113606 (2020).
  * (16) N. Gnan and E. Zaccarelli, The microscopic role of deformation in the dynamics of soft colloids, Nat. Phys. 15, 683–688 (2019).
  * (17) M. Wyart and M. E. Cates, Discontinous shear thickening without inertia in dense non-Brownian suspensions, Phys. Rev. Lett. 112, 098302 (2014).
  * (18) L. Isa, R. Besseling, A. N. Morozov and W. C. K. Poon, Velocity oscillations in microfluidic flows of concentrated colloidal suspensions, Phys. Rev. Lett. 102, 058302 (2009).
  * (19) L. Fei, A. Scagliarini, K. H. Luo and S. Succi, Discrete fluidization of dense monodisperse emulsions in neutral wetting microchannels, Soft Matter, 16, 651 (2020).
  * (20) N.-Q. Nguyen and A. J. C. Ladd, Lubrication corrections for lattice-Boltzmann simulations of particle suspensions, Phys. Rev. E 66, 046702 (2002).
  * (21) I. Pagonabarraga, A. J. Wagner and M. E. Cates, Bindary fluid demixing: the crossover region, J. Stat. Phys. 107, 39 (2002).
  * (22) M. E. Cates and E. Tjhung, Theories of binary fluid mixtures: from phase-separation kinetics to active emulsions, J. Fluid Mech. 836, P1 (2017).
  * (23) A. Tiribocchi, N. Stella, G. Gonnella and A. Lamura, Hybrid lattice Boltzmann model for binary fluid mixtures, Phys. Rev. E 80, 026701 (2009).
  * (24) D. Marenduzzo, E. Orlandini, M. E. Cates and J. M. Yeomans, Steady-state hydrodynamic instabilities of active liquid crystals: hybrid lattice Boltzmann simulations, Phys. Rev. E 76, 031921 (2007).
  * (25) L. N. Carenza, G. Gonnella, A. Lamura, G. Negro and A. Tiribocchi, Lattice Boltzmann methods and active fluids, Eur. Phys. J. E 42, 81 (2019).
  * (26) See online Supplemental Material at XXX, which includes more simulation details, additional results and Refs. DBSCAN ; carenza2020lattice ; carenza2020lattice1 .
  * (27) V. L. Popov, Contact Mechanics and Friction. Springer, Berlin, Heidelberg (2017).
  * (28) L. Isa, R. Besseling and W. C. K. Poon, Shear zones and wall slip in the capillary flow of concentrated colloidal suspensions, Phys. Rev. Lett. 98, 198305 (2007).
  * (29) D. Bonn, J. Eggers, J. Indekeu, J. Meunier and E. Rolley, Wetting and spreading, Rev. Mod. Phys. 81, 739 (2009).
  * (30) W. Helfrich, Capillary flow of cholesteric and smectic liquid crystals, Phys. Rev. Lett. 23, 372 (1969).
  * (31) J. Prost, Y. Pomeau and E. Guyon, Stability of permeative flows in 1 dimensionally ordered systems, J. Phys. II 1, 289 (1991).
  * (32) D. Bi, X. Yang, M. C. Marchetti and M. L. Manning, Motility-driven glass and jamming transitions in biological tissues, Phys. Rev. X 6, 021011 (2016).
  * (33) M. Chiang and D. Marenduzzo, Glass transitions in the cellular Potts model, EPL 116, 28009 (2016).
  * (34) B. Loewe, M. Chiang, D. Marenduzzo and M. C. Marchetti, Solid-liquid transition of deformable and overlapping active particles, Phys. Rev. Lett. 125, 038003 (2020).
  * (35) G. R. Lázaro, A. Hernández-Machado and I. Pagonabarraga, Rheology of red blood cells under flow in highly confined microchannels: I. effect of elasticity, Soft Matter 10, 7195 (2014).
  * (36) J. S. Lintuvuori, K. Stratford, M. E. Cates and D. Marenduzzo, Mixtures of Blue Phase Liquid Crystal with Simple Liquids: Elastic Emulsions, Phys. Rev. Lett. 121, 037802 (2018).
  * (37) M. Ester, H.-P. Kriegel, J. Sander, and X. Xu, A density-based algorithm for discovering clusters in large spatial databases with noise, Proceedings of the Second International Conference on Knowledge Discovery and Data Mining (1996).
  * (38) L. N. Carenza, G. Gonnella, D. Marenduzzo, and G. Negro, Rotation and propulsion in 3D active chiral droplets, Proc. Natl. Acad. Sci. U.S.A. 116, 22065 (2019).
  * (39) L. N. Carenza, G. Gonnella, D. Marenduzzo, and G. Negro, Chaotic and periodical dynamics of active chiral droplets, Physica (Amsterdam) 559A, 125025 (2020).


# Ad Hoc Teamwork in the Presence of Adversaries

Ted Fujimoto  Samrat Chatterjee  Auroop Ganguly

###### Abstract

Advances in ad hoc teamwork have the potential to create agents that
collaborate robustly in real-world applications. Agents deployed in the real
world, however, are vulnerable to adversaries with the intent to subvert them.
There has been little research in ad hoc teamwork that assumes the presence of
adversaries. We explain the importance of extending ad hoc teamwork to include
the presence of adversaries and clarify why this problem is difficult. We then
propose some directions for new research opportunities in ad hoc teamwork that
leads to more robust multi-agent cyber-physical infrastructure systems.

Machine Learning, ICML

##  1 Introduction

Ad hoc teamwork (AHT) is a step towards general multi-agent intelligence where
decision-making agents (human and computational) are trained to collaborate
with other previously unseen agents. That is, AHT assumes that each agent will
learn to collaborate without prior coordination or direct control of the other
agents (Mirsky et al., 2022). Although the challenge of AHT was first
formalized by Stone et al. (2010), there has been a recent resurgence of AHT
research that incorporates advances in deep reinforcement learning (RL).

Established works have focused on multi-agent consensus among a network of
agents typically with centralized training and decentralized execution, where
a fixed set of agents are known to all (Lowe et al., 2017; Rashid et al.,
2018). For example, this may be useful in real-world infrastructure protection
settings where a team of collaborative defense agents might be tasked with
securing an asset of interest. This might not be practical, however, in
dynamic, adversarial, and contested environments. In these settings, new
agents may need to enter the environment and start competing or collaborating
without access to a centralized policy. Even if we assume such a policy, it
would likely lead to brittle agents that are tethered to their training
partners for optimal performance (Li et al., 2019). This leads the group to be
unable to perform optimally in the presence of adversaries. One appeal of AHT
is the potential for multiple people or organizations to train agents that
collaborate with each other without the need to retrain with previously unseen
agents. Outside of Shafipour & Fallah (2021), AHT in the presence of
adversaries has not been explored due to the assumption of full cooperation
between agents. Adversarial machine learning (ML) researchers have an
opportunity to chart the next path of truly robust multi-agent learning. In
this paper, we contribute (1) a summary of past research in AHT, (2) why it
could be a necessary next step in high-stakes, real-world multi-agent RL, (3)
a list of some challenges in AHT in the presence of adversaries, and (4) a
proposal for some new directions for researchers in both adversarial ML and
AHT.

##  2 Ad Hoc Teamwork

###  2.1 Background

AHT is a subfield of cooperative AI, which is defined as AI methods that help
individuals, humans and machines, to find ways that improve joint welfare
(Dafoe et al., 2020). As mentioned previously, AHT methods train agents to
collaborate without prior coordination or direct control of other agents. In
this paper, AHT will mostly be synonymous with RL-based methods in AHT, called
zero-shot coordination (ZSC). The main reason is because ZSC has enjoyed some
recent success. These successes include a general and differentiable
information theoretic objective for training a diverse population of optimal
policies (Lupu et al., 2021), training methods that prevent agents from
exchanging information through arbitrary conventions (Hu et al., 2021), and
graph neural networks that learn agent models and joint-action value models
under varying team compositions (Rahman et al., 2021). Some approaches even
outperform self-play and behavioral cloning when collaborating with humans in
video games (Strouse et al., 2021). ZSC is an evolving field and there is
currently some disagreement over how to precisely formulate it (Treutlein et
al., 2021). We recommend Mirsky et al. (2022) as an up-to-date survey of AHT.

###  2.2 Potential Applications

While AHT has shown some success in robot soccer (Genter et al., 2017), there
are still many potential applications yet to be explored. Ideally, AHT can
perform more robustly than centralized multi-agent RL methods, so any multi-
agent environment can potentially train an AHT application. In the context of
secure cyber-physical infrastructure operations (such as transportation and
energy distribution), multi-agent learning and control is essential for
achieving system resilience goals (Januário et al., 2019; Hou et al., 2021;
Phan et al., 2021). Some future AHT applications worth mentioning include
important automated infrastructure, like energy and power distribution
(Biagioni et al., 2021; Pigott et al., 2021; Wang et al., 2021). AHT would be
particularly useful here because of the possibility that different buildings
or power plants might need to coordinate with previously unseen agents from
other allied power plants. Another application would be human-automated
vehicle interaction. It is still difficult for automated cars to learn and
appropriately react to humans while driving (Dommès et al., 2021).

##  3 Potential Challenges with Adversaries in Ad Hoc Teamwork

We now briefly describe some challenges that agents trained using AHT methods
could face unless it is framed as an adversarial ML problem.

###  3.1 Past Research in Adversarial RL

It is possible that some attacks described in past research on adversarial RL
(Huang et al., 2017; Gleave et al., 2019; Lin et al., 2020; Sun et al., 2020;
Figura et al., 2021; Fujimoto et al., 2021) will have similar negative effects
in AHT. Given that there are currently no universally effective defenses
against adversarial attacks in supervised learning (Athalye et al., 2018;
Short et al., 2019), the same is likely true for RL.

###  3.2 Balancing the Presence of Allies and Adversaries

Even if we assume adversaries might exist in the environment, it is not clear
how to train the agents to collaborate. It might be that the adversary is
secretly controlling the policy of one of the agents. We would need new
methods that assume some of the training agent’s partners might not be aligned
with the objectives of the group.

###  3.3 Methods of Evaluation

One aspect of AHT that might be overlooked is the evaluation of the agents at
test time. An AHT agent might have to act differently in the presence of an
adversary instead of a low-skill partner. In general, we would need AHT
methods that are more risk-aware and define what success looks like when there
is a high probability that some of the other agents will try to subvert the
group’s goal.

##  4 Future Research Directions

AHT was established with robust multi-agent behavior as the main goal. Some
past works in AHT have made steps in that direction. Type-based methods, like
in Rahman et al. (2021), attempt to infer the agent’s behavior from its
actions. Agents can also use communication to help each other quickly adapt to
new scenarios (Macke et al., 2021), which as been successful in defending
against attacks in multi-agent deep learning (Tu et al., 2021). Extending past
work in AHT can show us which methods have the most potential in being robust
against adversaries. Past works in adversarial training in RL (Pinto et al.,
2017; Pan et al., 2019) show that robustness can be achieved if the agent is
trained with adversarial examples.

Extending past methods alone, however, will likely not be enough. Introducing
new agents to other previously trained agents exacerbate the brittleness of
current RL methods (Gleave et al., 2019). The success of adversarial AHT, and
AHT in general, will depend on a foundation of risk-awareness for multi-agent
RL in the presence of adversaries. This foundation needs to account for the
“openness” of the environment, which models how new agents enter the group.
The presence of adversaries compels us to carefully understand how openness
impacts the training of an AHT agent. Openness that includes adversarial
agents will require researchers to assume additional formal descriptions of
the environment, which include (1) the distribution of new adversarial agents
to introduce to the group, (2) the diversity of the adversaries, and (3) the
reward function that accounts for the number and types of adversaries. As
multi-agent RL is starting to be applied to more real-world tasks (Wang et
al., 2021; Bae & Koumoutsakos, 2022), there is now more motivation to train
agents that are robust and risk-aware. This creates an incentive to train AHT
agents in high-stakes environments. An adversarial ML perspective can make
these endeavors more successful.

##  5 Conclusion

As RL continues to be pursued to accomplish complex, real-world tasks, we will
inevitably need to train RL agents that robustly cooperate with others. Since
AHT seems to be a promising path toward that direction, we propose accounting
for the possible presence of adversaries as a necessary component of that
field. We put forward that adversarial ML researchers will be critical to
ensuring AHT will advance safe collaboration between both humans and RL
agents.

## References

  * Athalye et al. (2018) Athalye, A., Carlini, N., and Wagner, D.  Obfuscated gradients give a false sense of security: Circumventing defenses to adversarial examples.  In _International conference on machine learning_ , pp. 274–283. PMLR, 2018.
  * Bae & Koumoutsakos (2022) Bae, H. J. and Koumoutsakos, P.  Scientific multi-agent reinforcement learning for wall-models of turbulent flows.  _Nature Communications_ , 13(1):1–9, 2022.
  * Biagioni et al. (2021) Biagioni, D., Zhang, X., Wald, D., Vaidhynathan, D., Chintala, R., King, J., and Zamzam, A. S.  Powergridworld: A framework for multi-agent reinforcement learning in power systems.  _arXiv preprint arXiv:2111.05969_ , 2021.
  * Dafoe et al. (2020) Dafoe, A., Hughes, E., Bachrach, Y., Collins, T., McKee, K. R., Leibo, J. Z., Larson, K., and Graepel, T.  Open problems in cooperative ai.  _arXiv preprint arXiv:2012.08630_ , 2020.
  * Dommès et al. (2021) Dommès, A., Merlhiot, G., Lobjois, R., Dang, N.-T., Vienne, F., Boulo, J., Oliver, A.-H., Cretual, A., and Cavallo, V.  Young and older adult pedestrians’ behavior when crossing a street in front of conventional and self-driving cars.  _Accident Analysis & Prevention_, 159:106256, 2021.
  * Figura et al. (2021) Figura, M., Kosaraju, K. C., and Gupta, V.  Adversarial attacks in consensus-based multi-agent reinforcement learning.  In _2021 American Control Conference (ACC)_ , pp. 3050–3055. IEEE, 2021.
  * Fujimoto et al. (2021) Fujimoto, T., Doster, T., Attarian, A., Brandenberger, J., and Hodas, N.  Reward-free attacks in multi-agent reinforcement learning.  _arXiv preprint arXiv:2112.00940_ , 2021.
  * Genter et al. (2017) Genter, K., Laue, T., and Stone, P.  Three years of the robocup standard platform league drop-in player competition.  _Autonomous Agents and Multi-Agent Systems_ , 31(4):790–820, 2017.
  * Gleave et al. (2019) Gleave, A., Dennis, M., Wild, C., Kant, N., Levine, S., and Russell, S.  Adversarial policies: Attacking deep reinforcement learning.  In _International Conference on Learning Representations_ , 2019.
  * Hou et al. (2021) Hou, J., Wang, F., Wang, L., and Chen, Z.  Reinforcement learning based multi-agent resilient control: From deep neural networks to an adaptive law.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 35, pp. 7737–7745, 2021.
  * Hu et al. (2021) Hu, H., Lerer, A., Cui, B., Pineda, L., Brown, N., and Foerster, J.  Off-belief learning.  In _International Conference on Machine Learning_ , pp. 4369–4379. PMLR, 2021.
  * Huang et al. (2017) Huang, S., Papernot, N., Goodfellow, I., Duan, Y., and Abbeel, P.  Adversarial attacks on neural network policies.  _arXiv preprint arXiv:1702.02284_ , 2017.
  * Januário et al. (2019) Januário, F., Cardoso, A., and Gil, P.  A distributed multi-agent framework for resilience enhancement in cyber-physical systems.  _IEEE Access_ , 7:31342–31357, 2019.
  * Li et al. (2019) Li, S., Wu, Y., Cui, X., Dong, H., Fang, F., and Russell, S.  Robust multi-agent reinforcement learning via minimax deep deterministic policy gradient.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 33, pp. 4213–4220, 2019.
  * Lin et al. (2020) Lin, J., Dzeparoska, K., Zhang, S. Q., Leon-Garcia, A., and Papernot, N.  On the robustness of cooperative multi-agent reinforcement learning.  In _2020 IEEE Security and Privacy Workshops (SPW)_ , pp. 62–68. IEEE, 2020.
  * Lowe et al. (2017) Lowe, R., Wu, Y. I., Tamar, A., Harb, J., Pieter Abbeel, O., and Mordatch, I.  Multi-agent actor-critic for mixed cooperative-competitive environments.  _Advances in neural information processing systems_ , 30, 2017.
  * Lupu et al. (2021) Lupu, A., Cui, B., Hu, H., and Foerster, J.  Trajectory diversity for zero-shot coordination.  In _International Conference on Machine Learning_ , pp. 7204–7213. PMLR, 2021.
  * Macke et al. (2021) Macke, W., Mirsky, R., and Stone, P.  Expected value of communication for planning in ad hoc teamwork.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 35, pp. 11290–11298, 2021.
  * Mirsky et al. (2022) Mirsky, R., Carlucho, I., Rahman, A., Fosong, E., Macke, W., Sridharan, M., Stone, P., and Albrecht, S. V.  A survey of ad hoc teamwork: Definitions, methods, and open problems.  _arXiv preprint arXiv:2202.10450_ , 2022.
  * Pan et al. (2019) Pan, X., Seita, D., Gao, Y., and Canny, J.  Risk averse robust adversarial reinforcement learning.  In _2019 International Conference on Robotics and Automation (ICRA)_ , pp. 8522–8528. IEEE, 2019.
  * Phan et al. (2021) Phan, T., Belzner, L., Gabor, T., Sedlmeier, A., Ritz, F., and Linnhoff-Popien, C.  Resilient multi-agent reinforcement learning with adversarial value decomposition.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 35, pp. 11308–11316, 2021.
  * Pigott et al. (2021) Pigott, A., Crozier, C., Baker, K., and Nagy, Z.  Gridlearn: Multiagent reinforcement learning for grid-aware building energy management.  _arXiv preprint arXiv:2110.06396_ , 2021.
  * Pinto et al. (2017) Pinto, L., Davidson, J., Sukthankar, R., and Gupta, A.  Robust adversarial reinforcement learning.  In _International Conference on Machine Learning_ , pp. 2817–2826. PMLR, 2017.
  * Rahman et al. (2021) Rahman, M. A., Hopner, N., Christianos, F., and Albrecht, S. V.  Towards open ad hoc teamwork using graph-based policy learning.  In _International Conference on Machine Learning_ , pp. 8776–8786. PMLR, 2021.
  * Rashid et al. (2018) Rashid, T., Samvelyan, M., Schroeder, C., Farquhar, G., Foerster, J., and Whiteson, S.  Qmix: Monotonic value function factorisation for deep multi-agent reinforcement learning.  In _International Conference on Machine Learning_ , pp. 4295–4304. PMLR, 2018.
  * Shafipour & Fallah (2021) Shafipour, E. and Fallah, S.  Task-based ad-hoc teamwork with adversary.  In _Annual Conference Towards Autonomous Robotic Systems_ , pp. 76–87. Springer, 2021.
  * Short et al. (2019) Short, A., La Pay, T., and Gandhi, A.  Defending against adversarial examples.  Technical report, Sandia National Lab.(SNL-NM), Albuquerque, NM (United States), 2019.
  * Stone et al. (2010) Stone, P., Kaminka, G. A., Kraus, S., and Rosenschein, J. S.  Ad hoc autonomous agent teams: Collaboration without pre-coordination.  In _Twenty-Fourth AAAI Conference on Artificial Intelligence_ , 2010\.
  * Strouse et al. (2021) Strouse, D., McKee, K., Botvinick, M., Hughes, E., and Everett, R.  Collaborating with humans without human data.  _Advances in Neural Information Processing Systems_ , 34, 2021.
  * Sun et al. (2020) Sun, J., Zhang, T., Xie, X., Ma, L., Zheng, Y., Chen, K., and Liu, Y.  Stealthy and efficient adversarial attacks against deep reinforcement learning.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 34, pp. 5883–5891, 2020.
  * Treutlein et al. (2021) Treutlein, J., Dennis, M., Oesterheld, C., and Foerster, J.  A new formalism, method and open issues for zero-shot coordination.  In _International Conference on Machine Learning_ , pp. 10413–10423. PMLR, 2021.
  * Tu et al. (2021) Tu, J., Wang, T., Wang, J., Manivasagam, S., Ren, M., and Urtasun, R.  Adversarial attacks on multi-agent communication.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pp. 7768–7777, 2021.
  * Wang et al. (2021) Wang, J., Xu, W., Gu, Y., Song, W., and Green, T. C.  Multi-agent reinforcement learning for active voltage control on power distribution networks.  _Advances in Neural Information Processing Systems_ , 34:3271–3284, 2021.