2\mathsf{S}_{X}(x,\nabla_{\textsf{{x}}}\hat{f}(x,y^{\circ}),\mathaccentV{bar}0\symboldoperators
16\lambda_{1})$.

###### Proof.

Inequality 66 follows from [54, Proposition 5.5] after recovering the constant
factors from the proof (see [54, Appendix A]), bounding the second term in the
right-hand side of the resulting inequality (cf. [54, Eq. (5.6)]) via Cauchy-
Schwarz, and extracting the square root. For $y^{\circ}=y^{\circ}(x)$ the
corresponding term vanishes, and we get
$\|\nabla\hat{\varphi}_{2\mathaccentV{bar}0\symboldoperators
16\lambda_{1}}(x)\|\leqslant
2\mathsf{S}_{X}(x,\nabla_{\textsf{{x}}}\hat{f}(x,y^{\circ}),\mathaccentV{bar}0\symboldoperators
16\lambda_{1})$ directly from [54, Eq. (5.6)]. We obtain 65 by the same
argument, as in this case $\hat{f}\equiv\hat{f}_{0}$ is formally maximized at
$\hat{y}$ for any $x\in X$. Note that the argument cannot be extended to
$k\geqslant 2$: in this case $\hat{f}_{k}(x,\cdot)$ is not concave, so [54,
Proposition 5.5] cannot be applied anymore. ∎

###  C.1 Proof of Theorem 5.1

We first observe that, with $\gamma_{\textsf{{x}}}={1}/{\lambda}$, the
quantity $\varepsilon_{t}$ computed in 7 of Algorithm 1 is nothing else but
$\mathsf{S}_{X}(x_{t},\nabla_{\textsf{{x}}}\hat{f}(x_{t},\hat{y}),\lambda)$.
Indeed:

$\displaystyle\mathsf{S}_{X}^{2}(x_{t},\nabla_{\textsf{{x}}}f(x_{t},\hat{y}),\lambda)$
$\displaystyle=2\lambda\max_{x\in
X}\big{\\{}-\langle\nabla_{\textsf{{x}}}f(x_{t},\hat{y}),x-x_{t}\rangle-\tfrac{\lambda}{2}\|x-x_{t}\|^{2}\big{\\}}$
(67) $\displaystyle=\tfrac{1}{\gamma_{\textsf{{x}}}^{2}}\max_{x\in
X}\big{\\{}-2\langle\gamma_{\textsf{{x}}}\nabla_{\textsf{{x}}}f(x_{t},\hat{y}),x-x_{t}\rangle-\|x-x_{t}\|^{2}\big{\\}}$
$\displaystyle=\|\nabla_{\textsf{{x}}}f(x_{t},\hat{y})\|^{2}-\tfrac{1}{\gamma_{\textsf{{x}}}^{2}}\min_{x\in
X}\|x_{t}-\gamma_{\textsf{{x}}}\nabla_{\textsf{{x}}}f(x_{t},\hat{y})-x\|^{2}=\varepsilon_{t}^{2}.$

Thus, Algorithm 1 maintains $x^{*}$ such that
$\mathsf{S}_{X}(x^{*},\nabla_{\textsf{{x}}}f(x^{*},\hat{y}),\lambda)=\min_{0\leqslant
t\leqslant
T-1}\mathsf{S}_{X}(x_{t},\nabla_{\textsf{{x}}}f(x_{t},\hat{y}),\lambda).$ On
the other hand, by the descent lemma (cf. Assumption 1), for each
$t\in\\{0,...,T-1\\}$ we have that

$\displaystyle f(x_{t+1},\hat{y})$ $\displaystyle\leqslant
f(x_{t},\hat{y})+\langle\nabla_{\textsf{{x}}}f(x_{t},\hat{y}),x_{t+1}-x_{t}\rangle+\tfrac{\lambda}{2}\|x_{t+1}-x_{t}\|^{2}$
$\displaystyle=f(x_{t},\hat{y})+\min_{x\in
X}\left\\{\langle\nabla_{\textsf{{x}}}f(x_{t},\hat{y}),x-x_{t}\rangle+\tfrac{\lambda}{2}\|x-x_{t}\|^{2}\right\\}$
$\displaystyle=f(x_{t},\hat{y})-\tfrac{1}{2\lambda}\mathsf{S}_{X}^{2}(x_{t},\nabla_{\textsf{{x}}}f(x_{t},\hat{y}),\lambda),$

where the first equality follows from the definition of $x_{t+1}$. Whence, via
telescoping and 29 we get

$\displaystyle\mathsf{S}_{X}^{2}(x^{*},\nabla_{\textsf{{x}}}f(x^{*},\hat{y}),\lambda)$
$\displaystyle=\min_{0\leqslant t\leqslant
T-1}\mathsf{S}_{X}^{2}(x_{t},\nabla_{\textsf{{x}}}f(x_{t},\hat{y}),\lambda)\leqslant\frac{1}{T}\sum_{t=0}^{T-1}\mathsf{S}_{X}^{2}(x_{t},\nabla_{\textsf{{x}}}f(x_{t},\hat{y}),\lambda)$
$\displaystyle\leqslant\frac{2\lambda\left[f(x_{0},\hat{y})-f(x_{T},\hat{y})\right]}{T}\leqslant\frac{2\lambda\left[\varphi(x_{0})-\psi(\hat{y})\right]}{T}<\frac{\varepsilon^{2}}{144}.$

By the first claim of Lemma C.1 (cf. 65) this gives
$\|\nabla\hat{\varphi}_{2\lambda}(x^{*})\|\leqslant
2\varepsilon^{*}\leqslant\varepsilon/6$. Finally, 30 implies the premise of
Theorem 3.1 with $k=0$; applying it we arrive at
$\|\nabla\varphi_{2\lambda}(x^{*})\|\leqslant\varepsilon$. ∎

###  C.2 Proof of Theorem 5.2

$\boldsymbol{{1}^{o}}$. We first consider the case
$\mu\leqslant(\mathaccentV{bar}0\symboldoperators
16\lambda_{1}\rho_{1})^{1/2}$, so that $\textsf{{Coupled}}=0$ and condition 32
reduces to

$200\mu\mathsf{D}\leqslant{\varepsilon}.$ (68)

As $\hat{f}:=\hat{f}_{1}$, 8 of Algorithm 2 reads
$y_{t}\in\operatorname*{Argmax}_{y\in Y}\hat{f}(x_{t},y)$; then the second
claim of Lemma C.1 (with $y_{t}=y^{\circ}$) guarantees that

$\|\nabla\hat{\varphi}_{2\mathaccentV{bar}0\symboldoperators
16\lambda_{1}}(x_{t})\|\leqslant
2\mathsf{S}_{X}(x_{t},\nabla_{\textsf{{x}}}\hat{f}(x_{t},y_{t}),\mathaccentV{bar}0\symboldoperators
16\lambda_{1}).$ (69)

Let us now upper-bound
$\mathsf{S}_{X}(x_{t},\nabla_{\textsf{{x}}}\hat{f}(x_{t},y_{t}),\mathaccentV{bar}0\symboldoperators
16\lambda_{1})$. To this end, for $\delta>0$ to be chosen later, consider

$\displaystyle
f_{\textsf{{reg}}}(x,y):=f(x,y)-\frac{\delta}{2}\|y-\hat{y}\|^{2},\quad\varphi_{\textsf{{reg}}}(x):=\max_{y\in
Y}f_{\textsf{{reg}}}(x,y).$ (70)

Since $f_{\textsf{{reg}}}(x,\cdot)$ is $\delta$-strongly concave, by Danskin’s
theorem (cf. [53, Lemma 24]) $\varphi_{\textsf{{reg}}}$ is differentiable;
$\nabla\varphi_{\textsf{{reg}}}(x)=\nabla_{\textsf{{x}}}f(x,y_{\textsf{{reg}}}(x))$
with $y_{\textsf{{reg}}}(x):=\operatorname*{argmax}_{y\in
Y}f_{\textsf{{reg}}}(x,y)$, and is $\lambda_{\textsf{{reg}}}$-Lipschitz with
$\lambda_{\textsf{{reg}}}:=\lambda+\frac{\mu^{2}}{\delta}.$ Now, let us choose
$\delta=\frac{\mu^{2}}{\mathaccentV{bar}0\symboldoperators 16\lambda_{1}},$ so
that $\lambda_{\textsf{{reg}}}=\lambda+\mathaccentV{bar}0\symboldoperators
16\lambda_{1}\in[\mathaccentV{bar}0\symboldoperators
16\lambda_{1},2\mathaccentV{bar}0\symboldoperators 16\lambda_{1}]$. By the
descent lemma

$\displaystyle\varphi_{\textsf{{reg}}}(x_{t+1})-\varphi_{\textsf{{reg}}}(x_{t})$
$\displaystyle\leqslant\langle\nabla\varphi_{\textsf{{reg}}}(x_{t}),x_{t+1}-x_{t}\rangle+\frac{\lambda_{\textsf{{reg}}}}{2}\|x_{t+1}-x_{t}\|^{2}=\langle\nabla_{\textsf{{x}}}f(x_{t},y_{\textsf{{reg}}}(x_{t})),x_{t+1}-x_{t}\rangle+\frac{\lambda_{\textsf{{reg}}}}{2}\|x_{t+1}-x_{t}\|^{2}$
$\displaystyle\leqslant\langle\nabla_{\textsf{{x}}}f(x_{t},y_{t}),x_{t+1}-x_{t}\rangle+\frac{3\lambda_{\textsf{{reg}}}}{4}\|x_{t+1}-x_{t}\|^{2}+\frac{1}{\lambda_{\textsf{{reg}}}}\|\nabla_{\textsf{{x}}}f(x_{t},y_{\textsf{{reg}}}(x_{t}))-\nabla_{\textsf{{x}}}f(x_{t},y_{t})\|^{2}$
$\displaystyle\leqslant\langle\nabla_{\textsf{{x}}}f(x_{t},y_{t}),x_{t+1}-x_{t}\rangle+\frac{3\lambda_{\textsf{{reg}}}}{4}\|x_{t+1}-x_{t}\|^{2}+\frac{\mu^{2}\mathsf{D}^{2}}{\lambda_{\textsf{{reg}}}}$
$\displaystyle\leqslant\langle\nabla_{\textsf{{x}}}f(x_{t},y_{t}),x_{t+1}-x_{t}\rangle+\frac{1}{2}\left({3\mathaccentV{bar}0\symboldoperators
16\lambda_{1}}+\frac{\mu^{2}}{\rho_{1}}\right)\|x_{t+1}-x_{t}\|^{2}+\frac{\mu^{2}\mathsf{D}^{2}}{\mathaccentV{bar}0\symboldoperators
16\lambda_{1}}$ $\displaystyle=\min_{x\in
X}\left\\{\langle\nabla_{\textsf{{x}}}f(x_{t},y_{t}),x-x_{t}\rangle+\frac{1}{2}\left({3\mathaccentV{bar}0\symboldoperators
16\lambda_{1}}+\frac{\mu^{2}}{\rho_{1}}\right)\|x-x_{t}\|^{2}\right\\}+\frac{\mu^{2}\mathsf{D}^{2}}{\mathaccentV{bar}0\symboldoperators
16\lambda_{1}}$ $\displaystyle=-\frac{1}{6\mathaccentV{bar}0\symboldoperators
16\lambda_{1}+2\mu^{2}/\rho_{1}}\mathsf{S}_{X}^{2}\left(x_{t},\nabla_{\textsf{{x}}}f(x_{t},y_{t}),{3\mathaccentV{bar}0\symboldoperators
16\lambda_{1}}+\frac{\mu^{2}}{\rho_{1}}\right)+\frac{\mu^{2}\mathsf{D}^{2}}{\mathaccentV{bar}0\symboldoperators
16\lambda_{1}}$
$\displaystyle\leqslant-\frac{1}{8\mathaccentV{bar}0\symboldoperators
16\lambda_{1}}\mathsf{S}_{X}^{2}\left(x_{t},\nabla_{\textsf{{x}}}f(x_{t},y_{t}),{3\mathaccentV{bar}0\symboldoperators
16\lambda_{1}}+\frac{\mu^{2}}{\rho_{1}}\right)+\frac{\mu^{2}\mathsf{D}^{2}}{\mathaccentV{bar}0\symboldoperators
16\lambda_{1}}.$ (71)

Here the second inequality is by Cauchy-Schwarz, the next one is via
Assumption 1, and the subsequent identities are by the definitions of
$x_{t+1}$ and $\gamma_{\textsf{{x}}}$ (cf. 11 of Algorithm 2). (Note that the
factor $\mu^{2}/\rho_{1}$ can be upper-bounded with
$\mathaccentV{bar}0\symboldoperators 16\lambda_{1}$ by the standing
assumption, but we avoid this in order for our estimates to have a similar
form as in the strongly coupled case (cf. $\boldsymbol{{2}^{o}}$) which is yet
to be considered.) We now proceed as in 67: by telescoping Section C.2 we get

$\displaystyle\mathsf{S}_{X}^{2}\left(x^{*},\nabla_{\textsf{{x}}}f(x^{*},y^{*}),3\mathaccentV{bar}0\symboldoperators
16\lambda_{1}+\frac{\mu^{2}}{\rho_{1}}\right)=\min_{0\leqslant t\leqslant
T-1}\mathsf{S}_{X}^{2}\left(x_{t},\nabla_{\textsf{{x}}}f(x_{t},y_{t}),3\mathaccentV{bar}0\symboldoperators
16\lambda_{1}+\frac{\mu^{2}}{\rho_{1}}\right)$
$\displaystyle\leqslant\frac{1}{T}\sum_{t=0}^{T-1}\mathsf{S}_{X}^{2}\left(x_{t},\nabla_{\textsf{{x}}}f(x_{t},y_{t}),3\mathaccentV{bar}0\symboldoperators
16\lambda_{1}+\frac{\mu^{2}}{\rho_{1}}\right)\leqslant\frac{8\mathaccentV{bar}0\symboldoperators
16\lambda_{1}[\varphi_{\textsf{{reg}}}(x_{0})-\varphi_{\textsf{{reg}}}(x_{T})]}{T}+8\mu^{2}\mathsf{D}^{2}$
$\displaystyle\leqslant\frac{8\mathaccentV{bar}0\symboldoperators
16\lambda_{1}[\varphi(x_{0})-\varphi(x_{T})]}{T}+\frac{4\mathaccentV{bar}0\symboldoperators
16\lambda_{1}\delta\mathsf{D}^{2}}{T}+8\mu^{2}\mathsf{D}^{2}\leqslant\frac{8\mathaccentV{bar}0\symboldoperators
16\lambda_{1}[\varphi(x_{0})-\varphi(x_{T})]}{T}+12\mu^{2}\mathsf{D}^{2}.$
(72)

As a result,

$\displaystyle\mathsf{S}_{X}^{2}(x^{*},\nabla_{\textsf{{x}}}\hat{f}(x^{*},y^{*}),\mathaccentV{bar}0\symboldoperators
16\lambda_{1})$
$\displaystyle\leqslant\mathsf{S}_{X}^{2}(x^{*},\nabla_{\textsf{{x}}}\hat{f}(x^{*},y^{*}),4\mathaccentV{bar}0\symboldoperators
16\lambda_{1})$ $\displaystyle=8\mathaccentV{bar}0\symboldoperators
16\lambda_{1}\max_{x\in
X}\left\\{-\langle\nabla_{\textsf{{x}}}\hat{f}(x^{*},y^{*}),x-x^{*}\rangle-2\mathaccentV{bar}0\symboldoperators
16\lambda_{1}\|x-x^{*}\|^{2}\right\\}$ $\displaystyle\leqslant
8\mathaccentV{bar}0\symboldoperators 16\lambda_{1}\max_{x\in
X}\left\\{-\langle\nabla_{\textsf{{x}}}f(x^{*},y^{*}),x-x^{*}\rangle-\frac{3\mathaccentV{bar}0\symboldoperators
16\lambda_{1}\|x-x^{*}\|^{2}}{2}\right\\}$
$\displaystyle\quad\quad\quad+4\|\nabla_{\textsf{{x}}}f(x^{*},y^{*})-\nabla_{\textsf{{x}}}\hat{f}(x^{*},y^{*})\|^{2}$
$\displaystyle=\frac{4}{3}\mathsf{S}_{X}^{2}(x^{*},\nabla_{\textsf{{x}}}f(x^{*},y^{*}),3\mathaccentV{bar}0\symboldoperators
16\lambda_{1})+4\|\nabla_{\textsf{{x}}}f(x^{*},y^{*})-\nabla_{\textsf{{x}}}\hat{f}(x^{*},y^{*})\|^{2}$
$\displaystyle\leqslant\frac{4}{3}\mathsf{S}_{X}^{2}(x^{*},\nabla_{\textsf{{x}}}f(x^{*},y^{*}),3\mathaccentV{bar}0\symboldoperators
16\lambda_{1})+16\mu^{2}\mathsf{D}^{2}$
$\displaystyle\leqslant\frac{32\mathaccentV{bar}0\symboldoperators
16\lambda_{1}[\varphi(x_{0})-\varphi(x_{T})]}{3T}+32\mu^{2}\mathsf{D}^{2}\leqslant
11\varepsilon^{2}\left(\frac{1}{2100}+\frac{3}{40000}\right)<\frac{\varepsilon^{2}}{144}.$

Here the first estimate relies on the proximal Polyak-Łojasiewicz (PL)
inequality that ensures that
$\mathsf{S}_{X}(x,\xi,\mathaccentV{bar}0\symboldoperators 16\lambda)$ is non-
decreasing in $\mathaccentV{bar}0\symboldoperators 16\lambda$ (cf. [93, Lemma
1]); the second one is by Cauchy-Schwarz, the third one is by Lemma 3.2, the
fourth one is by 72, and the last one is by (33)–(68). Finally, returning to
69 we get $\|\nabla\hat{\varphi}_{2\mathaccentV{bar}0\symboldoperators
16\lambda_{1}}(x^{*})\|\leqslant{\varepsilon}/{6}$, and by Theorem 3.1 this
results in $\|\nabla\varphi_{2\mathaccentV{bar}0\symboldoperators
16\lambda_{1}}(x^{*})\|\leqslant\varepsilon$.

$\boldsymbol{{2}^{o}}$. We now consider the case
$\mu\geqslant(\mathaccentV{bar}0\symboldoperators
16\lambda_{1}\rho_{1})^{1/2}$, where Algorithm 2 is run with
$\textsf{{Coupled}}=1$. By 32,

$200\sqrt{\mathaccentV{bar}0\symboldoperators
16\lambda_{1}\rho_{1}}\mathsf{D}\leqslant\varepsilon.$ (73)

Casting the update in 5 of Algorithm 2 as $y_{t}=\operatorname*{argmax}_{y\in
Y}\\{\langle\nabla_{\textsf{{y}}}f(x_{t},\hat{y}),y-\hat{y}\rangle-\frac{\rho_{1}}{2}\|y-\hat{y}\|^{2}\\}$
gives

$\left\langle\nabla_{\textsf{{y}}}f(x_{t},\hat{y})-\rho_{1}(y_{t}-\hat{y}),y-y_{t}\right\rangle\leqslant
0,\quad\forall y\in Y,$

from to the first-order optimality condition. As a result, for any $0\leqslant
t\leqslant T-1$ we have

$\displaystyle\mathsf{S}_{Y}^{2}(y_{t},-\nabla_{\textsf{{y}}}\hat{f}(x_{t},y_{t}),\rho_{1})=\mathsf{S}_{Y}^{2}(y_{t},-\nabla_{\textsf{{y}}}f(x_{t},\hat{y}),\rho_{1})$
$\displaystyle=2\rho_{1}\max_{y\in
Y}\left\\{\langle\nabla_{\textsf{{y}}}f(x_{t},\hat{y}),y-y_{t}\rangle-\frac{\rho_{1}}{2}\|y-y_{t}\|^{2}\right\\}$
$\displaystyle\leqslant 2\rho_{1}^{2}\max_{y\in
Y}\left\\{-\frac{1}{2}\|y-y_{t}\|^{2}+\langle
y_{t}-\hat{y},y-y_{t}\rangle\right\\}\leqslant\rho_{1}^{2}\mathsf{D}^{2}.$

where in the final step we maximized over the whole $E_{{\textsf{{y}}}}$.
Whence by Lemma C.1 (cf. 66) we get

$\displaystyle\|\nabla\hat{\varphi}_{2\mathaccentV{bar}0\symboldoperators
16\lambda_{1}}(x_{t})\|$
$\displaystyle\;\leqslant\;2\mathsf{S}_{X}(x_{t},\nabla_{\textsf{{x}}}\hat{f}(x_{t},y_{t}),\mathaccentV{bar}0\symboldoperators
16\lambda_{1})+4\sqrt{\lambda\rho_{1}}\mathsf{D}.$ (74)

Our next goal is to estimate
$\mathsf{S}_{X}(x_{t},\nabla_{\textsf{{x}}}\hat{f}(x_{t},y_{t}),\mathaccentV{bar}0\symboldoperators
16\lambda_{1})$ via a telescoping argument. To this end, let

$\displaystyle\widetilde{f}_{\textsf{{reg}}}(x,y):=\hat{f}_{1}(x,y)-\frac{\rho_{1}}{2}\|y-\hat{y}\|^{2},\quad\widetilde{\varphi}_{\textsf{{reg}}}(x):=\max_{y\in
Y}\widetilde{f}_{\textsf{{reg}}}(x,y).$ (75)

Recall that $\nabla_{\textsf{{x}}}\hat{f}_{1}(\cdot,y)$ is
$\mathaccentV{bar}0\symboldoperators 16\lambda_{1}$-Lipschitz for any $y\in Y$
(cf. Lemma 3.3); moreover, $\widetilde{f}_{\textsf{{reg}}}(x,\cdot)$ is
$\rho_{1}$-strongly concave. Therefore by Danskin’s theorem (cf. [53, Lemma
24]),
$\nabla\widetilde{\varphi}_{\textsf{{reg}}}(x)=\nabla_{\textsf{{x}}}\hat{f}(x,\widetilde{y}_{\textsf{{reg}}}(x))$
with $\widetilde{y}_{\textsf{{reg}}}(x):=\operatorname*{argmax}_{y\in
Y}\widetilde{f}_{\textsf{{reg}}}(x,y)$, and is
$\widetilde{\lambda}_{\textsf{{reg}}}$-Lipschitz with
$\widetilde{\lambda}_{\textsf{{reg}}}:=\mathaccentV{bar}0\symboldoperators
16\lambda_{1}+{\mu^{2}}/{\rho_{1}}.$ Moreover, we have

$y_{t}=\widetilde{y}_{\textsf{{reg}}}(x_{t}),$

as seen by looking at 5 of Algorithm 2 again. Whence, by the descent lemma
(cf. Section C.2) we get

$\displaystyle\widetilde{\varphi}_{\textsf{{reg}}}(x_{t+1})-\widetilde{\varphi}_{\textsf{{reg}}}(x_{t})$
$\displaystyle\leqslant\langle\nabla\widetilde{\varphi}_{\textsf{{reg}}}(x_{t}),x_{t+1}-x_{t}\rangle+\frac{\widetilde{\lambda}_{\textsf{{reg}}}}{2}\|x_{t+1}-x_{t}\|^{2}$
$\displaystyle=\langle\nabla_{\textsf{{x}}}\hat{f}(x_{t},y_{t}),x_{t+1}-x_{t}\rangle+\frac{\widetilde{\lambda}_{\textsf{{reg}}}}{2}\|x_{t+1}-x_{t}\|^{2}$
$\displaystyle\leqslant\langle\nabla_{\textsf{{x}}}\hat{f}(x_{t},y_{t}),x_{t+1}-x_{t}\rangle+\frac{1}{2}\left(3\mathaccentV{bar}0\symboldoperators
16\lambda_{1}+\frac{\mu^{2}}{\rho_{1}}\right)\|x_{t+1}-x_{t}\|^{2}$
$\displaystyle=\min_{x\in
X}\left\\{\langle\nabla_{\textsf{{x}}}\hat{f}(x_{t},y_{t}),x-x_{t}\rangle+\frac{1}{2}\left({3\mathaccentV{bar}0\symboldoperators
16\lambda_{1}}+\frac{\mu^{2}}{\rho_{1}}\right)\|x-x_{t}\|^{2}\right\\}$
$\displaystyle\leqslant-\frac{1}{6\mathaccentV{bar}0\symboldoperators
16\lambda_{1}+2\mu^{2}/\rho_{1}}\mathsf{S}_{X}^{2}\left(x_{t},\nabla_{\textsf{{x}}}\hat{f}(x_{t},y_{t}),3\mathaccentV{bar}0\symboldoperators
16\lambda_{1}+\frac{\mu^{2}}{\rho_{1}}\right).$

Here the penultimate line follows by recasting the update in 6 of Algorithm 2
as

$\widetilde{x}_{t+1}\left[=x_{t}-\gamma_{\textsf{{x}}}\nabla_{\textsf{{x}}}f(x_{t},\hat{y})-\gamma_{\textsf{{x}}}\nabla^{2}_{{\textsf{{x}}}{\textsf{{y}}}}f(x_{t},\hat{y})(y_{t}-\hat{y})\right]=x_{t}-\gamma_{\textsf{{x}}}\nabla_{\textsf{{x}}}\hat{f}(x_{t},y_{t}).$

Moreover, for the same reason we have (cf. 67)

$\displaystyle\varepsilon_{t}^{2}$
$\displaystyle=\|\nabla_{\textsf{{x}}}\hat{f}(x_{t},y_{t})\|^{2}-\tfrac{1}{\gamma_{\textsf{{x}}}^{2}}\min_{x\in
X}\|x_{t}-\gamma_{\textsf{{x}}}\nabla_{\textsf{{x}}}\hat{f}(x_{t},y_{t})-x\|^{2}=\mathsf{S}_{X}^{2}\left(x_{t},\nabla_{\textsf{{x}}}\hat{f}(x_{t},y_{t}),3\mathaccentV{bar}0\symboldoperators
16\lambda_{1}+\frac{\mu^{2}}{\rho_{1}}\right).$ (76)

Whence, proceeding as in 72 we get

$\displaystyle\mathsf{S}_{X}^{2}\left(x^{*},\nabla_{\textsf{{x}}}\hat{f}(x^{*},y^{*}),3\mathaccentV{bar}0\symboldoperators
16\lambda_{1}+\frac{\mu^{2}}{\rho_{1}}\right)$ $\displaystyle=\min_{0\leqslant
t\leqslant
T-1}\mathsf{S}_{X}^{2}\left(x_{t},\nabla_{\textsf{{x}}}\hat{f}(x_{t},y_{t}),3\mathaccentV{bar}0\symboldoperators
16\lambda_{1}+\frac{\mu^{2}}{\rho_{1}}\right)$
$\displaystyle\leqslant\frac{1}{T}\sum_{t=0}^{T-1}\mathsf{S}_{X}^{2}\left(x_{t},\nabla_{\textsf{{x}}}\hat{f}(x_{t},y_{t}),3\mathaccentV{bar}0\symboldoperators
16\lambda_{1}+\frac{\mu^{2}}{\rho_{1}}\right)$
$\displaystyle\leqslant\left(6\mathaccentV{bar}0\symboldoperators
16\lambda_{1}+\frac{2\mu^{2}}{\rho_{1}}\right)\frac{\widetilde{\varphi}_{\textsf{{reg}}}(x_{0})-\widetilde{\varphi}_{\textsf{{reg}}}(x_{T})}{T}$
$\displaystyle\stackrel{{\scriptstyle(i)}}{{\leqslant}}\frac{2}{T}\left(3\mathaccentV{bar}0\symboldoperators
16\lambda_{1}+\frac{\mu^{2}}{\rho_{1}}\right)\left(\hat{\varphi}(x_{0})-\hat{\varphi}(x_{T})+\frac{\rho_{1}\mathsf{D}^{2}}{2}\right)$
$\displaystyle\stackrel{{\scriptstyle(ii)}}{{\leqslant}}\frac{2}{T}\left(3\mathaccentV{bar}0\symboldoperators
16\lambda_{1}+\frac{\mu^{2}}{\rho_{1}}\right)\left(\varphi(x_{0})-\varphi(x_{T})+\frac{3\rho_{1}\mathsf{D}^{2}}{2}\right).$

Here for $(i)$ we used the two inequalities
$\widetilde{\varphi}_{\textsf{{reg}}}(x_{0})-\widetilde{\varphi}_{\textsf{{reg}}}(x_{T})\leqslant\hat{\varphi}(x_{0})-\max_{y\in
Y}\\{\hat{f}_{1}(x_{T},y)-\frac{\rho_{1}}{2}\|y-\hat{y}\|^{2}\\}$ and
$\hat{\varphi}(x_{T})\leqslant\max_{y\in
Y}\\{\hat{f}_{1}(x_{T},y)-\frac{\rho_{1}}{2}\|y-\hat{y}\|^{2}\\}+\frac{\rho_{1}\mathsf{D}^{2}}{2}$
(cf. 75); for $(ii)$ we applied Lemma 3.1 to the right-hand side of
$|\hat{\varphi}(x)-\varphi(x)|\leqslant\max_{y\in
Y}|f(x,y)-\hat{f}_{1}(x,y)|$. Returning to 74 we now obtain that

$\displaystyle\|\nabla\hat{\varphi}_{2\mathaccentV{bar}0\symboldoperators
16\lambda_{1}}(x^{*})\|$ $\displaystyle\leqslant
2\mathsf{S}_{X}\left(x_{t},\nabla_{\textsf{{x}}}\hat{f}(x_{t},y_{t}),3\mathaccentV{bar}0\symboldoperators
16\lambda_{1}+\frac{\mu^{2}}{\rho_{1}}\right)+4\sqrt{\lambda\rho_{1}}\mathsf{D}$
$\displaystyle\leqslant\sqrt{\frac{8}{T}\left(3\mathaccentV{bar}0\symboldoperators
16\lambda_{1}+\frac{\mu^{2}}{\rho_{1}}\right)\left(\varphi(x_{0})-\varphi(x_{T})+\frac{3\rho_{1}\mathsf{D}^{2}}{2}\right)}+4\sqrt{\lambda\rho_{1}}\mathsf{D}$
$\displaystyle\stackrel{{\scriptstyle\lx@cref{creftype~refnum}{eq:first-order-
strong-coupling-
condition}}}{{\leqslant}}\sqrt{\frac{8}{T}\left(3\mathaccentV{bar}0\symboldoperators
16\lambda_{1}+\frac{\mu^{2}}{\rho_{1}}\right)\left(\varphi(x_{0})-\varphi(x_{T})+\frac{3\varepsilon^{2}}{80000\mathaccentV{bar}0\symboldoperators
16\lambda_{1}}\right)}+\frac{\varepsilon}{50}$
$\displaystyle\leqslant\sqrt{\frac{8}{T}\left(3\mathaccentV{bar}0\symboldoperators
16\lambda_{1}+\frac{\mu^{2}}{\rho_{1}}\right)\left(\varphi(x_{0})-\varphi(x_{T})\right)}+\frac{3\varepsilon}{100}\sqrt{\frac{1}{T}\left(1+\frac{\mu^{2}}{3\mathaccentV{bar}0\symboldoperators
16\lambda_{1}\rho_{1}}\right)}+\frac{\varepsilon}{50}$
$\displaystyle\stackrel{{\scriptstyle\lx@cref{creftype~refnum}{eq:first-order-
complexity}}}{{\leqslant}}\varepsilon\left(\sqrt{\frac{8}{700}}+\frac{2\sqrt{3}}{100}+\frac{1}{50}\right)<\frac{\varepsilon}{6}.$

Finally, by applying Theorem 3.1 we conclude that
$\|\nabla\varphi_{2\mathaccentV{bar}0\symboldoperators
16\lambda_{1}}(x^{*})\|\leqslant\varepsilon$ as required. ∎

###  C.3 Proof of Proposition 5.1

We first observe that $\hat{f}(\cdot,y)\equiv\hat{f}_{2}(\cdot,y)$ is
$(\sigma_{0}+\sigma_{2}\mathsf{D}^{2})$-Lipschitz for any $y\in Y$ as can be
seen from

$\|\nabla_{\textsf{{x}}}\hat{f}(x,y)\|\leqslant\|\nabla_{\textsf{{x}}}f(x,y)\|+\|\nabla_{\textsf{{x}}}\hat{f}(x,y)-\nabla_{\textsf{{x}}}f(x,y)\|\leqslant\sigma_{0}+\sigma_{2}\mathsf{D}^{2},$
(77)

where the last step is by Lemma 3.2. Moreover, by Lemma 3.3
$\nabla_{\textsf{{x}}}\hat{f}(\cdot,y)$ is
$\mathaccentV{bar}0\symboldoperators 16\lambda_{2}$-Lipschitz, and thus
$\hat{\varphi}$ is $\mathaccentV{bar}0\symboldoperators 16\lambda_{2}$-weakly
convex. These two observations allow us to adapt the analysis of the projected
subgradient method from [61, Theorem 31] (initially carried out in [62]) to
Algorithm 3 in order to establish that $x^{*}$ is an approximate FOSP for
${\mathsf{P}}_{2}$ – and after that allude to Theorem 3.1.

For brevity, let $\hat{f}\equiv\hat{f}_{2}$ and $\hat{\varphi}=\max_{y\in
Y}\hat{f}(x,y)$. Observe that, for any $x\in X$ and iterate $(x_{t},y_{t})$,

$\displaystyle\hat{\varphi}(x)\geqslant\hat{f}(x,y_{t})$
$\displaystyle\geqslant\hat{f}(x_{t},y_{t})+\langle\nabla_{\textsf{{x}}}\hat{f}(x_{t},y_{t}),x-x_{t}\rangle-\frac{\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}{2}\|x-x_{t}\|^{2}$
$\displaystyle\geqslant\hat{\varphi}(x_{t})-\delta+\langle\nabla_{\textsf{{x}}}\hat{f}(x_{t},y_{t}),x-x_{t}\rangle-\frac{\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}{2}\|x-x_{t}\|^{2},$ (78)

where the last step is by definition of ApproxMax. Moreover, let
$x_{t}^{+}=\operatorname*{argmin}_{x\in
X}\\{\hat{\varphi}(x)+\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\|x-x_{t}\|^{2}\\}$ be the proximal mapping of $x_{t}$, so that
(cf. 9)

$\nabla\hat{\varphi}_{2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}(x_{t}^{\vphantom{+}})=2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}(x_{t}^{\vphantom{+}}-x_{t}^{+}).$ (79)

$\boldsymbol{{1}^{o}}$. Let Algorithm 3 be run with $\textsf{{Naive}}=0$,
i.e.,
$x_{t+1}=\Pi_{X}[x_{t}-\gamma_{\textsf{{x}}}\nabla_{\textsf{{x}}}\hat{f}(x_{t},y_{t})]$
(cf. 7). Then

$\displaystyle\hat{\varphi}_{2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}(x_{t+1})$
$\displaystyle\stackrel{{\scriptstyle\phantom{\lx@cref{creftype~refnum}{eq:quad-
grad-
bound}}}}{{\leqslant}}\hat{\varphi}(x_{t}^{+})+\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\|x_{t+1}^{\vphantom{+}}-x_{t}^{+}\|^{2}$ (80)
$\displaystyle\stackrel{{\scriptstyle\phantom{\lx@cref{creftype~refnum}{eq:quad-
grad-
bound}}}}{{\leqslant}}\hat{\varphi}(x_{t}^{+})+\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\|x_{t}-\gamma_{\textsf{{x}}}\nabla_{\textsf{{x}}}\hat{f}(x_{t},y_{t})-x_{t}^{+}\|^{2}$
$\displaystyle\stackrel{{\scriptstyle\lx@cref{creftype~refnum}{eq:quad-grad-
bound}}}{{\leqslant}}\hat{\varphi}(x_{t}^{+})+\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\|x_{t}^{\vphantom{+}}-x_{t}^{+}\|^{2}+2\gamma_{\textsf{{x}}}\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\langle\nabla_{\textsf{{x}}}\hat{f}(x_{t},y_{t}),x_{t}^{+}-x_{t}^{\vphantom{+}}\rangle+\gamma_{\textsf{{x}}}^{2}\mathaccentV{bar}0\symboldoperators
16\lambda_{2}(\sigma_{0}+\sigma_{2}\mathsf{D}^{2})^{2}$
$\displaystyle\stackrel{{\scriptstyle\phantom{\lx@cref{creftype~refnum}{eq:quad-
grad-bound}}}}{{=}}\hat{\varphi}_{2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}(x_{t})+2\gamma_{\textsf{{x}}}\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\langle\nabla_{\textsf{{x}}}\hat{f}(x_{t},y_{t}),x_{t}^{+}-x_{t}^{\vphantom{+}}\rangle+\gamma_{\textsf{{x}}}^{2}\mathaccentV{bar}0\symboldoperators
16\lambda_{2}(\sigma_{0}+\sigma_{2}\mathsf{D}^{2})^{2}$
$\displaystyle\stackrel{{\scriptstyle\lx@cref{creftype~refnum}{eq:quad-fun-
decrease}}}{{\leqslant}}\hat{\varphi}_{2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}(x_{t})+2\gamma_{\textsf{{x}}}\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\left(\hat{\varphi}(x_{t}^{+})-\hat{\varphi}(x_{t}^{\vphantom{+}})+\delta+\frac{\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}{2}\|x_{t}^{+}-x_{t}^{\vphantom{+}}\|^{2}\right)+\gamma_{\textsf{{x}}}^{2}\mathaccentV{bar}0\symboldoperators
16\lambda_{2}(\sigma_{0}+\sigma_{2}\mathsf{D}^{2})^{2},$

where the second line relies on the projection lemma. Repeating this for
$t\in\\{0,...,T-1\\}$ we get

$\hat{\varphi}_{2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}(x_{T})\leqslant\hat{\varphi}_{2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}(x_{0})+2\gamma_{\textsf{{x}}}\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\sum_{t=0}^{T-1}\left(\hat{\varphi}(x_{t}^{+})-\hat{\varphi}(x_{t}^{\vphantom{+}})+\delta+\frac{\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}{2}\|x_{t}^{+}-x_{t}^{\vphantom{+}}\|^{2}\right)+\gamma_{\textsf{{x}}}^{2}T\mathaccentV{bar}0\symboldoperators
16\lambda_{2}(\sigma_{0}+\sigma_{2}\mathsf{D}^{2})^{2},$

whence, by rearranging and dividing over
$2\gamma_{\textsf{{x}}}\mathaccentV{bar}0\symboldoperators 16\lambda_{2}T$,

$\frac{1}{T}\sum_{t=0}^{T-1}\hat{\varphi}(x_{t}^{\vphantom{+}})-\hat{\varphi}(x_{t}^{+})-\frac{\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}{2}\|x_{t}^{+}-x_{t}^{\vphantom{+}}\|^{2}\leqslant\frac{\hat{\varphi}_{2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}(x_{0})-\hat{\varphi}_{2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}(x_{T})}{2\gamma_{\textsf{{x}}}\mathaccentV{bar}0\symboldoperators
16\lambda_{2}T}+\frac{\gamma_{\textsf{{x}}}(\sigma_{0}+\sigma_{2}\mathsf{D}^{2})^{2}}{2}+\delta.$
(81)

On the other hand, by $\mathaccentV{bar}0\symboldoperators
16\lambda_{2}$-strong convexity of
$\hat{\varphi}(\cdot)+\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\|\cdot-x_{t}\|^{2}$ and the definition of $x_{t}^{+}$ we get

$\displaystyle\hat{\varphi}(x_{t}^{\vphantom{+}})-\hat{\varphi}(x_{t}^{+})-\frac{\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}{2}\|x_{t}^{+}-x_{t}^{\vphantom{+}}\|^{2}$
$\displaystyle=\hat{\varphi}(x_{t}^{\vphantom{+}})+\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\|x_{t}-x_{t}\|^{2}-(\hat{\varphi}(x_{t}^{+})+\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\|x_{t}^{+}-x_{t}^{\vphantom{+}}\|^{2})+\frac{\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}{2}\|x_{t}^{+}-x_{t}^{\vphantom{+}}\|^{2}$ (82)
$\displaystyle=\hat{\varphi}(x_{t}^{\vphantom{+}})+\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\|x_{t}-x_{t}\|^{2}-\min_{x\in
X}\left\\{\hat{\varphi}(x_{t}^{+})+\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\|x_{t}^{+}-x_{t}^{\vphantom{+}}\|^{2}\right\\}+\frac{\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}{2}\|x_{t}^{+}-x_{t}^{\vphantom{+}}\|^{2}$
$\displaystyle\geqslant\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\|x_{t}^{+}-x_{t}^{\vphantom{+}}\|^{2}\stackrel{{\scriptstyle\lx@cref{creftype~refnum}{eq:quad-
moreau-
step}}}{{=}}\frac{\|\nabla\hat{\varphi}_{2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}(x_{t}^{\vphantom{+}})\|^{2}}{4\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}.$

Plugging this into 81 we arrive at

$\displaystyle\mathds{E}[\|\nabla\hat{\varphi}_{2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}(x_{s}^{\vphantom{+}})\|^{2}]=\frac{1}{T}\sum_{t=0}^{T-1}\|\nabla\hat{\varphi}_{2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}(x_{t}^{\vphantom{+}})\|^{2}$
$\displaystyle\leqslant\frac{2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\left[\hat{\varphi}_{2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}(x_{0})-\hat{\varphi}_{2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}(x_{T})\right]}{\gamma_{\textsf{{x}}}T}+2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\gamma_{\textsf{{x}}}(\sigma_{0}+\sigma_{2}\mathsf{D}^{2})^{2}+4\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\delta$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leqslant}}\frac{2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\left[\hat{\varphi}(x_{0})-\hat{\varphi}(x_{T})\right]}{\gamma_{\textsf{{x}}}T}+2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\gamma_{\textsf{{x}}}(\sigma_{0}+\sigma_{2}\mathsf{D}^{2})^{2}+4\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\delta$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leqslant}}\frac{2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}(\Delta+\rho_{2}\mathsf{D}^{3})}{\gamma_{\textsf{{x}}}T}+2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\gamma_{\textsf{{x}}}(\sigma_{0}+\sigma_{2}\mathsf{D}^{2})^{2}+4\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\delta;$ (83)

here for $(a)$ we used that
$\hat{\varphi}_{2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}(x_{0})=\min_{x\in
X}\\{\hat{\varphi}(x)+\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\|x-x_{0}\|^{2}\\}\leqslant\hat{\varphi}(x_{0})$; for $(b)$ we
used Lemma 3.1. Whence by Markov’s inequality, for any ${p}\in(0,1)$ with
probability at least $1-{p}$ one has

$\displaystyle\|\nabla\hat{\varphi}_{2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}(x_{s}^{\vphantom{+}})\|^{2}$
$\displaystyle\leqslant\frac{4}{{p}}\left(\frac{\mathaccentV{bar}0\symboldoperators
16\lambda_{2}(\Delta+\rho_{2}\mathsf{D}^{3})}{\gamma_{\textsf{{x}}}T}+\gamma_{\textsf{{x}}}\mathaccentV{bar}0\symboldoperators
16\lambda_{2}(\sigma_{0}+\sigma_{2}\mathsf{D}^{2})^{2}+2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\delta\right)$
$\displaystyle\leqslant\frac{8}{{p}}\left((\sigma_{0}+\sigma_{2}\mathsf{D}^{2})\sqrt{\frac{\mathaccentV{bar}0\symboldoperators
16\lambda_{2}(\Delta+\rho_{2}\mathsf{D}^{3})}{T}}+\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\delta\right)<\frac{\varepsilon^{2}}{150},$

where for the second line we substituted $\gamma_{\textsf{{x}}}$ from 35 (note
that this choice of $\gamma_{\textsf{{x}}}$ balances the two terms) and then
performed an explicit calculation by plugging in $T$ and $\delta$ from 35.
Finally, under this event 34 allows to apply Theorem 3.1 (replacing
$\varepsilon$ with $\varepsilon/6$) and conclude that
$\|\nabla\varphi_{2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}(x_{s})\|\leqslant\varepsilon$.

$\boldsymbol{{2}^{o}}$. Conforming with the second claim we intend to prove,
from now on we shall assume that

$24\sigma_{2}\mathsf{D}^{2}\leqslant\varepsilon\sqrt{{p}},$ (84)

and consider the iterates of Algorithm 3 with $\textsf{{Naive}}=1$—i.e.,
$x_{t+1}=\Pi_{X}[x_{t}-\gamma_{\textsf{{x}}}\nabla_{\textsf{{x}}}f(x_{t},y_{t})]$
(cf. 5). First of all, let us correct 80 for the discrepancy between
$\nabla_{\textsf{{x}}}f(x_{t},y_{t})$ and
$\nabla_{\textsf{{x}}}\hat{f}(x_{t},y_{t})$: to this end,

$\displaystyle\hat{\varphi}_{2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}(x_{t+1})$
$\displaystyle\leqslant\hat{\varphi}(x_{t}^{+})+\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\|x_{t+1}^{\vphantom{+}}-x_{t}^{+}\|^{2}$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leqslant}}\hat{\varphi}(x_{t}^{+})+\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\|x_{t}-\gamma_{\textsf{{x}}}\nabla_{\textsf{{x}}}f(x_{t},y_{t})-x_{t}^{+}\|^{2}$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leqslant}}\hat{\varphi}(x_{t}^{+})+\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\left(\|x_{t}^{+}-x_{t}^{\vphantom{+}}\|^{2}+2\gamma_{\textsf{{x}}}\langle\nabla_{\textsf{{x}}}f(x_{t},y_{t}),x_{t}^{+}-x_{t}^{\vphantom{+}}\rangle+\gamma_{\textsf{{x}}}^{2}\sigma_{0}^{2}\right)$
$\displaystyle=\hat{\varphi}_{2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}(x_{t})+2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\gamma_{\textsf{{x}}}\langle\nabla_{\textsf{{x}}}f(x_{t},y_{t}),x_{t}^{+}-x_{t}^{\vphantom{+}}\rangle+\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\gamma_{\textsf{{x}}}^{2}\sigma_{0}^{2}$
$\displaystyle\stackrel{{\scriptstyle(c)}}{{\leqslant}}\hat{\varphi}_{2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}(x_{t})+2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\gamma_{\textsf{{x}}}\left(\langle\nabla_{\textsf{{x}}}\hat{f}(x_{t},y_{t}),x_{t}^{+}-x_{t}^{\vphantom{+}}\rangle+\sigma_{2}\mathsf{D}^{2}\|x_{t}^{+}-x_{t}^{\vphantom{+}}\|\right)+\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\gamma_{\textsf{{x}}}^{2}\sigma_{0}^{2}$
$\displaystyle\stackrel{{\scriptstyle(d)}}{{\leqslant}}\hat{\varphi}_{2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}(x_{t})+2\gamma_{\textsf{{x}}}\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\left(\langle\nabla_{\textsf{{x}}}\hat{f}(x_{t},y_{t}),x_{t}^{+}-x_{t}^{\vphantom{+}}\rangle+\frac{\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}{2}\|x_{t}^{+}-x_{t}^{\vphantom{+}}\|^{2}+\frac{\sigma_{2}^{2}\mathsf{D}^{4}}{2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}\right)+\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\gamma_{\textsf{{x}}}^{2}\sigma_{0}^{2}$
$\displaystyle\stackrel{{\scriptstyle\lx@cref{creftype~refnum}{eq:quad-fun-
decrease}}}{{\leqslant}}\hat{\varphi}_{2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}(x_{t})+2\gamma_{\textsf{{x}}}\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\left(\hat{\varphi}(x_{t}^{+})-\hat{\varphi}(x_{t}^{\vphantom{+}})+\delta+\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\|x_{t}^{+}-x_{t}^{\vphantom{+}}\|^{2}+\frac{\sigma_{2}^{2}\mathsf{D}^{4}}{2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}\right)+\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\gamma_{\textsf{{x}}}^{2}\sigma_{0}^{2}.$

Here $(a)$ is by the projection lemma; $(b)$ by the Lipschitz assumption;
$(c)$ by Cauchy-Schwarz and Lemma 3.2; $(d)$ by $uv\leqslant(u^{2}+v^{2})/2$
for $u,v\in\mathds{R}$. Whence we arrive at a counterpart of 81:

$\frac{1}{T}\sum_{t=0}^{T-1}\hat{\varphi}(x_{t}^{\vphantom{+}})-\hat{\varphi}(x_{t}^{+})-\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\|x_{t}^{+}-x_{t}^{\vphantom{+}}\|^{2}\leqslant\frac{\hat{\varphi}_{2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}(x_{0})-\hat{\varphi}_{2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}(x_{T})}{2\gamma_{\textsf{{x}}}\mathaccentV{bar}0\symboldoperators
16\lambda_{2}T}+\frac{\gamma_{\textsf{{x}}}\sigma_{0}^{2}}{2}+\delta+\frac{\sigma_{2}^{2}\mathsf{D}^{4}}{2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}.$ (85)

On the other hand, in the same spirit as for 82 we get

$\displaystyle\hat{\varphi}(x_{t}^{\vphantom{+}})-\hat{\varphi}(x_{t}^{+})-\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\|x_{t}^{+}-x_{t}^{\vphantom{+}}\|^{2}$
$\displaystyle=\hat{\varphi}(x_{t}^{\vphantom{+}})+\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\|x_{t}-x_{t}\|^{2}-(\hat{\varphi}(x_{t}^{+})+\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\|x_{t}^{+}-x_{t}^{\vphantom{+}}\|^{2})$ (86)
$\displaystyle=\hat{\varphi}(x_{t}^{\vphantom{+}})+\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\|x_{t}-x_{t}\|^{2}-\min_{x\in
X}\left\\{\hat{\varphi}(x_{t}^{+})+\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\|x_{t}^{+}-x_{t}^{\vphantom{+}}\|^{2}\right\\}$
$\displaystyle\geqslant\frac{\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}{2}\|x_{t}^{+}-x_{t}^{\vphantom{+}}\|^{2}\stackrel{{\scriptstyle\lx@cref{creftype~refnum}{eq:quad-
moreau-
step}}}{{=}}\frac{\|\nabla\hat{\varphi}_{2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}(x_{t}^{\vphantom{+}})\|^{2}}{8\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}.$

Plugging 86 into 85 and proceeding in the same way as in Section C.3, we
conclude that w.p. $\geqslant 1-{p}$,

$\displaystyle\|\nabla\hat{\varphi}_{2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}}(x_{s}^{\vphantom{+}})\|^{2}$
$\displaystyle\leqslant\frac{8}{{p}}\left(\frac{\mathaccentV{bar}0\symboldoperators
16\lambda_{2}(\Delta+\rho_{2}\mathsf{D}^{3})}{\gamma_{\textsf{{x}}}T}+\gamma_{\textsf{{x}}}\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\sigma_{0}^{2}+2\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\delta+\sigma_{2}^{2}\mathsf{D}^{4}\right)$
$\displaystyle\leqslant\frac{16}{{p}}\left(\sigma_{0}\sqrt{\frac{\mathaccentV{bar}0\symboldoperators
16\lambda_{2}(\Delta+\rho_{2}\mathsf{D}^{3})}{T}}+\mathaccentV{bar}0\symboldoperators
16\lambda_{2}\delta\right)+\frac{8\sigma_{2}^{2}\mathsf{D}^{4}}{{p}}\stackrel{{\scriptstyle\lx@cref{creftype~refnum}{eq:quad-
extra-
asm}}}{{\leqslant}}\frac{\varepsilon^{2}}{75}+\frac{\varepsilon^{2}}{72}<\frac{\varepsilon^{2}}{36}.$

In order to verify the claim, it remains to apply Theorem 3.1. ∎

##  Appendix D Background on weak convexity and the Moreau envelope notion

###  D.1 Characterization of weakly convex functions

Let $E_{{\textsf{{x}}}}$ be a finite-dimensional Euclidean space with inner
product $\left\langle\cdot,\cdot\right\rangle$ and norm
$\|x\|=\sqrt{\left\langle x,x\right\rangle}$. Moreover, we identify with
$E_{{\textsf{{x}}}}$ with its dual space $E_{{\textsf{{x}}}}^{*}$ (in
particular $\|\cdot\|^{*}=\|\cdot\|$). Also we let
$\mathaccentV{bar}0\symboldoperators 16\mathds{R}=(\infty,+\infty]$ and
tacitly assume that all arising convex functions are lsc and proper (not
identically $+\infty$). Recall that any such function
$\phi:E_{{\textsf{{x}}}}\to\mathaccentV{bar}0\symboldoperators 16\mathds{R}$
admits a variational representation as the pointwise supremum over a non-empty
family of affine functions, and this is in fact a criterion (see, e.g., [99,
Thm. 8.13]). The set of all affine minorants of $\phi$ at given $x\in
E_{{\textsf{{x}}}}$ generates the corresponding set of linear functionals,
called the subdifferential of $\phi$ at $x$ and denoted
$\partial\phi(x)$—i.e.,

$\partial\phi(x):=\\{\xi\in
E_{{\textsf{{x}}}}:\;\phi(x^{\prime})\geqslant\phi(x)+\langle\xi,x^{\prime}-x\rangle\;\;\forall
x^{\prime}\in E_{{\textsf{{x}}}}\\}.$ (87)

This set is non-empty at any interior point of effective domain ${\rm
dom}(\phi):=\\{x\in E_{{\textsf{{x}}}}:\phi(x)<+\infty\\}$, and its elements
are called subgradients of $\phi$ at $x$. (For simplicity, we assume ${\rm
int}({\rm dom}(\phi))$ is non-empty.)

######  Definition 3.

$\phi:E_{{\textsf{{x}}}}\to\mathaccentV{bar}0\symboldoperators 16\mathds{R}$
is called $\lambda$-weakly convex (for $\lambda\geqslant 0$) if the function
$\phi(\cdot)+\frac{\lambda}{2}\|\cdot\|^{2}$ is convex.

From the above characterization of (proper, lsc) convex functions in terms of
subdifferentials it immediately follows that $\phi$ is $\lambda$-weakly convex
if and only if the function
$\phi_{\lambda,x}(\cdot):=\phi(\cdot)+\tfrac{\lambda}{2}\|\cdot-x\|^{2}$ is
convex for any $x\in E_{{\textsf{{x}}}}$. Indeed, by Definition 3
$\phi_{\lambda,0}$ is convex; meanwhile, for any $x\in E_{{\textsf{{x}}}}$ we
have

$\phi_{\lambda,x}(u)=\phi_{\lambda,0}(u)+\tfrac{\lambda}{2}\left(\|u-x\|^{2}-\|u\|^{2}\right).$

By applying 87 to $\phi_{\lambda,0}$ we see that $\phi_{\lambda,x}$ has
subdifferential
$\partial\phi_{\lambda,x}(u)=\partial\phi_{\lambda,0}(u)-\lambda x$ (the
difference being in the Minkowski sense), and thus is convex with ${\rm
dom}(\phi_{\lambda,x})={\rm dom}(\phi)$. It is then natural to define the weak
subdifferential of $\phi$ at $x$ by

$\partial\phi(x)=\partial\phi_{\lambda,x}(x).$

Equivalently, $\partial\phi(u)=\partial\phi_{\lambda,x}(u)+\lambda(x-u)$ for
any $u,x\in E_{{\textsf{{x}}}}$—or, in other words,

$\partial\phi(x)=\\{\xi\in
E_{{\textsf{{x}}}}:\;\phi(x^{\prime})\geqslant\phi(x)+\langle\xi,x^{\prime}-x\rangle-\tfrac{\lambda}{2}\|x^{\prime}-x\|^{2},\;\;\forall
x,x^{\prime}\in E_{{\textsf{{x}}}}\\}.$ (88)

Thus, $\partial\phi(x)$ belongs to the Fréchet (local) subdifferential at $x$,
defined as the set of $\xi\in E_{{\textsf{{x}}}}$ satisfying

$\phi(x^{\prime})\geqslant\phi(x)+\langle\xi,x^{\prime}-x\rangle+o(\|x^{\prime}-x\|)\quad\text{as}\;\;x^{\prime}\to
x.$ (89)

cf. [99, Def. 8.3]. Generally, the property in 89 is much weaker than that in
88; however, for weakly convex functions the weak and Fréchet subdifferentials
coincide (see [99, Theorem 12.17]).

##### Hessian-based criterion.

Recall Alexandrov’s theorem ([100], see also [101, Thm. 2.1]): a convex
$\phi:E_{{\textsf{{x}}}}\to\mathaccentV{bar}0\symboldoperators 16\mathds{R}$
is twice differentiable almost everywhere (and hence
$\nabla^{2}\phi\succcurlyeq 0$) in the interior of its domain. Under $C^{1}$
regularity it is a criterion: if
$\phi:E_{{\textsf{{x}}}}\to\mathaccentV{bar}0\symboldoperators 16\mathds{R}$
is continuously differentiable on ${\rm int}({\rm dom}(\phi))$, and
$\nabla^{2}\phi\succcurlyeq 0$ almost everywhere on ${\rm int}({\rm
dom}(\phi))$, then $\phi$ is convex. From Definition 3 we see that this
criterion extends to $\lambda$-weakly convex functions if we replace
$\nabla^{2}\phi\succcurlyeq 0$ with $\nabla^{2}\phi\succcurlyeq-\lambda I$.

###  D.2 Upper envelope as a weakly convex function

Let $X\subseteq E_{{\textsf{{x}}}}$ be convex with non-empty interior. If
$\phi\to X$ is $\lambda$-smooth, i.e., such that $\nabla\phi$ exists and

$\|\nabla\phi(x^{\prime})-\nabla\phi(x)\|\leqslant\lambda\|x^{\prime}-x\|\quad\forall
x,x^{\prime}\in X,$

then $\phi$ is $\lambda$-weakly convex. (This follows from the Hessian-vased
criterion above, as $\nabla\phi(x)$ exists almost everywhere on $X$ and its
eigenvalues are in $[-\lambda,\lambda]$.) More generally, the upper envelope
of a family of such functions is $\lambda$-weakly convex. That is, given
$f:X\times Y$ (in the setup of $\mathsf{P}$) satisfying Assumption 1 with
$\lambda<\infty$, the primal function $\varphi(x):=\max_{y\in Y}f(x,y)$ is
$\lambda$-weakly convex. Indeed, note that
$f_{\lambda,x^{\prime}}(\cdot,y)=f(\cdot,y)+\frac{\lambda}{2}\|\cdot-x^{\prime}\|^{2}$
is convex for any $y\in Y$ (cf. Definition 3), therefore

$\varphi(x)+\tfrac{\lambda}{2}\|x-x^{\prime}\|^{2}=\max_{y\in
Y}f_{\lambda,x^{\prime}}(x,y)$

is convex, i.e., $\varphi$ is $\lambda$-weakly convex. Moreover,
$\partial\varphi(x)$ can be expressed explicitly through
$\nabla_{\textsf{{x}}}f(x,y)$.

######  Lemma D.1.

Let $f:X\times Y\to\mathds{R}$ be continuous and satisfy Assumption 1 with
$\lambda<\infty$, $X$ be convex, and $Y$ be convex and compact. Then
$\partial\varphi(x)$ is the closed convex hull of the set of active gradients:

$\partial\varphi(x)=\overline{\textup{conv}}\\{\nabla_{\textsf{{x}}}f(x,y^{*}),\forall
y^{*}\in\textstyle\operatorname*{Argmax}_{y\in Y}f(x,y)\\}.$

###### Proof.

We apply Danskin’s theorem ([102, Sec. B.5]) to the function
$\varphi_{\lambda,x^{\prime}}(x)=\max_{y\in Y}f_{\lambda,x^{\prime}}(x,y),$

$\displaystyle\partial\varphi_{\lambda,x^{\prime}}(x)$
$\displaystyle=\overline{\textup{conv}}\\{\nabla_{\textsf{{x}}}f_{\lambda,x^{\prime}}(x,y^{*}),y^{*}\in\textstyle\operatorname*{Argmax}_{y\in
Y}f_{\lambda,x^{\prime}}(x,y)\\},$
$\displaystyle=\overline{\textup{conv}}\\{\nabla_{\textsf{{x}}}f(x,y^{*})+\lambda(x-x^{\prime}),y^{*}\in\textstyle\operatorname*{Argmax}_{y\in
Y}f(x,y)\\}.$

This holds because each $f_{\lambda,x^{\prime}}(\cdot,y)$ is convex, and the
set of maximizers for $f_{\lambda,x^{\prime}}(x,y)$ (over $y\in Y$) is the
same as for $f(x,y)$. It remains to note that
$\partial\varphi(x)=\partial\varphi_{\lambda,x^{\prime}}(x)-\lambda(x-x^{\prime})$
by Definition 3. ∎

###  D.3 Characterization of $(\varepsilon,2\lambda)$-FOSP in terms of the
primal function

The following property of $(\varepsilon,2\lambda)$-FOSP for $\mathsf{P}$ has
been leveraged in the proof Proposition 3.2. It extends [91, Lem. 2.2] to the
case $X\neq E_{{\textsf{{x}}}}$ with guarantee 90 in terms of
$\partial\varphi$ rather than $\partial(\varphi+I_{X})$.

######  Proposition D.1 ([54, Proposition 5.1]).

Let $\phi:X\to\mathds{R}$ be $\lambda$-weakly convex, then it holds that

$\nabla\phi_{2\lambda}(x)=2\lambda(x-x^{+}(x)),\;\;\text{where}\;\;x^{+}(x)=x^{+}_{\frac{1}{2\lambda}\phi}(x)\;\;\text{for
any}\;\;x\in X.$

where $\phi_{2\lambda}$ is the $2\lambda$-Moreau envelope of $\phi$ (cf.
Definition 1). Moreover, if
$\|\nabla\phi_{2\lambda}(x)\|\leqslant\varepsilon$, then

$\min_{\xi\in\partial\phi(x^{+})}\mathsf{S}_{X}(x^{+},\xi,2\lambda)\leqslant\varepsilon,$
(90)

where we define the functional
$\mathsf{S}_{X}^{2}(x,\xi,\mathaccentV{bar}0\symboldoperators
16\lambda):=2\mathaccentV{bar}0\symboldoperators 16\lambda\max_{u\in
X}\big{\\{}-\langle\xi,u-x\rangle-\tfrac{\mathaccentV{bar}0\symboldoperators
16\lambda}{2}\|u-x\|^{2}\big{\\}}$ for $x,\xi\in E_{{\textsf{{x}}}}$,
$\mathaccentV{bar}0\symboldoperators 16\lambda>0$.

##  Appendix E Solving 40: quadratic optimization on a joint Krylov subspace

Following [92, Appendix A], in this section we describe a procedure for
maximizing a general quadratic function $\Psi_{H,g}(y)$, cf. 36, over the
intersection of $B_{d}(2\mathsf{R})$ and the joint Krylov subspace

$\mathcal{K}_{2m}=\mathcal{K}_{2m}(H,\\{g,\xi\\}):=\textup{span}\left(\\{H^{j}g,H^{j}\xi\\}_{j\in\\{0,...,m-1\\}}\right)$

in $O(m)$ computations of the matrix-vector product $y\mapsto Hy$ and
elementwise vector operations on $E_{\textsf{{y}}}$ (including the inner
product). This procedure is given in pseudocode in Algorithm 4. It consists of
two steps which we are now about to outline and discuss.

##### Finding the right basis for $\mathcal{K}_{2m}$.

In the first stage of the algorithm (up to 21) we construct an orthonormal
basis of $\mathcal{K}_{2m}$ in which the associated linear operator
(corresponding to $H$ in the initial basis), when restricted to
$\mathcal{K}_{2m}$, is represented by a block tridiagonal matrix
$\widetilde{H}$ with blocks of size $2$ (which pertains to the two vectors
$\\{g,\xi\\}$ generating $\mathcal{K}_{2m}$). In other words, we construct
$\widetilde{H}\in\mathds{R}^{2m\times 2m}$ and $Q\in\mathds{R}^{d\times 2m}$
such that $\widetilde{H}=Q^{\top}HQ,$ where $Q$ has orthonormal columns (i.e.,
$Q^{\top}Q=I_{2m}$) and

$\widetilde{H}=\begin{bmatrix}\alpha_{1}^{\vphantom{\top}}&\beta_{2}^{\top}&&\\\
\beta_{2}&\ddots&\ddots&\\\ &\ddots&\ddots&\beta_{m}^{\top}\\\
&&\beta_{m}&\alpha_{m}\end{bmatrix},$

where each $\alpha_{i}\in\mathds{R}^{2\times 2}$ is symmetric and each
$\beta_{i}\in\mathds{R}^{2\times 2}$ is upper-triangular (so $\widetilde{H}$
is pentadiagonal). In the large-scale context, where $H$ is accessed through
the $y\mapsto Hy$ oracle, the one can construct such a decomposition in a
computationally feasible manner via the block Lanczos process [103, 104]:

  1. 1.

Form $q_{1}\in\mathds{R}^{d\times 2}$ as the Gram-Schmidt orthonormalization
of $\\{g,\xi\\}$.

  2. 2.

Compute $\alpha_{1}=q_{1}^{\top}Hq_{1}$ and put $\beta_{1}=0_{2\times 2}$. (We
let $0_{a\times b}$ be the $a\times b$ zero matrix.)

  3. 3.

For $t\in\\{1,m-1\\}$ iterate:

$\displaystyle q_{t+1}^{\prime}$
$\displaystyle=Hq_{t}^{\vphantom{\top}}-q_{t}^{\vphantom{\top}}\alpha_{t}^{\vphantom{\top}}-q_{t-1}^{\vphantom{\top}}\beta_{t}^{\top},$
$\displaystyle(q_{t+1},\beta_{t+1})$
$\displaystyle=\text{QR}(q_{t+1}^{\prime}),$
$\displaystyle\alpha_{t+1}^{\vphantom{\top}}$
$\displaystyle=q_{t+1}^{\top}Hq_{t+1}^{\vphantom{\top}}.$

Here $\text{QR}(A)$ is the QR decomposition of a squared matrix $A$—i.e., the
ordered pair $(U,R)$ of square matrices such that $A=UR$, $U$ is orthogonal,
and $R$ is upper-triangular. To see that the above process is sound, observe
that it incrementally builds up $q_{t},\alpha_{t},\beta_{t}$ that satisfy the
matrix relations

$H[q_{1}\;\cdots\;q_{t}]=[q_{1}\;\cdots\;q_{t}\,|\,q_{t+1}]\begin{bmatrix}\alpha_{1}^{\vphantom{\top}}&\beta_{2}^{\top}&&\\\
\beta_{2}&\ddots&\ddots&\\\ &\ddots&\ddots&\beta_{t}^{\top}\\\
&&\beta_{t}&\alpha_{t}\\\
\hline\cr&&&\beta_{t+1}\end{bmatrix}\;\;\text{for}\;\;t\in\\{1,...,m-1\\},\quad\text{and}\quad
Hq_{m}=[q_{1}\;\cdots\;q_{m}]\begin{bmatrix}0\\\ \vdots\\\ 0\\\
\beta_{m}^{\top}\\\ \alpha_{m}\end{bmatrix}$

while ensuring that $\beta_{t}$’s are upper-triangular and $q_{t}$’s form an
orthonormal system. These relations amount to the identity $HQ=Q\widetilde{H}$
with $Q=[q_{1}\;\cdots\;q_{m}]$—and thus to $Q^{\top}HQ=\widetilde{H}$ as
desired. The first part of Algorithm 4 (until 21) implements the block Lanczos
process with an explicitly rendered QR factorization step.

##### Passing to the equivalent problem and solving it.

Once $\widetilde{H}$ and $Q$ have been constructed, we compute
$\widetilde{g}=Q^{\top}g$ and the top eigenvalue $\omega_{0}$ of
$\widetilde{H}$ (cf. 24). The first of these two steps amounts to $O(m)$ inner
products, and the second one can be done in $O(m)$ arithmetic operations since
$\widetilde{H}$ is pentadiagonal. Using these data, we recast the initial
maximization problem (cf. 40) by changing the variable $y$ to
$z=Q^{\top}y\in\mathds{R}^{2m}$—i.e., we pass to

$\max_{z\in
B_{2m}(2\mathsf{R})}\tfrac{1}{2}z^{\top}\widetilde{H}z+\widetilde{g}^{\top}z.$
(91)

The new problem is equivalent to 40 in the following sense: if $z$ is
$\delta$-suboptimal in 91, then $y=Qz$ is $\delta$-suboptimal in 40. Thus we
can focus on 91 and then recover a candidate solution to 40 by multiplying
with $Q$ (which has essentially the same computational cost as forming
$\widetilde{g}=Q^{\top}g$ from $g$).

The final part of Algorithm 4 (after 24) consists of solving 91. Here we
proceed as follows.

  1. 1.

We first test if $\widetilde{H}$ is negative-semidefinite (by inspecting
$\mathop{\rm sign}(\omega_{0})$) and, if yes, whether $B_{2m}(2\mathsf{R})$
contains an unconstrained maximizer $z_{0}$—an optimal solution to the linear
system $\widetilde{H}z+\widetilde{g}=0$. If the answers to both questions are
positive, we in fact found a maximizer on $B_{2m}(2\mathsf{R})$, so we output
$Qz_{0}$ and terminate. Note that here if $\widetilde{H}$ is negative-
definite, then $z_{0}$ exists, is unique, and is given by
$-\widetilde{H}^{-1}\widetilde{g}$; otherwise (i.e., when $\omega_{0}=0$),
$\widetilde{H}$ is rank-deficient, so $\widetilde{H}z+\widetilde{g}=0$ is
solvable whenever $\widetilde{g}$ is in the range of $\widetilde{H}$, and then
the least-norm solution is given by
$z_{0}=-\widetilde{H}^{\dagger}\widetilde{g}$, where $\widetilde{H}^{\dagger}$
is the generalized inverse of $\widetilde{H}$ (i.e., $\widetilde{H}^{\dagger}$
has the same eigenspaces as $\widetilde{H}$ and reciprocal nonzero
eigenvalues). Note that $\widetilde{H}^{\dagger}$ can be computed in $O(m)$
arithmetic operations (a.o.).

  2. 2.

If the above test failed (and thus we have not terminated), then
$\Psi_{\widetilde{H},\widetilde{g}}$ has a unique maximizer at the boundary of
$B_{2m}(2\mathsf{R})$, and this maximizer also solves the strictly concave
problem

$z_{\omega}=\operatorname*{argmax}_{z\in
B_{2m}(2\mathsf{R})}\tfrac{1}{2}z^{\top}(\widetilde{H}-\omega
I)z+\widetilde{g}^{\top}z$

for the unique value of $\omega>\max\\{\underline{\omega},0\\}$ such that
$\|z_{\omega}\|=\mathsf{R}$ (see, e.g., [105, Corollary 7.2.2]). Thus, we can
proceed by a root-finding method of choice: evaluation of the mapping
$\omega\mapsto\|z_{\omega}\|$ amounts to solving a well-defined linear system
with a pentadiagonal matrix and thus can be done in $O(m)$ a.o.; thus, we find
a high-accuracy maximizer—exact one in floating point arithmetic—in $O(m)$
a.o. In particular, a good option is practice is to seek for the root of
$\psi(\omega)-{1}/{\mathsf{R}}$ with $\psi(\omega)=1/\|z_{\omega}\|$ by
Newton’s method as discussed in [105, Chapter 7].

Algorithm 4 Approximate maximization of a quadratic
$\Psi_{H,g}(y)=\frac{1}{2}y^{\top}Hy+g^{\top}y$ on a ball
$Y=B_{d}(2\mathsf{R})$

1:$g\in\mathds{R}^{d}$, oracle $y\mapsto Hy$, $\mathsf{R}>0$, $m\in\mathds{N}$

2:$\triangleright$ Restrict search to the joint Krylov subspace
$\mathcal{K}_{2m}(H,\\{g,\xi\\})$, cf. 38, by using Proposition 5.2.

3:

4:$\triangleright$ Blockwise Lanczos iterations for
$\mathcal{K}_{2m}=\mathcal{K}_{2m}(H,\\{g,\xi\\})$ with blocks of size $2$
([92, Appendix A]); results in a block tridiagonal matrix
$\widetilde{H}=Q^{\top}HQ\in\mathds{R}^{2m\times 2m}$ with
$Q\in\mathds{R}^{d\times 2m}$: $Q^{\top}Q=I$.

5:$\triangleright$ $q_{t}\in\mathds{R}^{d\times 2}$,
$\alpha_{t},\beta_{t}\in\mathds{R}^{2\times 2}$

6:

7:$q_{0}=0_{d\times 2}$ $\triangleright$ $d\times 2$ zero matrix

8: Draw $\xi\sim\textup{{Uniform}}(\mathds{S}^{d-1})$; $u=\frac{1}{\|g\|}g$;
$w=\xi-\langle\xi,u\rangle u$; $v=\frac{1}{\|w\|}w$

9: $q_{1}=[u,v]$; $\alpha_{1}=q_{1}^{\top}Hq_{1}^{\vphantom{\top}}$
$\triangleright$ $q_{1}$ is an orthonormalization of $[g,\xi]$

10:

11:for $t\in\\{1,...,m-1\\}$ do

12:
$q_{t+1}^{\prime}=Hq_{t}-q_{t}\alpha_{t}-q_{t-1}^{\vphantom{\top}}\beta_{t}^{\top}$
$\triangleright$ $\beta_{1}$ never used since $q_{0}=0$

13: $\triangleright$ Compute
$(q_{t+1}^{\vphantom{{}^{\prime}}},\beta_{t+1}^{\vphantom{{}^{\prime}}})=\textup{QR}(q_{t+1}^{\prime})$
via Gram-Schmidt:

14: $[u^{\prime},v^{\prime}]=q_{t+1}^{\prime}$ $\triangleright$ extract
columns from $q_{t+1}^{\prime}$

15: $u=\frac{1}{\|u^{\prime}\|}u^{\prime}$; $w=v^{\prime}-\left\langle
v^{\prime},u\right\rangle u$; $v=\frac{1}{\|w\|}w$

16: $q_{t+1}=[u,v]$; $\beta_{t+1}=\begin{bmatrix}\left\langle
u^{\prime},u\right\rangle&\left\langle v^{\prime},u\right\rangle\\\
0&\left\langle v^{\prime},v\right\rangle\end{bmatrix}$ $\triangleright$
$q_{t+1}$ is an orthonormalization of
$q_{t+1}^{\prime}=q_{t+1}^{\vphantom{{}^{\prime}}}\beta_{t+1}^{\vphantom{{}^{\prime}}}$

17: $\alpha_{t+1}^{\vphantom{\top}}=q_{t+1}^{\top}Hq_{t+1}^{\vphantom{\top}}$

18:end for

19:

20:$Q=\begin{bmatrix}q_{1}\;\cdots\;q_{m}\end{bmatrix}$ $\triangleright$
columns of $Q$ form an orthonormal basis for $\mathcal{K}_{2m}$

21:$\widetilde{H}=\begin{bmatrix}\alpha_{1}^{\vphantom{\top}}&\beta_{2}^{\top}&&\\\
\beta_{2}&\ddots&\ddots&\\\ &\ddots&\ddots&\beta_{m}^{\top}\\\
&&\beta_{m}&\alpha_{m}\end{bmatrix}$ $\triangleright$
$\widetilde{H}=Q^{\top}HQ\in\mathds{R}^{2m\times 2m}$ is block tridiagonal
with $2\times 2$ blocks

22:

23:$\widetilde{g}=Q^{\top}g$

24:$\omega_{0}=\lambda_{\max}(\widetilde{H})$

25:

26:$\triangleright$ Eliminate the interior case: $\Psi_{H,g}$ restricted to
$\mathcal{K}_{2m}$ is concave and maximized inside the ball:

27:if $\omega_{0}<0$ then $\triangleright$ $\widetilde{H}$ is full-rank, so
the unconstrained maximizer $z_{0}$ exists and is unique

28: $z_{0}=-\widetilde{H}^{-1}\widetilde{g}$

29:end if

30:if $\omega_{0}==0$ and
$\widetilde{H}(\widetilde{H}^{2})^{\dagger}\widetilde{H}\widetilde{g}==\widetilde{g}$
then $\triangleright$ $\widetilde{H}$ is rank-deficient, but $\widetilde{g}$
is in its range

31: $z_{0}=-\widetilde{H}^{\dagger}\widetilde{g}$

32:end if

33:if $z_{0}$ is defined and $\|z_{0}\|\leqslant\mathsf{R}$ then

34: return $Qz_{0}$

35:end if

36:

37:$\triangleright$ If we have not terminated yet, the maximizer is on the
boundary ([105, Corollary 7.2.2])

38:Find $\omega>\max\\{\omega_{0},\,0\\}$ such that the norm of
$z_{\omega}=-(\widetilde{H}-\omega I)^{-1}\widetilde{g}$ is $\mathsf{R}$

39:return $Qz_{\omega}$

##### Time and memory expenses.

By carefully inspecting Algorithm 4 one may verify that the Lanczos process
takes $O(1)$ matrix-vector products and elementwise vector operations in
$E_{{\textsf{{y}}}}$ per iteration, so $O(m)$ such operations in total.
Computation of $\widetilde{g}$ takes $O(m)$ inner products, and those of
$\omega_{0}$, $\widetilde{H}^{\dagger}$, and $(\widetilde{H}-\omega I)^{-1}$
each take $O(m)$ a.o. Finally, as we previously discussed, a root finder
performs $O(1)$ evaluations of $\omega$ in the floating-point model of
computation. Hence, the total computational cost of Algorithm 4 in a.o. is

$O(m(d_{\textsf{{y}}}+\textsf{time$($HVP$)$})),$ (92)

where time$($HVP$)$ is the cost of computing $y\mapsto Hy$, and
$d_{\textsf{{y}}}=\dim(E_{{\textsf{{y}}}})$. As for the memory cost, it is
dominated by

$O(md_{\textsf{{y}}}+\textsf{mem$($HVP$)$}),$ (93)

where the first term allows to explicitly store $Q$ (which we need in order to
report the final solution $Qz_{0}$ or $Qz_{\omega}$), and the second term is
the memory needed for computing $y\mapsto Hy$ (i.e.,
$\textsf{mem$($HVP$)$}\lesssim d_{\textsf{{y}}}^{2}$ or less if $H$ has
special structure). Note also that we can replace the term $md_{\textsf{{y}}}$
in 93 with $m+d_{\textsf{{y}}}=O(d_{\textsf{{y}}})$ if Algorithm 4 in the
situation where there is no need to report the solution in
$E_{{\textsf{{y}}}}$ (e.g., if only the optimal value is of interest). Indeed,
apart from reporting the final solution, we only need $Q$ when computing
$\widetilde{g}=Q^{\top}g=[q_{1}^{\top}g;\;...;\;q_{m}^{\top}g]$—but this can
be done incrementally during the Lanczos process, without ever having to
memorize more than $O(1)$ columns of $Q$ at once.

## References

  * [1] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio.  Generative adversarial nets.  Advances in Neural Information Processing Systems, 27, 2014\.
  * [2] M. Sanjabi, M. Razaviyayn, J. Ba, and J. D. Lee.  On the convergence and robustness of training gans with regularized optimal transport.  Advances in Neural Information Processing Systems, 2018:7091–7101, 2018.
  * [3] S. Nowozin, B. Cseke, and R. Tomioka.  f-gan: Training generative neural samplers using variational divergence minimization.  In Proceedings of the 30th International Conference on Neural Information Processing Systems, pages 271–279, 2016.
  * [4] I. Gulrajani, F. Ahmed, M. Arjovsky, V. Dumoulin, and A. Courville.  Improved training of Wasserstein GANs.  In Proceedings of the 31st International Conference on Neural Information Processing Systems, pages 5769–5779, 2017.
  * [5] A. Madry, A. Makelov, L. Schmidt, D. Tsipras, and A. Vladu.  Towards deep learning models resistant to adversarial attacks.  In International Conference on Learning Representations, 2018.
  * [6] F. Pourbabaee.  Robust experimentation in the continuous time bandit problem.  Economic Theory, pages 1–31, 2020.
  * [7] H. Zhang, Y. Yu, J. Jiao, E. Xing, L. El Ghaoui, and M. Jordan.  Theoretically principled trade-off between robustness and accuracy.  In International Conference on Machine Learning, pages 7472–7482. PMLR, 2019.
  * [8] S. Sagawa, P. W. Koh, T. B. Hashimoto, and P. Liang.  Distributionally robust neural networks for group shifts: On the importance of regularization for worst-case generalization.  arXiv preprint arXiv:1911.08731, 2019.
  * [9] W. Hu, G. Niu, I. Sato, and M. Sugiyama.  Does distributionally robust supervised learning give robust classifiers?  In International Conference on Machine Learning, pages 2029–2037. PMLR, 2018.
  * [10] A. Sinha, H. Namkoong, and J. Duchi.  Certifying some distributional robustness with principled adversarial training.  In International Conference on Learning Representations, 2018.
  * [11] A. Rezaei, R. Fathony, O. Memarrast, and B. Ziebart.  Fairness for robust log loss classification.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 34, pages 5511–5518, 2020.
  * [12] S. Baharlouei, M. Nouiehed, A. Beirami, and M. Razaviyayn.  Rényi fair inference.  In 8th International Conference on Learning Representations (ICLR), 2020.
  * [13] B. H. Zhang, B. Lemoine, and M. Mitchell.  Mitigating unwanted biases with adversarial learning.  In Proceedings of the 2018 AAAI/ACM Conference on AI, Ethics, and Society, pages 335–340, 2018.
  * [14] D. Xu, S. Yuan, L. Zhang, and X. Wu.  Fairgan: Fairness-aware generative adversarial networks.  In 2018 IEEE International Conference on Big Data, pages 570–575. IEEE, 2018.
  * [15] A. Lowy, R. Pavan, S. Baharlouei, M. Razaviyayn, and A. Beirami.  Fermi: Fair empirical risk minimization via exponential Rényi mutual information.  arXiv preprint arXiv:2102.12586, 2021.
  * [16] T. Adel, I. Valera, Z. Ghahramani, and A. Weller.  One-network adversarial fairness.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 33, pages 2412–2420, 2019.
  * [17] B. Dai, A. Shaw, L. Li, L. Xiao, N. He, Z. Liu, J. Chen, and L. Song.  Sbeed: Convergent reinforcement learning with nonlinear function approximation.  In International Conference on Machine Learning, pages 1125–1134. PMLR, 2018.
  * [18] M. Rabbat and R. Nowak.  Distributed optimization in sensor networks.  In Proceedings of the 3rd international symposium on Information processing in sensor networks, pages 20–27, 2004.
  * [19] N. Li and J. R. Marden.  Designing games for distributed optimization.  IEEE Journal of Selected Topics in Signal Processing, 7(2):230–242, 2013.
  * [20] D. Hajinezhad, M. Hong, T. Zhao, and Z. Wang.  Nestt: A nonconvex primal-dual splitting method for distributed and stochastic optimization.  Advances in neural information processing systems, 29:3215–3223, 2016.
  * [21] M. Hong, D. Hajinezhad, and M.-M. Zhao.  Prox-pda: The proximal primal-dual algorithm for fast distributed nonconvex optimization and learning over networks.  In International Conference on Machine Learning, pages 1529–1538. PMLR, 2017.
  * [22] M. Hong, M. Razaviyayn, and J. Lee.  Gradient primal-dual algorithm converges to second-order stationary solution for nonconvex distributed optimization over networks.  In International Conference on Machine Learning, pages 2009–2018. PMLR, 2018.
  * [23] G. Scutari, D. P. Palomar, and S. Barbarossa.  Optimal linear precoding strategies for wideband noncooperative systems based on game theory—part i: Nash equilibria.  IEEE Transactions on Signal Processing, 56(3):1230–1249, 2008.
  * [24] G. Scutari and J.-S. Pang.  Joint sensing and power allocation in nonconvex cognitive radio games: Nash equilibria and distributed algorithms.  IEEE Transactions on Information Theory, 59(7):4626–4661, 2013\.
  * [25] J. F. Nash.  Equilibrium points in $n$-person games.  Proceedings of the National Academy of Sciences, 36(1):48–49, 1950.
  * [26] J. von Neumann.  Zur theorie der gesellschaftsspiele.  Mathematische annalen, 100(1):295–320, 1928.
  * [27] J. von Neumann and O. Morgenstern.  Theory of games and economic behavior.  Princeton university press, 1953.
  * [28] A. Nemirovski.  Prox-method with rate of convergence $O(1/t)$ for variational inequalities with Lipschitz continuous monotone operators and smooth convex-concave saddle point problems.  SIAM Journal on Optimization, 15(1):229–251, 2004.
  * [29] G. M. Korpelevich.  Extragradient method for finding saddle points and other problems.  Matekon, 13(4):35–49, 1977.
  * [30] A. Juditsky, A. Nemirovski, and C. Tauvel.  Solving variational inequalities with stochastic mirror-prox algorithm.  Stochastic Systems, 1(1):17–58, 2011.
  * [31] N. He, A. Juditsky, and A. Nemirovski.  Mirror prox algorithm for multi-term composite minimization and semi-separable problems.  Computational Optimization and Applications, 61(2):275–319, 2015\.
  * [32] N. He and Z. Harchaoui.  Semi-proximal mirror-prox for nonsmooth composite minimization.  In Proceedings of the 28th International Conference on Neural Information Processing Systems, volume 2, pages 3411–3419, 2015.
  * [33] A. Juditsky and A. Nemirovski.  Solving variational inequalities with monotone operators on domains given by linear minimization oracles.  Mathematical Programming, 156(1-2):221–256, 2016.
  * [34] B. Cox, A. Juditsky, and A. Nemirovski.  Decomposition techniques for bilinear saddle point problems and variational inequalities with affine monotone operators.  Journal of Optimization Theory and Applications, 172(2):402–435, 2017.
  * [35] L. Xiao.  Dual averaging methods for regularized stochastic learning and online optimization.  Journal of Machine Learning Research, 11(88):2543–2596, 2010.
  * [36] A. Juditsky and A. Nemirovski.  First order methods for nonsmooth convex large-scale optimization, ii: utilizing problem’s structure.  Optimization for Machine Learning, pages 149–183, 2011.
  * [37] Y. Nesterov and A. Nemirovski.  On first-order algorithms for $l_{1}$/nuclear norm minimization.  Acta Numerica, 22:509–575, 2013.
  * [38] Y. Wang and J. Li.  Improved algorithms for convex-concave minimax optimization.  arXiv preprint arXiv:2006.06359, 2020.
  * [39] J. Zhang, M. Hong, and S. Zhang.  On lower iteration complexity bounds for the convex concave saddle point problems.  Mathematical Programming, pages 1–35, 2021.
  * [40] C. Daskalakis and I. Panageas.  Last-iterate convergence: Zero-sum games and constrained min-max optimization.  arXiv preprint arXiv:1807.04252, 2018.
  * [41] J. Abernethy, K. A. Lai, and A. Wibisono.  Last-iterate convergence rates for min-max optimization.  arXiv preprint arXiv:1906.02027, 2019.
  * [42] N. Golowich, S. Pattathil, C. Daskalakis, and A. Ozdaglar.  Last iterate is slower than averaged iterate in smooth convex-concave saddle point problems.  In Conference on Learning Theory, pages 1758–1784. PMLR, 2020.
  * [43] C.-Y. Wei, C.-W. Lee, M. Zhang, and H. Luo.  Linear last-iterate convergence in constrained saddle-point optimization.  arXiv preprint arXiv:2006.09517, 2020.
  * [44] A. Juditsky and A. Nemirovski.  On well-structured convex-concave saddle point problems and variational inequalities with monotone operators.  arXiv preprint arXiv:2102.01002, 2021.
  * [45] M. Grant and S. Boyd.  CVX: Matlab software for disciplined convex programming, version 2.1.
  * [46] D. Kinderlehrer and G. Stampacchia.  An introduction to variational inequalities and their applications.  SIAM, 2000.
  * [47] L. Rice, E. Wong, and Z. Kolter.  Overfitting in adversarially robust deep learning.  In International Conference on Machine Learning, pages 8093–8104. PMLR, 2020.
  * [48] M. Razaviyayn, T. Huang, S. Lu, M. Nouiehed, M. Sanjabi, and M. Hong.  Nonconvex min-max optimization: Applications, challenges, and recent theoretical advances.  IEEE Signal Processing Magazine, 37(5):55–66, 2020.
  * [49] J. Zhang, P. Xiao, R. Sun, and Z.-Q. Luo.  A single-loop smoothed gradient descent-ascent algorithm for nonconvex-concave min-max problems.  arXiv preprint arXiv:2010.15768, 2020.
  * [50] S. Sagawa, P. W. Koh, T. B. Hashimoto, and P. Liang.  Distributionally robust neural networks for group shifts: On the importance of regularization for worst-case generalization.  International Conference on Learning Representations (ICLR), 2020\.
  * [51] S. Lu, I. Tsaknakis, M. Hong, and Y. Chen.  Hybrid block successive approximation for one-sided non-convex min-max problems: algorithms and applications.  IEEE Transactions on Signal Processing, 68:3676–3691, 2020.
  * [52] Q. Qian, S. Zhu, J. Tang, R. Jin, B. Sun, and H. Li.  Robust optimization over multiple domains.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 33, pages 4739–4746, 2019.
  * [53] M. Nouiehed, M. Sanjabi, J. D. Lee, and M. Razaviyayn.  Solving a class of non-convex min-max games using iterative first order methods.  Advances in Neural Information Processing Systems, 32, 2019.
  * [54] D. M. Ostrovskii, A. Lowy, and M. Razaviyayn.  Efficient search of first-order Nash equilibria in nonconvex-concave smooth min-max problems.  arXiv:2002.07919, 2020.
  * [55] K. K. Thekumparampil, P. Jain, P. Netrapalli, and S. Oh.  Efficient algorithms for smooth minimax optimization.  arXiv:1907.01543, 2019.
  * [56] T. Lin, C. Jin, and M. Jordan.  Near-optimal algorithms for minimax optimization.  arXiv preprint arXiv:2002.02417, 2020.
  * [57] R. Zhao.  A primal dual smoothing framework for max-structured nonconvex optimization.  arXiv preprint arXiv:2003.04375, 2020.
  * [58] Y. Han, G. Xie, and Z. Zhang.  Lower complexity bounds of finite-sum optimization problems: The results and construction.  arXiv preprint arXiv:2103.08280, 2021.
  * [59] S. Zhang, J. Yang, C. Guzmán, N. Kiyavash, and N. He.  The complexity of nonconvex-strongly-concave minimax optimization.  arXiv preprint arXiv:2103.15888, 2021.
  * [60] H. Li, Y. Tian, J. Zhang, and A. Jadbabaie.  Complexity lower bounds for nonconvex-strongly-concave min-max optimization.  arXiv preprint arXiv:2104.08708, 2021.
  * [61] C. Jin, P. Netrapalli, and M. Jordan.  What is local optimality in nonconvex-nonconcave minimax optimization?  In International Conference on Machine Learning, pages 4880–4889. PMLR, 2020.
  * [62] D. Davis and D. Drusvyatskiy.  Stochastic subgradient method converges at the rate $O(k^{-1/4})$ on weakly convex functions.  arXiv preprint arXiv:1802.02988, 2018.
  * [63] H. von Stackelberg.  Marktform und gleichgewicht.  J. Springer, 1934.
  * [64] C. Daskalakis, S. Skoulakis, and M. Zampetakis.  The complexity of constrained min-max optimization.  arXiv preprint arXiv:2009.09623, 2020.
  * [65] C. H. Papadimitriou.  On the complexity of the parity argument and other inefficient proofs of existence.  Journal of Computer and system Sciences, 48(3):498–532, 1994.
  * [66] C. Daskalakis, P. W. Goldberg, and C. H. Papadimitriou.  The complexity of computing a Nash equilibrium.  SIAM Journal on Computing, 39(1):195–259, 2009.
  * [67] F. Facchinei and J.-S. Pang.  Finite-dimensional variational inequalities and complementarity problems.  Springer Science & Business Media, 2007.
  * [68] J. S. Pang and M. Razaviyayn.  A unified distributed algorithm for non-cooperative games with non-convex and non-differentiable objectives.  In Big Data over Networks, chapter 4, pages 101–134. Cambridge University Press, 2016.
  * [69] C. D. Dang and G. Lan.  On the convergence properties of non-Euclidean extragradient methods for variational inequalities with generalized monotone operators.  Computational Optimization and applications, 60(2):277–310, 2015\.
  * [70] P. Mertikopoulos, H. Zenati, B. Lecouat, C.-S. Foo, V. Chandrasekhar, and G. Piliouras.  Mirror descent in saddle-point problems: Going the extra (gradient) mile.  arXiv preprint arXiv:1807.02629, 2018.
  * [71] H. Rafique, M. Liu, Q. Lin, and T. Yang.  Weakly-convex–concave min–max optimization: provable algorithms and applications in machine learning.  Optimization Methods and Software, pages 1–35, 2021.
  * [72] C. Song, Z. Zhou, Y. Zhou, Y. Jiang, and Y. Ma.  Optimistic dual extrapolation for coherent non-monotone variational inequalities.  arXiv preprint arXiv:2103.04410, 2021.
  * [73] J. Diakonikolas, C. Daskalakis, and M. I. Jordan.  Efficient methods for structured nonconvex-nonconcave min-max optimization.  In Proceedings of the 24th International Conference on Artificial Intelligence and Statistics, volume 130, pages 2746–2754. PMLR, 2021\.
  * [74] S. Lee and D. Kim.  Semi-anchored multi-step gradient descent ascent method for structured nonconvex-nonconcave composite minimax problems.  arXiv preprint arXiv:2105.15042, 2021.
  * [75] W. Liu, A. Mokhtari, A. Ozdaglar, S. Pattathil, Z. Shen, and N. Zheng.  A decentralized proximal point-type method for saddle point problems.  arXiv preprint arXiv:1910.14380, 2019.
  * [76] B. Barazandeh, D. A. Tarzanagh, and G. Michailidis.  Solving a class of non-convex min-max games using adaptive momentum methods.  In IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 3625–3629. IEEE, 2021.
  * [77] B. Grimmer, H. Lu, P. Worah, and V. Mirrokni.  The landscape of nonconvex-nonconcave minimax optimization.  arXiv preprint arXiv:2006.08667, 2020.
  * [78] T. Fiez, C. Jin, P. Netrapalli, and L. J. Ratliff.  Minimax optimization with smooth algorithmic adversaries.  arXiv preprint arXiv:2106.01488, 2021.
  * [79] O. Mangoubi and N. K. Vishnoi.  Greedy adversarial equilibrium: an efficient alternative to nonconvex-nonconcave min-max optimization.  In Proceedings of the 53rd Annual ACM SIGACT Symposium on Theory of Computing, pages 896–909, 2021.
  * [80] V. Keswani, O. Mangoubi, S. Sachdeva, and N. K. Vishnoi.  A convergent and dimension-independent first-order algorithm for min-max optimization.  arXiv preprint arXiv:2006.12376, 2020.
  * [81] Y. Wang, G. Zhang, and J. Ba.  On solving minimax optimization locally: A follow-the-ridge approach.  arXiv preprint arXiv:1910.07512, 2019.
  * [82] C. Domingo-Enrich, S. Jelassi, A. Mensch, G. Rotskoff, and J. Bruna.  A mean-field analysis of two-player zero-sum games.  arXiv preprint arXiv:2002.06277, 2020.
  * [83] E. V. Mazumdar, M. I. Jordan, and S. S. Sastry.  On finding local Nash equilibria (and only local Nash equilibria) in zero-sum games.  arXiv preprint arXiv:1901.00838, 2019.
  * [84] T. Fiez and L. Ratliff.  Gradient descent-ascent provably converges to strict local minmax equilibria with a finite timescale separation.  arXiv preprint arXiv:2009.14820, 2020.
  * [85] Y.-P. Hsieh, P. Mertikopoulos, and V. Cevher.  The limits of min-max optimization algorithms: convergence to spurious non-critical sets.  arXiv preprint arXiv:2006.09065, 2020.
  * [86] M. Heusel, H. Ramsauer, T. Unterthiner, B. Nessler, and S. Hochreiter.  GANs trained by a two time-scale update rule converge to a local Nash equilibrium.  In Advances in Neural Information Processing Systems, pages 6626–6637, 2017.
  * [87] L. Mescheder, S. Nowozin, and A. Geiger.  The numerics of GANs.  In Advances in Neural Information Processing Systems, pages 1825–1835, 2017.
  * [88] T. Liang and J. Stokes.  Interaction matters: A note on non-asymptotic local convergence of generative adversarial networks.  In The 22nd International Conference on Artificial Intelligence and Statistics, pages 907–915, 2019.
  * [89] F. Farnia and A. Ozdaglar.  GANs may have no Nash equilibria.  arXiv preprint arXiv:2002.09124, 2020.
  * [90] P. Foret, A. Kleiner, H. Mobahi, and B. Neyshabur.  Sharpness-aware minimization for efficiently improving generalization.  In International Conference on Learning Representations, 2020.
  * [91] D. Davis and D. Drusvyatskiy.  Stochastic model-based minimization of weakly convex functions.  SIAM Journal on Optimization, 29(1):207–239, 2019.
  * [92] Y. Carmon and J. C. Duchi.  First-order methods for nonconvex quadratic minimization.  SIAM Review, 62(2):395–436, 2020.
  * [93] H. Karimi, J. Nutini, and M. Schmidt.  Linear convergence of gradient and proximal-gradient methods under the Polyak-Łojasiewicz condition.  In Joint European Conference on Machine Learning and Knowledge Discovery in Databases, pages 795–811. Springer, 2016.
  * [94] Y. Nesterov.  Random walk in a simplex and quadratic optimization over convex polytopes.  Technical report, CORE, 2003.
  * [95] W. Kong and R. D. C. Monteiro.  An accelerated inexact proximal point method for solving nonconvex-concave min-max problems.  arXiv:1905.13433, 2019.
  * [96] T. Lin, C. Jin, and M. I. Jordan.  On gradient descent ascent for nonconvex-concave minimax problems.  arXiv preprint arXiv:1906.00331, 2019.
  * [97] B. Barazandeh and M. Razaviyayn.  Solving non-convex non-differentiable min-max games using proximal gradient method.  In Proceedings of the 45th International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 3162–3166. IEEE, 2020\.
  * [98] N. I. M. Gould, S. Lucidi, M. Roma, and P. L. Toint.  Solving the trust-region subproblem using the Lanczos method.  SIAM Journal on Optimization, 9(2):504–525, 1999.
  * [99] R. T. Rockafellar and R. J.-B. Wets.  Variational analysis, volume 317.  Springer Science & Business Media, 2009.
  * [100] A. Alexandrov.  An application of the theorem on the invariance of domains to existence proofs.  Izvestiya Rossiiskoi Akademii Nauk. Seriya Matematicheskaya, 3(3):243–256, 1939.
  * [101] R. T. Rockafellar.  Second-order convex analysis.  Journal of Nonlinear Convex Analysis, 1(1-16):84, 1999.
  * [102] D. P. Bertsekas.  Nonlinear programming.  Journal of the Operational Research Society, 48(3):334–334, 1997\.
  * [103] J. Cullum and W. E. Donath.  A block Lanczos algorithm for computing the $q$ algebraically largest eigenvalues and a corresponding eigenspace of large, sparse, real symmetric matrices.  In 1974 IEEE Conference on Decision and Control, pages 505–509. IEEE, 1974.
  * [104] G. H. Golub and R. Underwood.  The block Lanczos method for computing eigenvalues.  In Mathematical software, pages 361–377. Elsevier, 1977.
  * [105] A. R. Conn, N. I. M. Gould, and P. L. Toint.  Trust-region methods.  SIAM, 2000.