Extragradient with Positive Momentum is Optimal
for Games with Cross-Shaped Jacobian Spectrum

Junhyung Lyle Kim⋆ Gauthier Gidel Anastasios Kyrillidis Fabian Pedregosa

Rice University Mila & Université de Montréal Canada CIFAR AI Chair Rice
University Google Research

###### Abstract

The extragradient method has recently gained increasing attention, due to its
convergence behavior on smooth games. In $n$-player differentiable games, the
eigenvalues of the Jacobian of the vector field are distributed on the complex
plane, exhibiting more convoluted dynamics compared to classical (i.e., single
player) minimization. In this work, we take a polynomial-based analysis of the
extragradient with momentum for optimizing games with _cross-shaped_ Jacobian
spectrum on the complex plane. We show two results. First, based on the
hyperparameter setup, the extragradient with momentum exhibits three different
modes of convergence: when the eigenvalues are distributed $i)$ on the real
line, $ii)$ both on the real line along with complex conjugates, and $iii)$
only as complex conjugates. Then, we focus on the case $ii)$, i.e., when the
eigenvalues of the Jacobian have _cross-shaped_ structure, as observed in
training generative adversarial networks. For this problem class, we derive
the optimal hyperparameters of the momentum extragradient method, and show
that it achieves an accelerated convergence rate.

##  1 INTRODUCTION

While most of the machine learning models are formulated as the minimization
of a single loss involving data, a growing number of works adopt game
formulations that involve multiple players and objectives. For example,
generative adversarial networks (GANs) [14], actor-critic models [37], and
automatic curricula [43] all can be formulated as a two-player game, to name a
few. Naturally, there is increasing interest in the theoretical understanding
of differentiable games.

Due to the interaction of multiple players and objectives, there are issues in
optimizing games that are not present in minimizing a single objective. Most
notably, the eigenvalues of the game Jacobian lie on the complex plane, unlike
single-objective minimization, where the eigenvalues of the Hessian lie on the
real line. Further, even for simple bilinear games, standard algorithms—like
the gradient method—fail to converge due to the presence of rotational
dynamics [25, 4, 13, 5].

On the other hand, the extragradient method (EG), originally introduced in
[19] for saddle-point problems, was proven to converge for bilinear games
[44]. As such, many recent works theoretically analyze EG from different
angles, such as variational inequality [12, 15], stochastic [27, 22], and
distributed [24, 7] perspectives.

Most existing works, including some of the aforementioned ones, analyze EG by
assuming some structure on the objectives, such as (strong) monotonicity or
Lipschitzness [42, 44, 39]. Such assumptions, in the context of differentiable
games, subsequently confine the distribution of the eigenvalues of the game
Jacobian. For instance, the strong monotonicity and Lipschitzness assumption
translates to the game Jacobians with their eigenvalues in the intersection
between a circle and a half plane [3]. However, such assumptions rarely hold
in applications such as training GANs [14], which is often an empirical
motivation to study differentiable games and convergent algorithms therein,
including EG. Thus, a significant gap between theory and practice has emerged.

In this work, we take a different approach, and analyze the momentum
extragradient method (EGM) through the lens of (residual) polynomials,
inspired by recent works that take a similar approach to analyze (single-
objective) optimization algorithms [20, 36, 40, 17]. Such approach enables us
to consider a setting with a linear vector field where the distribution of the
eigenvalues of the game Jacobian has a _cross-shape_ , as illustrated in
Figure 1, and empirically observed in [5] in training GANs.

Our contributions can be summarized as follows:

  * •

We derive the residual polynomials of EGM, and identify its three different
modes of convergence depending on the hyperparameter setup, that can be
applied to different eigenvalue structures of the game Jacobian (c.f., Theorem
2).

  * •

We then focus on the case where the eigenvalues of the game Jacobian have a
_cross-shaped_ structure, and derive the (asymptotic) accelerated convergence
rate of EGM, as well as its optimal hyperparameters for the considered class
of games (c.f., Theorem 3 and Corollary 1).

  * •

We compare the convergence rate with other first order methods including the
gradient method (GD), the momentum gradient method (GDM), and the
extragradient method (EG). For the class of games we consider with the cross-
shaped eigenvalue structure, we show that none of the aformentioned methods
(asymptotically) converge at an accelerated rate (c.f., Remarks 2 and 3).

###  1.1 Problem Setup

Following [21, 4], we define the $n$-player differentiable game as a family of
twice continuously differentiable losses
$\ell_{i}:\mathbb{R}^{d_{i}}\rightarrow\mathbb{R}$, for $i=1,\dots,n.$ The
player $i$ controls the parameter $w^{(i)}\in\mathbb{R}^{d_{i}}$. We denote
the concatenated parameters by $w=[w^{(1)},\dots,w^{(n)}]\in\mathbb{R}^{d}$,
where $d=\sum_{i=1}^{n}d_{i}.$

For this problem, a Nash equilibrium satisfies:

$w^{(i)^{\star}}\in\operatorname*{arg\,min}_{w^{(i)}\in\mathbb{R}^{d_{i}}}\ell_{i}\big{(}w^{(i)},w^{(\neg
i)^{\star}}\big{)}\quad\forall i\in\\{1,\dots,n\\},$

where the notation $\cdot^{(\neg i)}$ denotes all indices except for $i$. We
also define the vector field of the game $v$ as the concatenation of the
individual gradients:

$v(w)=[\nabla_{w^{(1)}}\ell_{1}(w),\dots,\nabla_{w^{(n)}}\ell_{n}(w)]^{\top},$

and its associated Jacobian $\nabla v$ as:

$\nabla
v(w)=\begin{bmatrix}\nabla^{2}_{w^{(1)}}\ell_{1}(w)&\dots&\nabla_{w^{(n)}}\nabla_{w^{(1)}}\ell_{1}(w)\\\
\vdots&&\vdots\\\
\nabla_{w^{(1)}}\nabla_{w^{(n)}}\ell_{n}(w)&\dots&\nabla^{2}_{w^{(n)}}\ell_{1}(w)\end{bmatrix}.$

Unfortunately, finding Nash equilibria for general games is impractical to
solve [41, 21].111Finding Nash equilibira can be refomulated as a nonlinear
complementarity problem, which is PPAD hard [9, 21]. Instead, we focus on
finding a stationary point of the vector field $v$, since a Nash equilibrium
is always a stationary point of the gradient dynamics. That is, we want to
solve:

$\text{Find}\quad w^{\star}\in\mathbb{R}^{d}\quad\text{such that}\quad
v(w^{\star})=0.$ (1)

Notation. We denote the spectrum of a matrix $M$ by $\text{Sp}(M),$ and its
spectral radius by $\rho(M):=\max\\{|\lambda|:\lambda\in\text{Sp}(M)\\}$.
$M\succ 0$ denotes that $M$ is a positive-definite matrix. $\mathfrak{R}(z)$
and $\mathfrak{I}(z)$ respectively denote the real and the imaginary part of a
complex number $z.$

###  1.2 Related Work

The above game formulation is strictly more general than minimization, and
includes many problems such as saddle-point problems [33], bilinear games [23,
28], and variational inequality problems [44, 11].

For minimization, the vector field $v$ and the Jacobian $\nabla v$ are,
respectively, the gradient and the Hessian of the objective function. The
Hessian is a symmetric matrix, so its eigenvalues lie on the real line; the
ratio of the largest and the smallest eigenvalues is called the condition
number. It is now a common wisdom in numerical linear algebra and optimization
theory how the notion of condition number characterizes the difficulty of a
(strongly convex) minimization problem [31]. In conjunction with provable
lower bounds within various function classes [30], an algorithm is considered
“optimal” for a class of functions when its attainable upper bound on the
convergence rate matches the lower bound of that function class [32, 31].

Differentiable games, on the other hand, have much more convoluted dynamics,
as the eigenvalues of the Jacobian at the solution are distributed over the
complex plane [26, 4]. To address this, recent line of works proposes new
algorithms with better upper bounds on the convergence rate [29, 34, 26, 10,
12, 13, 28]. Vast majority of the aforementioned literature proves convergence
based on two approaches: either bounding $i)$ the magnitude, or $ii)$ the real
part of the eigenvalues of the (submatrices of the) Jacobian. In this paper we
make instead a more fine-grained assumption on these eigenvalues, to obtain
faster algorithms and convergence rates.

To understand the difficulty of a differentiable game, [3] proposed a
geometric interpretation of the conditioning of a game via the notion of
spectral shape, defined as the set containing all eigenvalues of the Jacobian.
Specifically, let $\mathcal{A}_{\mathcal{K}}$ be the set of matrices $A$ whose
spectrum, denoted by $\text{Sp}(A)$, belong to a set $\mathcal{K}$ on the
complex plane with positive real part:

$\mathcal{A}_{\mathcal{K}}:=\\{A\in\mathbb{R}^{d\times
d}:\text{Sp}(A)\subset\mathcal{K}\subset\mathbb{C}_{+}\\}.$

When $\mathcal{K}$ is “simple” (e.g., a real line segment or a disc on the
complex plane), [3] characterized the lower bound and the optimality of some
first-order methods [30, 1] (see also Section 3.2).

###  1.3 Games with Cross-Shaped Jacobian Spectrum

We consider a more specific class of games, where the Jacobian $\nabla v$ is a
liner operator and its eigenvalues are contained in a (shifted) cross-shape on
the complex plane. That is, we consider the following spectrum model:

$\begin{split}&\text{Sp}(\nabla v)\in\mathcal{S}^{\star}\\\
&=[\mu,L]\cup\\{a+bi\in\mathbb{C}:a=c^{\prime}>0,~{}b\in[-c,c]\\}.\end{split}$
(2)

The first set $[\mu,L]$ denotes a segment on the real line. This interval also
arises when minimizing $\mu$-strongly convex and $L$-smooth functions. The
second set has a fixed real component ($c^{\prime}>0$), along with imaginary
components symmetric across the real line (i.e., complex conjugates), as the
Jacobian is real. We illustrate (2) in Figure 1.

Figure 1: Illustration of the cross-shaped spectrum model in (2). A special
case considered in Section 3.2 is when $c=\frac{L-\mu}{2},$ i.e., when the
length of the two intervals coincide.

This is a strict generalization of the purely imaginary interval $[-bi,bi]$
commonly considered in the bilinear games literature [23, 3, 28]. While many
recent papers on bilinear games cite GANs [14] as a motivation, the work in
[5, Figure 4] empirically shows that the spectrum of GANs is not contained in
the imaginary axis. Instead, we note that the cross-shaped spectrum model
above would include some of the observed GAN spectra.

The cross-shaped spectrum in (2) is also observed in a class of quadratic
games. Specifically, consider the following two player quadratic game, where
$x\in\mathbb{R}^{d_{1}}$ and $y\in\mathbb{R}^{d_{2}}$ are the parameters
controlled by each player, whose loss function respectively is:

$\displaystyle\ell_{1}(x,y)$
$\displaystyle=\frac{1}{2}x^{\top}S_{1}x+x^{\top}M_{12}y+x^{\top}b_{1}\quad\text{and}$
(3) $\displaystyle\ell_{2}(x,y)$
$\displaystyle=\frac{1}{2}y^{\top}S_{2}y+y^{\top}M_{21}x+y^{\top}b_{2},$ (4)

where $S_{1},S_{2}\succ 0.$ Then, the vector field can be written as:

$\displaystyle v(x,y)=\begin{bmatrix}S_{1}x+M_{12}y+b_{1}\\\
M_{21}x+S_{2}y+b_{2}\end{bmatrix}=Aw+b,~{}~{}\text{where}$ $\displaystyle
A=\begin{bmatrix}S_{1}&M_{12}\\\
M_{21}&S_{2}\end{bmatrix},~{}~{}w=\begin{bmatrix}x\\\
y\end{bmatrix},~{}~{}\text{and}~{}~{}b=\begin{bmatrix}b_{1}\\\
b_{2}\end{bmatrix}.$ (5)

If $S_{1}=S_{2}=0$ and $M_{12}=-M_{21}^{\top},$ the game Jacobian $\nabla v=A$
has only purely imaginary eigenvalues [3, Lemma 7], recovering bilinear games.

Further, when $M_{12}=-M_{21}^{\top},$ and they share a common basis with
$S_{1}$ and $S_{2}$ as in the below proposition, then $A$ has a cross-shaped
spectrum. We utilize this proposition in designing our experiments in Section
4.

######  Proposition 1.

Let $A$ be a matrix of the form $\begin{bmatrix}S_{1}&B\\\
-B^{\top}&S_{2}\end{bmatrix}$, where $S_{1},S_{2}\succ 0$. Without loss of
generality, assume that $dim(S_{1})>dim(S_{2})=d$. Then, $\text{Sp}(A)$ has a
cross-shape if there exists orthonormal matrices $U,V$ and diagonal matrices
$D_{1},D_{2}$ such that

$\displaystyle S_{1}$
$\displaystyle=U\text{diag}(a,\ldots,a,D_{1})U^{\top},~{}S_{2}=\text{diag}(a,\ldots,a),~{}\text{and}$
$\displaystyle B$ $\displaystyle=UD_{2}V^{\top}.$

Note that all proofs can be found in the supplementary material. To solve (1)
for games with Jacobian eigenvalue structure in (2), we use the extragradient
with momentum, which we expound in the next section.

##  2 ALGORITHM

###  2.1 Momentum Extragradient via Chebyshev Polynomials

To solve (1) for games with Jacobian eigenvalue structure in (2), we focus on
the extragradient method with momentum (EGM):

$w_{t+1}=w_{t}-{\color[rgb]{0,0.75,0.16}h}v(w_{t}-{\color[rgb]{0.6,0.2,0.8}\gamma}v(w_{t}))+{\color[rgb]{.75,0,.25}m}(w_{t}-w_{t-1}),$
(6)

where ${\color[rgb]{0,0.75,0.16}h}$ is the step size,
${\color[rgb]{0.6,0.2,0.8}\gamma}$ is the extrapolation step size, and
${\color[rgb]{.75,0,.25}m}$ is the momentum parameter. The extragradient (EG)
method (without momentum) was first proposed in [19] for saddle point
problems, and recently received much attention due to its convergence behavior
for some class of differentiable games such as bilinear games, for which the
standard gradient method diverges [13, 3, 2]. For completeness, we remind the
gradient method (GD) and the momentum gradient method (GDM):

$\displaystyle\text{(GD)}~{}~{}w_{t+1}$
$\displaystyle=w_{t}-{\color[rgb]{0,0.75,0.16}h}v(w_{t}),\quad\text{and}\quad$
$\displaystyle\text{(GDM)}~{}~{}w_{t+1}$
$\displaystyle=w_{t}-{\color[rgb]{0,0.75,0.16}h}v(w_{t})+{\color[rgb]{.75,0,.25}m}(w_{t}-w_{t-1}).$
(7)

We take a polynomial-based approach to analyze EGM. Such analysis dates back
to the development of the conjugate gradient method [18], and is still
actively used, for instance, to derive lower bounds [1], develop accelerated
decentralized algorithms [6], and analyze average-case performance [36].

On that end, we use the following lemma [8], which elucidates the intimate
connection between first-order methods222We consider first-order methods as
those that only use first-order information in each iteration. and (residual)
polynomials, when the vector field satisfies $v(w)=Aw+b$:

######  Lemma 1 ([8]).

There exists a real polynomial $p_{t}$ of degree at most $t$ satisfying

$w_{t}-w^{\star}=p_{t}(A)(w_{0}-w^{\star})\,,$ (8)

where $p_{t}(0)=1$, and $v(w^{\star})=Aw^{\star}+b=0$.

By taking norms, (8) further implies the following worst-case convergence
rate:

$\displaystyle\|w_{t}-w^{\star}\|$
$\displaystyle=\|p_{t}(A)(w_{0}-w^{\star})\|$
$\displaystyle\leqslant\|p_{t}(A)\|\cdot\|w_{0}-w^{\star}\|.$ (9)

As EGM in (6) is a first-order method [1, 3], we can study the associated
residual polynomials, as we summarize in the following theorem:

######  Theorem 1 (Residual polynomials of EGM).

Consider the momentum extragradient method in (6). Its associated residual
polynomials are:

$\begin{split}\tilde{P}_{0}(\lambda)&=1,~{}~{}\tilde{P}_{1}(\lambda)=1-\tfrac{{\color[rgb]{0,0.75,0.16}h}\lambda(1-{\color[rgb]{0.6,0.2,0.8}\gamma}\lambda)}{1+{\color[rgb]{.75,0,.25}m}},~{}~{}\text{and}~{}\\\
\tilde{P}_{t+1}(\lambda)&=(1+{\color[rgb]{.75,0,.25}m}-{\color[rgb]{0,0.75,0.16}h}\lambda(1-{\color[rgb]{0.6,0.2,0.8}\gamma}\lambda))\tilde{P}_{t}(\lambda)-{\color[rgb]{.75,0,.25}m}\tilde{P}_{t-1}(\lambda).\end{split}$

Further, let $T_{t}(\cdot)$ and $U_{t}(\cdot)$ be the Chebyshev polynomials of
the first and the second kind, respectively. Then, the above expressions
simplify to:

$\displaystyle
P_{t}(\lambda)={\color[rgb]{.75,0,.25}m}^{t/2}\left({\tfrac{2{\color[rgb]{.75,0,.25}m}}{1+{\color[rgb]{.75,0,.25}m}}}T_{t}(\sigma(\lambda))+{\tfrac{1-{\color[rgb]{.75,0,.25}m}}{1+{\color[rgb]{.75,0,.25}m}}}U_{t}(\sigma(\lambda))\right)$
(10) $\displaystyle\text{with
}\sigma(\lambda)=\tfrac{1+{\color[rgb]{.75,0,.25}m}-{\color[rgb]{0,0.75,0.16}h}\lambda(1-{\color[rgb]{0.6,0.2,0.8}\gamma}\lambda)}{2\sqrt{{\color[rgb]{.75,0,.25}m}}},$
(11)

where we refer to the term $\sigma(\lambda)$ as the link function.

As can be seen in (11), the link function for EGM is quadratic with respect to
$\lambda$. As we show below, GDM admits an almost identical expression, but
with a different link function. For completeness, we write below the residual
polynomials of GDM, expressed in terms of the Chebyshev polynomials [35]:

$\displaystyle P_{t}^{\text{GDM}}(\lambda)$
$\displaystyle={\color[rgb]{.75,0,.25}m}^{t/2}\left({\tfrac{2{\color[rgb]{.75,0,.25}m}}{1+{\color[rgb]{.75,0,.25}m}}}T_{t}(\xi(\lambda)){+}{\tfrac{1-{\color[rgb]{.75,0,.25}m}}{1+{\color[rgb]{.75,0,.25}m}}}U_{t}(\xi(\lambda))\right)$
(12) $\displaystyle\text{ with }\xi(\lambda)$
$\displaystyle=\tfrac{1+{\color[rgb]{.75,0,.25}m}-{\color[rgb]{0,0.75,0.16}h}\lambda}{2\sqrt{{\color[rgb]{.75,0,.25}m}}}.$
(13)

Notice that the Chebyshev polynomial-based expression of EGM in (10) and that
of GDM in (12) are identical except for the link functions $\sigma(\lambda)$
and $\xi(\lambda)$, which enter as arguments in $T_{t}(\cdot)$ and
$U_{t}(\cdot)$. This difference in link functions is crucial, as the Chebyshev
polynomials behave very differently based on the domain, as shown in the
following lemma.

######  Lemma 2 ([16]).

Let $z$ be a complex number. The sequence

$\displaystyle\left(\left|\tfrac{2{\color[rgb]{.75,0,.25}m}}{1+{\color[rgb]{.75,0,.25}m}}T_{t}(z)+\tfrac{1-{\color[rgb]{.75,0,.25}m}}{1+{\color[rgb]{.75,0,.25}m}}U_{t}(z)\right|\right)_{t\geqslant
0}$

grows exponentially in $t$ for $z\notin[-1,1]$, while in that interval, the
following bounds hold:

$|T_{t}(z)|\leqslant 1\quad\text{and}\quad|U_{t}(z)|\leqslant t+1.\\\ $ (14)

Therefore, to study the convergence behavior of EGM, we are interested in the
case where the set of step sizes and the momentum parameters lead to
$|\sigma(\lambda)|\leqslant 1$, so that we can use the bounds in (14). We call
the set of hyperparameters such that the image of the link function lies in
the interval $[-1,1]$ as the robust region, i.e., $\sigma^{-1}([-1,1]).$

###  2.2 Three Modes of the Momentum Extragradient

Within the robust region of EGM, we can compute its worst-case convergence
rate as follows:

$\displaystyle r_{t}$
$\displaystyle:=\max_{\lambda\in\mathcal{S}^{\star}}|P_{t}(\lambda)|$
$\displaystyle=\max_{\lambda\in\mathcal{S}^{\star}}\bigg{|}{\color[rgb]{.75,0,.25}m}^{t/2}\bigg{(}\tfrac{2m}{1+{\color[rgb]{.75,0,.25}m}}T_{t}(\sigma(\lambda))+\tfrac{1-{\color[rgb]{.75,0,.25}m}}{1+{\color[rgb]{.75,0,.25}m}}U_{t}(\sigma(\lambda))\bigg{)}\bigg{|}$
$\displaystyle\leqslant{\color[rgb]{.75,0,.25}m}^{t/2}\bigg{(}\tfrac{2m}{1+{\color[rgb]{.75,0,.25}m}}\max_{\lambda\in\mathcal{S}^{\star}}|T_{t}(\sigma(\lambda))|+\tfrac{1-{\color[rgb]{.75,0,.25}m}}{1+{\color[rgb]{.75,0,.25}m}}\max_{\lambda\in\mathcal{S}^{\star}}|U_{t}(\sigma(\lambda))|\bigg{)}$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:chebyshev-
complex}}}{{\leqslant}}{\color[rgb]{.75,0,.25}m}^{t/2}\bigg{(}\tfrac{2m}{1+{\color[rgb]{.75,0,.25}m}}+\tfrac{1-{\color[rgb]{.75,0,.25}m}}{1+{\color[rgb]{.75,0,.25}m}}(t+1)\bigg{)}\leqslant{\color[rgb]{.75,0,.25}m}^{t/2}(t+2),$
(15)

where the last inequality uses ${\color[rgb]{.75,0,.25}m}\in(0,1).$ Note that
$r_{t}$ denotes the worst-case convergence rate by definition. Indeed, in
conjunction with (9), it follows
$\frac{\|w_{t}-w^{\star}\|}{\|w_{0}-w^{\star}\|}\leqslant r_{t}$. To get the
asymptotic rate, we take the limit supremum of the $2t$-th root of $r_{t}$:

$\displaystyle\limsup_{t\rightarrow\infty}\sqrt[2t]{r_{t}}=\limsup_{t\rightarrow\infty}\Big{(}{\color[rgb]{.75,0,.25}m}^{t/2}(t+2)\Big{)}^{\frac{1}{2t}}=\sqrt[4]{{\color[rgb]{.75,0,.25}m}}.$
(16)

The reason why we take the $2t$-th root of $r_{t}$ is to normalize by the
number of vector field computations, as we compare in Section 3 the asymptotic
rate of EGM in (16) with other gradient methods that use a single vector field
computation in the recurrences, such as GD and GDM.

Recall that the Chebyshev polynomial expressions of EGD in (10) and that of
GDM333Asymptotically, GDM enjoys $\sqrt{{\color[rgb]{.75,0,.25}m}}$
convergence rate instead of the $\sqrt[4]{{\color[rgb]{.75,0,.25}m}}$ of EGM,
as it uses a single vector field computation per iteration instead of the two.
are identical except for the link functions. Hence, the convergence rate in
(15) applies to both EGM and GDM, as long as the link functions
$|\sigma(\lambda)|$ and $|\xi(\lambda)|$ are bounded by $1$.

Yet, due to the difference in link functions, the robust region of EGM and GDM
are drastically different. Specifically, due to the quadratic link function of
EGM, the robust region is more flexible; in particular, this link function can
still satisfy the upper bounds (14) in Lemma 2 for a subset of complex values,
resulting in the convergence rate of (15). Hence, compared to GDM whose link
function is linear in $\lambda$ and so can only satisfy (14) for a subset of
real values, EGM enjoys the convergence rate in (15) for a much wider set of
$\lambda$ values.

The robust region of EGM can be described with the four extreme points below:

$\displaystyle\sigma^{-1}(-1)$
$\displaystyle=\tfrac{1}{2{\color[rgb]{0.6,0.2,0.8}\gamma}}\pm\sqrt{\tfrac{1}{4{\color[rgb]{0.6,0.2,0.8}\gamma}^{2}}-\tfrac{(1+\sqrt{{\color[rgb]{.75,0,.25}m}})^{2}}{{\color[rgb]{0,0.75,0.16}h}{\color[rgb]{0.6,0.2,0.8}\gamma}}},\quad\text{and}\quad$
(17) $\displaystyle\sigma^{-1}(1)$
$\displaystyle=\tfrac{1}{2{\color[rgb]{0.6,0.2,0.8}\gamma}}\pm\sqrt{\tfrac{1}{4{\color[rgb]{0.6,0.2,0.8}\gamma}^{2}}-\tfrac{(1-\sqrt{{\color[rgb]{.75,0,.25}m}})^{2}}{{\color[rgb]{0,0.75,0.16}h}{\color[rgb]{0.6,0.2,0.8}\gamma}}}.$
(18)

Depending on the choice of hyperparameters, the above four extreme points can
be distributed differently on three different modes, which we summarize in the
next theorem.

######  Theorem 2.

Consider the extragradient with momentum in (6), expressed with the Chebyshev
polynomials as in (10). Then, the robust region summarized in (17) and (18)
have the following three modes:

  1. 1.

If
$\frac{{\color[rgb]{0,0.75,0.16}h}}{4{\color[rgb]{0.6,0.2,0.8}\gamma}}\geqslant(1+\sqrt{{\color[rgb]{.75,0,.25}m}})^{2}$,
then $\sigma^{-1}(-1)$ and $\sigma^{-1}(1)$ are all real numbers;

  2. 2.

If
$(1-\sqrt{{\color[rgb]{.75,0,.25}m}})^{2}\leqslant\frac{{\color[rgb]{0,0.75,0.16}h}}{4{\color[rgb]{0.6,0.2,0.8}\gamma}}<(1+\sqrt{{\color[rgb]{.75,0,.25}m}})^{2}$,
then $\sigma^{-1}(-1)$ are complex, and $\sigma^{-1}(1)$ are real;

  3. 3.

If
$(1-\sqrt{{\color[rgb]{.75,0,.25}m}})^{2}>\frac{{\color[rgb]{0,0.75,0.16}h}}{4{\color[rgb]{0.6,0.2,0.8}\gamma}},$
then $\sigma^{-1}(-1)$ and $\sigma^{-1}(1)$ are all complex numbers.

We empirically verify Theorem 2 in Figure 2. Note that the hyperparameters are
set up so that each case in Theorem 2 is covered. In the remainder of the
paper, we specifically focus on the second case, which can be modeled via the
Jacobian spectrum model in (2).

Figure 2: Red points represent the input arguments to the link function
$\sigma(\lambda)$ in (11) that outputs the range in $[-1,1]$. Black solid
lines are the real and the imaginary axes. Left: Case 1; Middle: Case 2;
Right: Case 3 of Theorem 2.

##  3 OPTIMAL PARAMETERS AND CONVERGENCE RATES

###  3.1 Optimal Parameters for Cross-Shaped Games

We now focus on solving (1) for games with cross-shaped Jacobian eigenvalue
structure as in (2). Therefore, in the remainder of the paper, we focus on the
second case of Theorem 2, which is illustrated in the middle panel of Figure
2. As mentioned in Section 1.3, the cross-shaped eigenvalue structure of game
Jacobians have been empirically observed in [5, Figure 4] while training GANs
in various settings. Therefore, we analyze the more realistic and strictly
more general setting than the purely imaginary interval $[-bi,bi]$ commonly
considered in the bilinear game literature [23, 3, 28].

We now connect the cross-shaped spectrum model in (2) and the robust region of
EGM described with four extra points in (17) and (18). As we focus on the
second case of Theorem 2, the hyperparameters satisfy the condition
$(1-\sqrt{{\color[rgb]{.75,0,.25}m}})^{2}\leqslant\tfrac{{\color[rgb]{0,0.75,0.16}h}}{4{\color[rgb]{0.6,0.2,0.8}\gamma}}<(1+\sqrt{{\color[rgb]{.75,0,.25}m}})^{2}.$
In this case, we can write the robust region of EGM as:

$\displaystyle\sigma^{-1}([-1,1])$
$\displaystyle=\bigg{[}\tfrac{1}{2{\color[rgb]{0.6,0.2,0.8}\gamma}}{-}\sqrt{\tfrac{1}{4{\color[rgb]{0.6,0.2,0.8}\gamma}^{2}}{-}\tfrac{(1-\sqrt{{\color[rgb]{.75,0,.25}m}})^{2}}{{\color[rgb]{0,0.75,0.16}h}{\color[rgb]{0.6,0.2,0.8}\gamma}}},\tfrac{1}{2{\color[rgb]{0.6,0.2,0.8}\gamma}}{+}\sqrt{\tfrac{1}{4{\color[rgb]{0.6,0.2,0.8}\gamma}^{2}}{-}\tfrac{(1-\sqrt{{\color[rgb]{.75,0,.25}m}})^{2}}{{\color[rgb]{0,0.75,0.16}h}{\color[rgb]{0.6,0.2,0.8}\gamma}}}\bigg{]}\bigcup$
$\displaystyle~{}~{}~{}~{}~{}\bigg{[}\tfrac{1}{2{\color[rgb]{0.6,0.2,0.8}\gamma}}{-}\sqrt{\tfrac{1}{4{\color[rgb]{0.6,0.2,0.8}\gamma}^{2}}{-}\tfrac{(1+\sqrt{{\color[rgb]{.75,0,.25}m}})^{2}}{{\color[rgb]{0,0.75,0.16}h}{\color[rgb]{0.6,0.2,0.8}\gamma}}},\tfrac{1}{2{\color[rgb]{0.6,0.2,0.8}\gamma}}{+}\sqrt{\tfrac{1}{4{\color[rgb]{0.6,0.2,0.8}\gamma}^{2}}{-}\tfrac{(1+\sqrt{{\color[rgb]{.75,0,.25}m}})^{2}}{{\color[rgb]{0,0.75,0.16}h}{\color[rgb]{0.6,0.2,0.8}\gamma}}}\bigg{]}.$
(19)

Here, the first interval lies on $\mathbb{R},$ as the square root term is
real; conversely, in the second interval, the square root term is imaginary,
with the fixed real component: $\frac{1}{2{\color[rgb]{0.6,0.2,0.8}\gamma}}.$

The optimal hyperparameters of EGM in terms of the worst-case asymptotic
convergence rate in (16) for this problem occurs when the robust region in
(3.1) and the spectrum model in (2) coincide. This condition can be summarized
in the following three equalities:

$\begin{split}\tfrac{1}{2{\color[rgb]{0.6,0.2,0.8}\gamma}}-\sqrt{\tfrac{1}{4{\color[rgb]{0.6,0.2,0.8}\gamma}^{2}}-\tfrac{(1-\sqrt{{\color[rgb]{.75,0,.25}m}})^{2}}{{\color[rgb]{0,0.75,0.16}h}{\color[rgb]{0.6,0.2,0.8}\gamma}}}&=\mu,\\\
\tfrac{1}{2{\color[rgb]{0.6,0.2,0.8}\gamma}}+\sqrt{\tfrac{1}{4{\color[rgb]{0.6,0.2,0.8}\gamma}^{2}}-\tfrac{(1-\sqrt{{\color[rgb]{.75,0,.25}m}})^{2}}{{\color[rgb]{0,0.75,0.16}h}{\color[rgb]{0.6,0.2,0.8}\gamma}}}&=L,~{}~{}\text{and}~{}~{}\\\
\sqrt{\tfrac{(1+\sqrt{{\color[rgb]{.75,0,.25}m}})^{2}}{{\color[rgb]{0,0.75,0.16}h}{\color[rgb]{0.6,0.2,0.8}\gamma}}-\tfrac{1}{4{\color[rgb]{0.6,0.2,0.8}\gamma}^{2}}}&=c.\end{split}$

Based on the above equalities, we can compute the optimal hyperparameters, as
summarized in the next theorem:

######  Theorem 3 (Optimal hyperparameters for EGM).

Consider solving (1) for games where the Jacobian has cross-shaped spectrum as
in (2). For this problem, the optimal hyperparameters for the momentum
extragradient method in (6) can be set as follows:

$\begin{split}{\color[rgb]{0,0.75,0.16}h}&=\tfrac{16(\mu+L)}{(\sqrt{4c^{2}+(\mu+L)^{2}}+\sqrt{4\mu
L})^{2}},\quad{\color[rgb]{0.6,0.2,0.8}\gamma}=\tfrac{1}{\mu+L},\quad\text{and}\quad\\\
{\color[rgb]{.75,0,.25}m}&=\left(\tfrac{\sqrt{4c^{2}+(\mu+L)^{2}}-\sqrt{4\mu
L}}{\sqrt{4c^{2}+(\mu+L)^{2}}+\sqrt{4\mu L}}\right)^{2}.\end{split}$

We summarize EGM with the above optimal hyperparameters in Algorithm 1.
Recalling (16), we immediately get the asymptotic convergence rate from
Theorem 3. Further, this formula can be simplified in the ill-conditioned
regime, where the inverse condition number $\tau:=\mu/L\rightarrow 0$:

$\displaystyle\sqrt[4]{{\color[rgb]{.75,0,.25}m}}$
$\displaystyle=\left(\tfrac{\sqrt{4c^{2}+(\mu+L)^{2}}-\sqrt{4\mu
L}}{\sqrt{4c^{2}+(\mu+L)^{2}}+\sqrt{4\mu L}}\right)^{1/2}$
$\displaystyle=\left(\tfrac{\sqrt{(2c/L)^{2}+(1+\tau)^{2}}-2\sqrt{\tau}}{\sqrt{(2c/L)^{2}+(1+\tau)^{2}}+2\sqrt{\tau}}\right)^{1/2}$
$\displaystyle\underset{\tau\rightarrow
0}{=}1-\tfrac{2\sqrt{\tau}}{\sqrt{(2c/L)^{2}+1}}+o(\sqrt{\tau}).$ (20)

From (20), we can see that EGM achieves accelerated convergence rate
$1-O(\sqrt{\mu/L})$, as long as $c=O(L).$

A special case of the spectrum model in (2) is when the length of the first
set and the second equals, i.e., when $c=\tfrac{L-\mu}{2}$. Then, the optimal
hyperparameters of EGM in Theorem 3 further simplify, as in the next
corollary:

######  Corollary 1.

Consider solving (1) for games where the Jacobian has cross-shaped spectrum as
in (2). Further, assume that $c=\tfrac{L-\mu}{2}$. Then, the optimal
hyperparameters for the momentum extragradient method in (6) can be set as
follows:

$\begin{split}{\color[rgb]{0,0.75,0.16}h}&=\tfrac{8(\mu+L)}{(\sqrt{\mu^{2}+L^{2}}+\sqrt{2\mu
L})^{2}},\quad{\color[rgb]{0.6,0.2,0.8}\gamma}=\tfrac{1}{\mu+L},\quad\text{and}\quad\\\
{\color[rgb]{.75,0,.25}m}&=\left(\tfrac{\sqrt{\mu^{2}+L^{2}}-\sqrt{2\mu
L}}{\sqrt{\mu^{2}+L^{2}}+\sqrt{2\mu L}}\right)^{2}.\end{split}$

Again, from Corollary 1, we can immediately get the asymptotic convergence
rate via (16), which can be interpreted further in the ill-conditioned regime
where $\tau\rightarrow 0$:

$\displaystyle\sqrt[4]{{\color[rgb]{.75,0,.25}m}}$
$\displaystyle=\left(\tfrac{\sqrt{\mu^{2}+L^{2}}-\sqrt{2\mu
L}}{\sqrt{\mu^{2}+L^{2}}+\sqrt{2\mu L}}\right)^{1/2}$
$\displaystyle\underset{\tau\rightarrow
0}{=}1-\sqrt{2}\sqrt{\tau}+o(\sqrt{\tau})\approx
1-\sqrt{2}\sqrt{\tfrac{\mu}{L}}.$ (21)

Similarly to (20), in the simplified setting of the Jacobian spectrum model in
(2) with $c=\frac{L-\mu}{2}$, EGM with optimal hyperparameters achieve an
accelerated convergence rate.

######  Remark 1.

Notice that the optimal momentum ${\color[rgb]{.75,0,.25}m}$ in both Theorem 3
and Corollary 1 are positive. This is in contrast with [13] where the gradient
method with negative momentum is studied. This difference elucidates the
different dynamics of how momentum interacts with the gradient and the
extragradient methods.

Input: Initialization $w_{0}$, problem constants $\mu,L,c$.

Set:
${\color[rgb]{0,0.75,0.16}h}=\tfrac{16(\mu+L)}{(\sqrt{4c^{2}+(\mu+L)^{2}}+\sqrt{4\mu
L})^{2}},\quad{\color[rgb]{0.6,0.2,0.8}\gamma}=\tfrac{1}{\mu+L},~{}~{}\text{and}~{}~{}$
${\color[rgb]{.75,0,.25}m}=\left(\tfrac{\sqrt{4c^{2}+(\mu+L)^{2}}-\sqrt{4\mu
L}}{\sqrt{4c^{2}+(\mu+L)^{2}}+\sqrt{4\mu L}}\right)^{2}.$

$w_{1}=w_{0}-\tfrac{{\color[rgb]{0,0.75,0.16}h}}{1+{\color[rgb]{.75,0,.25}m}}v(w_{0}-{\color[rgb]{0.6,0.2,0.8}\gamma}v(w_{0}))$

for  _$t=1,2,\dots$_ do

$w_{t+1}=w_{t}-{\color[rgb]{0,0.75,0.16}h}v(w_{t}-{\color[rgb]{0.6,0.2,0.8}\gamma}v(w_{t}))+{\color[rgb]{.75,0,.25}m}(w_{t}-w_{t-1})$

end for

Algorithm 1 Momentum extragradient with optimal hyperparameters for cross-
shaped spectrum in (2)

###  3.2 Comparison with Other Methods

We now compare the asymptotic rates of EGM we obtained in the previous
subsection with other first order methods, including GD, GDM, and EG. Note
that for simplicity we focus on the setting when $c=\frac{L-\mu}{2}$ in (2),
i.e., when the two intervals in (2) have the same length. The optimal
hyperparameters for this setting was presented in Corollay 1, and the
asymptotic convergence rate was obtained in (21). The analysis of the general
case with any $c$ can be extended with similar approach, which we leave for
future work.

In [2], GD and EG are interpreted as fixed-point iterations, following [38].
In this framework, the iterates of a method are considered to be generated by
the following scheme:

$w_{t+1}=F(w_{t}),\quad\forall t\geqslant 0,$ (22)

where $F:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ is an operator representing
the method. Analyzing the above scheme in general settings is usually
challenging as $F$ can be any nonlinear operator. However, if $F$ is twice
differentiable and $w$ is near a stationary point $w^{\star}$, the analysis
can be simplified by linearizing $F$ at the stationary point, i.e.,

$F(w)\approx F(w^{\star})+\nabla F(w^{\star})(w-w^{\star}).$

Then, for $w_{0}$ in a neighborhood of $w^{\star}$, [38] shows one can obtain
an asymptotic convergence rate of (22) by $\rho(\nabla
F(w^{\star}))\leqslant\rho^{\star}<1,$ meaning that (22) converges linearly to
$w^{\star}$ at the rate $O((\rho^{\star}+\varepsilon)^{t})$ for
$\varepsilon\geqslant 0$. Further, if $F$ is linear, $\varepsilon=0$.

Under this scheme, the fixed point operators
$F_{{\color[rgb]{0,0.75,0.16}h}}^{\text{GD}}$ and
$F_{{\color[rgb]{0,0.75,0.16}h}}^{\text{EG}}$ of GD and EG,444[2] assumes that
EG uses the same step size ${\color[rgb]{0,0.75,0.16}h}$ for both the main and
the extrapolation steps. respectively, can be written as follows:

$\displaystyle\text{(GD)}~{}~{}w_{t+1}$
$\displaystyle=w_{t}-{\color[rgb]{0,0.75,0.16}h}v(w_{t})=F_{{\color[rgb]{0,0.75,0.16}h}}^{\text{GD}}(w_{t}),~{}~{}\text{and}~{}~{}$
(23) $\displaystyle\text{(EG)}~{}~{}w_{t+1}$
$\displaystyle=w_{t}-{\color[rgb]{0,0.75,0.16}h}v(w_{t}-{\color[rgb]{0,0.75,0.16}h}v(w_{t}))=F_{{\color[rgb]{0,0.75,0.16}h}}^{\text{EG}}(w_{t}).$
(24)

The local convergence rate can now be obtained by bounding the spectral radius
of the Jacobian of the operators under certain assumptions. We summarize the
relevant theorems from [13] and [2] below:

######  Theorem 4 ([2]).

Let $w^{\star}$ be a stationary point of $v$. Denote by
$\Delta^{\star}:=\text{Sp}(\nabla v(w^{\star}))$ for notational brevity.
Further, assume the eigenvalues of $\nabla v(w^{\star})$ all have positive
real parts. Then,

  1. 1.

[13] For the gradient method in (23) with step size
${\color[rgb]{0,0.75,0.16}h}=\min_{\lambda\in\Delta^{\star}}\mathfrak{R}(1/\lambda)$,
it satisfies:555Note that the spectral radius $\rho$ is squared, but
asymptotically is almost the same as $\sqrt{1-x}\leqslant 1-x/2$.

$\rho(\nabla
F_{{\color[rgb]{0,0.75,0.16}h}}^{\text{GD}}(w^{\star}))^{2}\leqslant
1-\min_{\lambda\in\Delta^{\star}}\mathfrak{R}\left(\frac{1}{\lambda}\right)\min_{\lambda\in\Delta^{\star}}\mathfrak{R}(\lambda).$
(25)

  2. 2.

For the extragradient method in (24) with step size
${\color[rgb]{0,0.75,0.16}h}=(4\max_{\lambda\in\Delta^{\star}}|\lambda|)^{-1}$,
it satisfies:

$\displaystyle\rho_{\text{EG}}^{2}\leqslant
1-\tfrac{1}{4}\left(\tfrac{\min_{\lambda\in\Delta^{\star}}\mathfrak{R}(\lambda)}{\max_{\lambda\in\Delta^{\star}}|\lambda|}+\tfrac{\min_{\lambda\in\Delta^{\star}}|\lambda|^{2}}{\max_{\lambda\in\Delta^{\star}}|\lambda|^{2}}\right),$
(26)

where we denoted $\rho_{\text{EG}}:=\rho(\nabla
F_{{\color[rgb]{0,0.75,0.16}h}}^{\text{EG}}(w^{\star}))$ for brevity.

Since our spectrum model in (2) satisfies the condition that the eigenvalues
of $\nabla v(w^{\star})$ all have positive real parts, we can obtain the
convergence rate of GD and EG, based on Theorem 4. With our (simplified)
spectrum model with $c=\frac{L-\mu}{2},$ for GD, it follows that
$\min_{\lambda\in\Delta^{\star}}\mathfrak{R}(\lambda)=\mu$, and
$\min_{\lambda\in\Delta^{\star}}\mathfrak{R}(1/\lambda)=1/L$. Similarly, for
EG, it follows that $\max_{\lambda\in\Delta^{\star}}|\lambda|=L,$ and
$\max_{\lambda\in\Delta^{\star}}|\lambda|^{2}=L^{2}.$ Furthermore,
$\min_{\lambda\in\Delta^{\star}}|\lambda|^{2}=\mu^{2}$ if
$L\geqslant(\sqrt{2}+1)\mu$, and
$\min_{\lambda\in\Delta^{\star}}|\lambda|^{2}=(L-\mu)^{2}/2$ otherwise.

With the above analysis, we can translate the convergence rates of GD and EG
in (25) and (26), respectively, for solving (1) with the Jacobian spectrum
structure in (2). We summarize this in the next remark.

######  Remark 2.

With the conditions in Theorem 4, for games with cross-shaped Jacobian
spectrum as in (2) with $c=\frac{L-\mu}{2}$, the above analysis leads to the
following local convergence rates of GD and EG:

$\displaystyle\rho(\nabla
F_{{\color[rgb]{0,0.75,0.16}h}}^{\text{GD}}(w^{\star}))^{2}\leqslant
1-\tfrac{\mu}{L},~{}~{}\text{and}~{}~{}$
$\displaystyle\rho_{\text{EG}}^{2}\leqslant\begin{cases}1-\tfrac{1}{4}\left(\tfrac{\mu}{L}+\tfrac{\mu^{2}}{16L^{2}}\right)&\text{if}~{}L\geqslant(\sqrt{2}+1)\mu,\\\
1-\tfrac{1}{4}\left(\tfrac{\mu}{L}+\tfrac{(L-\mu)^{2}}{16L^{2}}\right)&\text{otherwise}.\end{cases}$

We see that both GD and EG have non-accelerated convergence rates, unlike the
rate of EGM we obtained in (21), for the same problem class.

We now compare the convergence rate of EGM with that of GDM, which iterates as
in (7). In [3], it was shown that GD is optimal for the class of games where
the eigenvalues of the Jacobian is included in a _disc_ in the complex plane,
indicating acceleration is not possible for this problem class. However, it is
well-known that GDM achieves the optimal convergence rate when the the the
eigenvalues of the Jacobian lie on the (strictly positive) real line segment,
which amounts to (strongly-convex) minimization. Hence, [3] studies the
intermediate case, when the eigenvalues of the Jacobian are included in an
ellipse, which can be thought of as the real segment $[\mu,L]$ perturbed with
$\epsilon$ in an elliptic way. That is, they consider the following spectral
shape:666A visual illustration of this ellipse can be found in [3, Figure 2].

$K_{\epsilon}=\left\\{z\in\mathbb{C}:\left(\tfrac{\mathfrak{R}z-\tfrac{\mu+L}{2}}{\tfrac{L-\mu}{2}}\right)^{2}+\left(\tfrac{\mathfrak{I}z}{\epsilon}\right)^{2}\right\\}.$

Similarly to GD and EG above, in [3], GDM is interpreted as a fixed point
iteration,777As the GDM recursion for $w_{t+1}$ depends on $w_{t}$ and
$w_{t-1}$, [3] uses an augmented operator; see Lemma 2 in that work for
details. as in

$\displaystyle w_{t+1}$
$\displaystyle=w_{t}-{\color[rgb]{0,0.75,0.16}h}v(w_{t})+{\color[rgb]{.75,0,.25}m}(w_{t}-w_{t-1})$
$\displaystyle=F^{\text{GDM}}(w_{t},w_{t-1}).$ (27)

To study the convergence rate of GDM, we use the following theorem from [3]:

######  Theorem 5 ([3]).

Define $\epsilon(\mu,L)$ as $\epsilon(\mu,L)/L=(\mu/L)^{\theta}=\tau^{\theta}$
with $\theta>0$ and $a\land b=\min(a,b).$ If $\text{Sp}(\nabla
F^{\text{GDM}}(w^{\star},w^{\star}))\subset K_{\epsilon}$, and when
$\tau\rightarrow 0,$ it satisfies:

$\rho_{\text{GDM}}\leqslant\begin{cases}1-2\sqrt{\tau}+O\left(\tau^{\theta\land
1}\right),&\text{if}~{}~{}\theta>\frac{1}{2}\\\
1-2(\sqrt{2}-1)\sqrt{\tau}+O\left(\tau\right),&\text{if}~{}~{}\theta=\frac{1}{2}\\\
1-\tau^{1-\theta}+O\left(\tau^{1\land(2-3\theta)}\right),&\text{if}~{}~{}\theta<\frac{1}{2},\end{cases}$

where we denoted $\rho_{\text{GDM}}:=\rho(\nabla
F^{\text{GDM}}(w^{\star},w^{\star}))$ for notational brevity. Note that the
hyperparameters ${\color[rgb]{0,0.75,0.16}h}$ and ${\color[rgb]{.75,0,.25}m}$
are functions of $\mu,L,$ and $\epsilon$ only.

From Theorem 5, we see that GDM achieves the accelerated rate, i.e.,
$1-O(\sqrt{\tau}),$ until $\theta=\frac{1}{2}$. In other words, the biggest
elliptic perturbation $\epsilon$ where GDM permits the accelerated rate is
$\epsilon=\sqrt{\mu L}.$888Observe that
$(\mu/L)^{1/2}=\epsilon(\mu,L)/L\implies\epsilon(\mu,L)=\sqrt{\mu L}$. We
interpret Theorem 5 for games with cross-shaped Jacobian spectrum in the
following remark:

######  Remark 3.

With the conditions in Theorem 5, for games with cross-shaped Jacobian
spectrum in (2) with $c=\frac{L-\mu}{2}$, GDM cannot achieve an accelerated
rate when

$\displaystyle\tfrac{L-\mu}{2}=c>\epsilon=\sqrt{\mu L}.$

Using $L>\mu,$ the above inequality implies $\frac{L}{\mu}>\sqrt{5}.$ That is,
when the condition number exceeds $\sqrt{5}\approx 2.236$, GDM cannot be
accelerated; on the contrary, as we showed in (20) and (21), EGM converges at
the accelerated rate in the ill-conditioned regime.

##  4 EXPERIMENTS

In this section, we perform synthetic experiments on optimizing a game where
the Jacobian has a cross-shaped spectrum in (2). Specifically, we focus on
two-player quadratic games, where player 1 controls $x\in\mathbb{R}^{d_{1}}$
with loss function (3), and player 2 controls $y\in\mathbb{R}^{d_{2}}$ with
loss function (4).

In our setting, the corresponding vector field in (5) satisfies
$M_{12}=-M_{21}^{\top}$, but $S_{1}$ and $S_{2}$ can be nonzero symmetric
matrices. Further, the Jacobian $\nabla v=A$ has the cross-shape eigenvalue
structure in (2), with $c=\frac{L-\mu}{2}$ to reflect the analyses we provided
in Section 3.2.

We utilize Proposition 1 to design our experiments. We briefly describe how we
generate $A$ with the desired spectrum. Notice that
$\text{Sp}\left(\begin{bmatrix}a&-b\\\ b&a\end{bmatrix}\right)=\\{a\pm bi\\}.$
Hence, we generate several $2\times 2$ matrices of the above form with the
desired eigenvalues that align with (2), and put them in a block diagonal
matrix. We then perform a change of basis to keep the same eigenvalues, but
make the experiments more realistic.

For the problem constants, we use $\mu=1$, and $L=200.$ The optimum
$[x^{\star}~{}y^{\star}]^{\top}=w^{\star}\in\mathbb{R}^{200}$ is generated
using the standard normal distribution. For simplicity, we assume
$b=[b_{1}~{}b_{2}]^{\top}=[0~{}0]^{\top}.$

We compare GD in (23), GDM in (7), EG in (24), and EGM in (6). All algorithms
are initialized with $0$. We plot the experimental results in Figure 3.

For EGM (optimal), we set the hyperparameters using Corollary 1. For GD
(theory) and EG (theory), we set the hyperparameters using Theorem 4. For GDM
(grid search), as Theorem 5 does not give a specific form for hyperparater
setup, we perform grid search of ${\color[rgb]{0,0.75,0.16}h}^{\text{GDM}}$
and ${\color[rgb]{.75,0,.25}m}^{\text{GDM}}$, and choose the best performing
ones. Specifically, we consider
$0.005\leqslant{\color[rgb]{0,0.75,0.16}h}^{\text{GDM}}\leqslant 0.015$ with
$10^{-3}$ increment, and
$0.01\leqslant{\color[rgb]{.75,0,.25}m}^{\text{GDM}}\leqslant 0.99$ with
$10^{-2}$ increment. Further, as Theorem 4 might be conservative, we perform
grid searches for GD and EG as well. For GD (grid search), we use the same
setup as ${\color[rgb]{0,0.75,0.16}h}^{\text{GDM}}$. For EG (grid search), we
use $0.001\leqslant{\color[rgb]{0,0.75,0.16}h}^{\text{EG}}\leqslant 0.05$ with
$10^{-4}$ increment.

There are several remarks to make. First, EGM (optimal) with the
hyperparameters obtained in Corollary 1 indeed exhibit orders of magnitude
faster rate of convergence, even compared to other methods with grid-search
hyperparameter tuning, corroborating our theoretical findings in Section 3.
Second, while EG (theory) is slower than GD (theory), confirming Remark 2, GD
(grid search) can be tuned to converge faster via grid search. Third, even
though the best performance of GDM (grid search) is obtained through grid
search, one can see the GD (grid search) obtains slightly faster convergence
rate than GDM (grid search), confirming Remark 3.

Figure 3: Comparison of the first-order methods considered in this work. All
algorithms are initialized with $0$, and run for 2000 iterations.

##  5 CONCLUSION

In this work, we focused on finding a stationary point of games. Via
polynomial-based analysis, we showed the extragradient with momentum converges
in three different modes of eigenvalue structure of the game Jacobian,
depending on the hyperparameters. Then, we focused on a more specific problem
class, where the eigenvalues of the game Jacobian are distributed in the
complex plane with a cross-shape. For this problem class, we obtained the
optimal hyperparameters for the momentum extragradient method, which enjoys an
accelerated linear asymptotic convergence rate. We compared the obtained rates
of the momentum extragradient method with other first order methods, and
showed that the considered methods do not achieve the accelerated convergence
rate.

## References

  * [1] Yossi Arjevani and Ohad Shamir  “On the iteration complexity of oblivious first-order optimization algorithms”  In _International Conference on Machine Learning_ , 2016, pp. 908–916  PMLR
  * [2] Waı̈ss Azizian, Ioannis Mitliagkas, Simon Lacoste-Julien and Gauthier Gidel  “A tight and unified analysis of gradient-based methods for a whole spectrum of differentiable games”  In _International Conference on Artificial Intelligence and Statistics_ , 2020, pp. 2863–2873  PMLR
  * [3] Waı̈ss Azizian et al.  “Accelerating smooth games by manipulating spectral shapes”  In _International Conference on Artificial Intelligence and Statistics_ , 2020, pp. 1705–1715  PMLR
  * [4] David Balduzzi et al.  “The mechanics of n-player differentiable games”  In _International Conference on Machine Learning_ , 2018, pp. 354–363  PMLR
  * [5] Hugo Berard et al.  “A Closer Look at the Optimization Landscapes of Generative Adversarial Networks”  In _International Conference on Learning Representations_ , 2020  URL: https://openreview.net/forum?id=HJeVnCEKwH
  * [6] Raphaël Berthier, Francis Bach and Pierre Gaillard  “Accelerated gossip in networks of given dimension using Jacobi polynomial iterations”  In _SIAM Journal on Mathematics of Data Science_ 2.1 SIAM, 2020, pp. 24–47
  * [7] Aleksandr Beznosikov et al.  “Decentralized local stochastic extra-gradient for variational inequalities”  In _arXiv preprint arXiv:2106.08315_ , 2021
  * [8] Theodore S Chihara  “An introduction to orthogonal polynomials”  Courier Corporation, 2011
  * [9] Constantinos Daskalakis, Paul W Goldberg and Christos H Papadimitriou  “The complexity of computing a Nash equilibrium”  In _SIAM Journal on Computing_ 39.1 SIAM, 2009, pp. 195–259
  * [10] Constantinos Daskalakis, Andrew Ilyas, Vasilis Syrgkanis and Haoyang Zeng  “Training gans with optimism”  In _arXiv preprint arXiv:1711.00141_ , 2017
  * [11] Francisco Facchinei and Jong-Shi Pang  “Finite-dimensional variational inequalities and complementarity problems”  Springer, 2003
  * [12] Gauthier Gidel et al.  “A variational inequality perspective on generative adversarial networks”  In _arXiv preprint arXiv:1802.10551_ , 2018
  * [13] Gauthier Gidel et al.  “Negative momentum for improved game dynamics”  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , 2019, pp. 1802–1811  PMLR
  * [14] Ian Goodfellow et al.  “Generative Adversarial Nets”  In _Advances in Neural Information Processing Systems_ 27 Curran Associates, Inc., 2014  URL: https://proceedings.neurips.cc/paper/2014/file/5ca3e9b122f61f8f06494c97b1afccf3-Paper.pdf
  * [15] Eduard Gorbunov, Nicolas Loizou and Gauthier Gidel  “Extragradient method: O (1/K) last-iterate convergence for monotone variational inequalities and connections with cocoercivity”  In _International Conference on Artificial Intelligence and Statistics_ , 2022, pp. 366–402  PMLR
  * [16] Baptiste Goujaud and Fabian Pedregosa  “Cyclical Step-sizes”, 2022  URL: http://fa.bianp.net/blog/2022/cyclical/
  * [17] Baptiste Goujaud et al.  “Super-acceleration with cyclical step-sizes”  In _International Conference on Artificial Intelligence and Statistics_ , 2022, pp. 3028–3065  PMLR
  * [18] Magnus R Hestenes and Eduard Stiefel  “Methods of conjugate gradients for solving”  In _Journal of research of the National Bureau of Standards_ 49.6, 1952, pp. 409
  * [19] Galina M Korpelevich  “The extragradient method for finding saddle points and other problems”  In _Matecon_ 12, 1976, pp. 747–756
  * [20] Jonathan Lacotte and Mert Pilanci  “Optimal randomized first-order methods for least-squares problems”  In _International Conference on Machine Learning_ , 2020, pp. 5587–5597  PMLR
  * [21] Alistair Letcher et al.  “Differentiable game mechanics”  In _The Journal of Machine Learning Research_ 20.1 JMLR. org, 2019, pp. 3032–3071
  * [22] Chris Junchi Li et al.  “On the convergence of stochastic extragradient for bilinear games with restarted iteration averaging”  In _arXiv preprint arXiv:2107.00464_ , 2021
  * [23] Tengyuan Liang and James Stokes  “Interaction matters: A note on non-asymptotic local convergence of generative adversarial networks”  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , 2019, pp. 907–915  PMLR
  * [24] Mingrui Liu et al.  “A decentralized parallel algorithm for training generative adversarial nets”  In _Advances in Neural Information Processing Systems_ 33, 2020, pp. 11056–11070
  * [25] Lars Mescheder, Andreas Geiger and Sebastian Nowozin  “Which training methods for GANs do actually converge?”  In _International conference on machine learning_ , 2018, pp. 3481–3490  PMLR
  * [26] Lars Mescheder, Sebastian Nowozin and Andreas Geiger  “The numerics of gans”  In _Advances in neural information processing systems_ 30, 2017
  * [27] Konstantin Mishchenko et al.  “Revisiting stochastic extragradient”  In _International Conference on Artificial Intelligence and Statistics_ , 2020, pp. 4573–4582  PMLR
  * [28] Aryan Mokhtari, Asuman Ozdaglar and Sarath Pattathil  “A unified analysis of extra-gradient and optimistic gradient methods for saddle point problems: Proximal point approach”  In _International Conference on Artificial Intelligence and Statistics_ , 2020, pp. 1497–1507  PMLR
  * [29] Arkadi 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”  In _SIAM Journal on Optimization_ 15.1 SIAM, 2004, pp. 229–251
  * [30] Arkadij Semenovič Nemirovskij and David Borisovich Yudin  “Problem complexity and method efficiency in optimization”  Wiley-Interscience, 1983
  * [31] Yurii Nesterov  “Introductory lectures on convex optimization: A basic course”  Springer Science & Business Media, 2003
  * [32] Yurii Evgen’evich Nesterov  “A method of solving a convex programming problem with convergence rate $O(k^{2})$”  In _Doklady Akademii Nauk_ 269.3, 1983, pp. 543–547  Russian Academy of Sciences
  * [33] Roy A Nicolaides  “Existence, uniqueness and approximation for generalized saddle point problems”  In _SIAM Journal on Numerical Analysis_ 19.2 SIAM, 1982, pp. 349–357
  * [34] Balamurugan Palaniappan and Francis Bach  “Stochastic variance reduction methods for saddle-point problems”  In _Advances in Neural Information Processing Systems_ 29, 2016
  * [35] Fabian Pedregosa  “Momentum: when Chebyshev meets Chebyshev”, 2020  URL: http://fa.bianp.net/blog/2020/momentum/
  * [36] Fabian Pedregosa and Damien Scieur  “Acceleration through spectral density estimation” ISSN: 2640-3498  In _Proceedings of the 37th International Conference on Machine Learning_ PMLR, 2020, pp. 7553–7562  URL: https://proceedings.mlr.press/v119/pedregosa20a.html
  * [37] David Pfau and Oriol Vinyals  “Connecting generative adversarial networks and actor-critic methods”  In _arXiv preprint arXiv:1610.01945_ , 2016
  * [38] Boris T Polyak  “Introduction to optimization. optimization software”  In _Inc., Publications Division, New York_ 1, 1987, pp. 32
  * [39] Ernest K Ryu, Kun Yuan and Wotao Yin  “Ode analysis of stochastic gradient methods with optimism and anchoring for minimax problems”  In _arXiv preprint arXiv:1905.10899_ , 2019
  * [40] Damien Scieur and Fabian Pedregosa  “Universal average-case optimality of polyak momentum”  In _International conference on machine learning_ , 2020, pp. 8565–8572  PMLR
  * [41] Yoav Shoham and Kevin Leyton-Brown  “Multiagent systems: Algorithmic, game-theoretic, and logical foundations”  Cambridge University Press, 2008
  * [42] Mikhail V Solodov and Benar F Svaiter  “A hybrid approximate extragradient–proximal point algorithm using the enlargement of a maximal monotone operator”  In _Set-Valued Analysis_ 7.4 Springer, 1999, pp. 323–345
  * [43] Sainbayar Sukhbaatar et al.  “Intrinsic motivation and automatic curricula via asymmetric self-play”  In _arXiv preprint arXiv:1703.05407_ , 2017
  * [44] Paul Tseng  “On linear convergence of iterative methods for the variational inequality problem”  In _Journal of Computational and Applied Mathematics_ 60.1-2 Elsevier, 1995, pp. 237–252

##  Appendix A Missing Proofs

###  A.1 Proof of Proposition 1

###### Proof.

Define $D_{3}=\text{diag}(a,\ldots,a)$ of dimensions $d\times d$. Let us prove
that if there exists $U,V$ orthonormal matrices and $D_{1},D_{2}$ matrices
with non-zeros coefficients only on the diagonal such that (with a slight
abuse of notation)

$\displaystyle
S_{1}=U\text{diag}(D_{3},D_{1})U^{\top},S_{2}=VD_{3}V^{\top},\quad\text{and}\quad
B=UD_{2}V^{\top},$

then the spectrum of $A$ is crossed shaped. In that case, we have

$\displaystyle A$
$\displaystyle=\begin{bmatrix}U[D_{3};D_{1}]U^{\top}&UD_{2}V^{\top},\\\
-VD_{2}^{\top}U^{\top}&VD_{3}V^{\top}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}U&0\\\
0&V\end{bmatrix}\begin{bmatrix}[D_{3};D_{1}]&D_{2}\\\
-D_{2}^{\top}&D_{3}\end{bmatrix}\begin{bmatrix}U&0\\\
0&V\end{bmatrix}^{\top}.$

Now by considering the basis
$W=((U_{1},0),(0,V_{1}),\ldots,(U_{d_{v}},0),(0,V_{d_{v}}),(U_{d_{v}+1},0),\ldots,(U_{d},0))$
we have that $A$ can be block diagonalized in that basis as

$A=W\text{diag}\left(\begin{bmatrix}a&[D_{2}]_{11}\\\
-[D_{2}]_{11}&a\end{bmatrix},\ldots,\begin{bmatrix}a&[D_{2}]_{d_{v},d_{v}}\\\
-[D_{2}]_{d_{v}d_{v}}&a\end{bmatrix},[D_{1}]_{1},\dots,[D_{1}]_{d_{u}-d_{v}d_{u}-d_{v}}\right)W^{\top}.$
(28)

Now, notice that

$\displaystyle\text{Sp}\left(\begin{bmatrix}a&-b\\\
b&a\end{bmatrix}\right)=\\{a\pm bi\\},$ (29)

since the associated characteristic polynomial of the above matrix is:

$\displaystyle(a-\lambda)^{2}+b^{2}$ $\displaystyle=0\implies a-\lambda=\pm
bi\implies\lambda=a\pm bi.$

Hence, using (29) in the formulation of $A$ in (28), we have that the spectrum
of $A$ is cross-shaped.

∎

###  A.2 Proof of Lemma 1

Proof of Lemma 1 can be found for example in [3, Section B (appendix)].

###  A.3 Proof of Theorem 1

#### Derivation of the first part

###### Proof.

We want to find the residual polynomial $\tilde{P}_{t}(A)$ of the
extragradient with momentum (EGM) in (6). That is, we want to find

$w_{t}-w^{\star}=\tilde{P}_{t}(A)(w_{0}-w^{\star}),$ (30)

where $\\{w_{t}\\}_{t=0}$ is the iterates generated by EGM, which is possible
by Lemma 1, as EGM is a first-order method [1, 3]. We now prove this is by
induction. To do so, we will use the following properties. First, note that as
we are looking for a stationary point, it holds that $v(w^{\star})=0.$
Further, as $v$ is linear by the assumption of Lemma 1, it holds that
$v(w)=A(w-w^{\star}).$

Base case. For $t=0$, $\tilde{P}_{0}(A)$ is a degree-zero polynomial, and
hence equals $\mathbb{I}$, which denotes the identity matrix. Thus,
$w_{0}-w^{\star}=\mathbb{I}(w_{0}-w^{\star})$ holds true.

For completeness, we also prove when $t=1.$ In that case, observe that EGM in
Algorithm 1 proceeds as $w_{1}=w_{0}-\frac{h}{1+m}v(w_{0}-\gamma v(w_{0}))$.
Subtracting $w^{\star}$ on both sides, we have:

$\displaystyle w_{1}-w^{\star}$
$\displaystyle=w_{0}-w^{\star}-\frac{h}{1+m}v(w_{0}-\gamma v(w_{0}))$
$\displaystyle=w_{0}-w^{\star}-\frac{h}{1+m}v(w_{0}-\gamma
A(w_{0}-w^{\star}))$
$\displaystyle=w_{0}-w^{\star}-\frac{h}{1+m}A(w_{0}-\gamma
A(w_{0}-w^{\star})-w^{\star})$
$\displaystyle=w_{0}-w^{\star}-\frac{h}{1+m}A(w_{0}-w^{\star})+\frac{h\gamma}{1+m}A^{2}(w_{0}-w^{\star})$
$\displaystyle=\left(\mathbb{I}-\frac{h}{1+m}A+\frac{h\gamma}{1+m}A^{2}\right)(w_{0}-w^{\star})$
$\displaystyle=\left(\mathbb{I}-\frac{h}{1+m}A(\mathbb{I}-\gamma
A)\right)(w_{0}-w^{\star})$ $\displaystyle=\tilde{P}_{1}(A)(w_{0}-w^{\star}).$

Induction step. As the induction hypothesis, assume $\tilde{P}_{t}$ satisfies
(30). We want to prove this holds for $t+1.$ We have:

$\displaystyle w_{t+1}$ $\displaystyle=w_{t}-hv(w_{t}-\gamma
v(w_{t}))+m(w_{t}-w_{t-1})$ $\displaystyle=w_{t}-hv(w_{t}-\gamma
A(w_{t}-w^{\star}))+m(w_{t}-w_{t-1})$ $\displaystyle=w_{t}-hA(w_{t}-\gamma
A(w_{t}-w^{\star})-w^{\star})+m(w_{t}-w_{t-1})$
$\displaystyle=w_{t}-hA(w_{t}-w^{\star})+h\gamma
A^{2}(w_{t}-w^{\star})+m(w_{t}-w_{t-1})$
$\displaystyle=w_{t}-hA(\mathbb{I}-\gamma
A)(w_{t}-w^{\star})+m(w_{t}-w_{t-1}).$

Subtracting $w^{\star}$ on both sides, we have:

$\displaystyle w_{t+1}-w^{\star}$
$\displaystyle=w_{t}-w^{\star}-hA(\mathbb{I}-\gamma
A)(w_{t}-w^{\star})+m(w_{t}-w_{t-1})$
$\displaystyle=\left(\mathbb{I}-hA(\mathbb{I}-\gamma
A)\right)(w_{t}-w^{\star})+m(w_{t}-w^{\star}-(w_{t-1}-w^{\star}))$
$\displaystyle\overset{\eqref{eq:ind-hypothesis-
egm}}{=}\left(\mathbb{I}-hA(\mathbb{I}-\gamma
A)\right)\tilde{P}_{t}(A)(w_{0}-w^{\star})+m(\tilde{P}_{t}(A)(w_{0}-w^{\star})-\tilde{P}_{t-1}(A)(w_{0}-w^{\star}))$
$\displaystyle=\left(\mathbb{I}+m\mathbb{I}-hA(\mathbb{I}-\gamma
A)\right)\tilde{P}_{t}(A)(w_{0}-w^{\star})-m\tilde{P}_{t-1}(A)(w_{0}-w^{\star})$
$\displaystyle=\tilde{P}_{t+1}(A)(w_{0}-w^{\star}),$

where in the third equality, we used the induction hypothesis in (30).

∎

#### Derivation of the second part in (10)

###### Proof.

We show $P_{t}=\tilde{P}_{t}$ for all $t$ via induction.

Base case. For $t=0,$ by the definition of Chebyshev polynomials of the first
and the second kinds, we have $T_{0}(\lambda)=U_{0}(\lambda)=1$. Thus,

$\displaystyle P_{0}(\lambda)$
$\displaystyle=m^{0}\left({\frac{2m}{1+m}}T_{0}(\sigma(\lambda))+{\frac{1-m}{1+m}}U_{0}(\sigma(\lambda))\right)$
$\displaystyle=\frac{2m}{1+m}+\frac{1-m}{1+m}=1=\tilde{P}_{0}(\lambda).$

Again, for completeness, we prove when $t=1$ as well. In that case, by the
definition of Chebyshev polynomials of the first and the second kinds, we have
$T_{1}(\lambda)=\lambda,$ and $U_{1}(\lambda)=2\lambda.$ Therefore,

$\displaystyle P_{1}(\lambda)$
$\displaystyle=m^{t/2}\left({\frac{2m}{1+m}}T_{1}(\sigma(\lambda))+{\frac{1-m}{1+m}}U_{1}(\sigma(\lambda))\right)$
$\displaystyle=m^{t/2}\left({\frac{2m}{1+m}}\,\sigma(\lambda)+{\frac{1-m}{1+m}}\cdot
2\cdot\sigma(\lambda)\right)$
$\displaystyle=m^{t/2}\left(\frac{2\sigma(\lambda)}{1+m}\right)$
$\displaystyle=1-\frac{h\lambda(1-\gamma\lambda)}{1+m}=\tilde{P}_{1}(\lambda).$

Induction step. As the induction hypothesis, assume that $P_{t}=\tilde{P}_{t}$
for $t$. We want to show this holds for $t+1.$

$\displaystyle P_{t+1}$
$\displaystyle=m^{(t+1)/2}\bigg{[}\frac{2m}{1+m}T_{t+1}(\sigma(\lambda))+\frac{1-m}{1+m}U_{t+1}(\sigma(\lambda))\bigg{]}$
$\displaystyle=m^{(t+1)/2}\bigg{[}\frac{2m}{1+m}\bigg{(}2\sigma(\lambda)T_{t}(\sigma(\lambda))-T_{t-1}(\sigma(\lambda))\bigg{)}$
$\displaystyle\hskip
48.36958pt+\frac{1-m}{1+m}\bigg{(}2\sigma(\lambda)U_{t}(\sigma(\lambda)-U_{t-1}(\sigma(\lambda)))\bigg{)}\bigg{]}$
$\displaystyle=2\sigma(\lambda)\cdot
m^{1/2}\cdot\underbrace{m^{t/2}\bigg{(}\frac{2m}{1+m}T_{t}(\sigma(\lambda))+\frac{1-m}{1+m}U_{t}(\sigma(\lambda))\bigg{)}}_{P_{t}(\lambda)}$
$\displaystyle\hskip
48.36958pt-m\cdot\underbrace{m^{(t-1)/2}\bigg{(}\frac{2m}{1+m}T_{t-1}(\sigma(\lambda))+\frac{1-m}{1+m}U_{t-1}(\sigma(\lambda))\bigg{)}}_{P_{t-1}(\lambda)}$
$\displaystyle=2\sigma(\lambda)\cdot\sqrt{m}\cdot\tilde{P}_{t}(\lambda)-m\cdot\tilde{P}_{t-1}(\lambda)$
$\displaystyle=(1+m-h\lambda(1-\gamma\lambda))\tilde{P}_{t}(\lambda)-m\tilde{P}_{t-1}(\lambda),$

where in the second to last equality we use the induction hypothesis. ∎

###  A.4 Proof of Lemma 2

Proof of Lemma 2 can be found in [16].

###  A.5 Proof of Theorem 2

###### Proof.

We analyze each case separately.

Case 1: There are two square roots:
$\sqrt{\frac{1}{4\gamma^{2}}-\frac{(1-\sqrt{m})^{2}}{h\gamma}}$ and
$\sqrt{\frac{1}{4\gamma^{2}}-\frac{(1+\sqrt{m})^{2}}{h\gamma}}$. The second
one is real if:

$\frac{1}{4\gamma^{2}}\geqslant\frac{(1+\sqrt{m})^{2}}{h\gamma}\implies\frac{h\gamma}{4\gamma^{2}}=\frac{h}{4\gamma}\geqslant(1+\sqrt{m})^{2},$

which implies the first is real, as
$(1+\sqrt{m})^{2}\geqslant(1-\sqrt{m})^{2}$.

Case 3: There are two square roots:
$\sqrt{\frac{1}{4\gamma^{2}}-\frac{(1-\sqrt{m})^{2}}{h\gamma}}$ and
$\sqrt{\frac{1}{4\gamma^{2}}-\frac{(1+\sqrt{m})^{2}}{h\gamma}}$. The first one
is complex if:

$\frac{1}{4\gamma^{2}}<\frac{(1-\sqrt{m})^{2}}{h\gamma}\implies\frac{h\gamma}{4\gamma^{2}}=\frac{h}{4\gamma}<(1-\sqrt{m})^{2},$

which implies the second is complex, as
$(1+\sqrt{m})^{2}\geqslant(1-\sqrt{m})^{2}$.

Case 2: This case follows automatically from the above two cases.

∎

###  A.6 Proof of Theorem 3

###### Proof.

We recall the conditions required for Theorem 3 below:

$\displaystyle\frac{1}{2\gamma}-\sqrt{\frac{1}{4\gamma^{2}}-\frac{(1-\sqrt{m})^{2}}{h\gamma}}$
$\displaystyle=\mu,$ (31)
$\displaystyle\frac{1}{2\gamma}+\sqrt{\frac{1}{4\gamma^{2}}-\frac{(1-\sqrt{m})^{2}}{h\gamma}}$
$\displaystyle=L,\quad\text{and}$ (32)
$\displaystyle\sqrt{\frac{(1+\sqrt{m})^{2}}{h\gamma}-\frac{1}{4\gamma^{2}}}$
$\displaystyle=c.$ (33)

First, by adding (31) and (32), we get:

$\displaystyle\frac{1}{\gamma}=\mu+L\implies\gamma=\frac{1}{\mu+L}.$ (34)

Plugging (34) back into (31), we have:

$\displaystyle\frac{1}{2\gamma}-\sqrt{\frac{1}{4\gamma^{2}}-\frac{(1-\sqrt{m})^{2}}{h\gamma}}$
$\displaystyle=\mu$ $\displaystyle\frac{\mu+L}{2}-\mu$
$\displaystyle=\sqrt{\left(\frac{\mu+L}{2}\right)^{2}-\frac{(1-\sqrt{m})^{2}(\mu+L)}{h}}$
$\displaystyle\left(\frac{L-\mu}{2}\right)^{2}$
$\displaystyle=\left(\frac{\mu+L}{2}\right)^{2}-\frac{(1-\sqrt{m})^{2}(\mu+L)}{h}$
$\displaystyle\frac{(1-\sqrt{m})^{2}(\mu+L)}{h}$
$\displaystyle=\left(\frac{\mu+L}{2}\right)^{2}-\left(\frac{L-\mu}{2}\right)^{2}=\mu
L$ $\displaystyle h$ $\displaystyle=\frac{(1-\sqrt{m})^{2}(\mu+L)}{\mu L}.$
(35)

Plugging (34) and (35) into (33), we have:

$\displaystyle\sqrt{\frac{(1+\sqrt{m})^{2}}{h\gamma}-\frac{1}{4\gamma^{2}}}$
$\displaystyle=c$ $\displaystyle\sqrt{\frac{(1+\sqrt{m})^{2}\cdot\mu
L}{(1-\sqrt{m})^{2}}-\left(\frac{\mu+L}{2}\right)^{2}}$ $\displaystyle=c$
$\displaystyle\frac{(1+\sqrt{m})^{2}\cdot\mu L}{(1-\sqrt{m})^{2}}$
$\displaystyle=c^{2}+\left(\frac{\mu+L}{2}\right)^{2}=\frac{4c^{2}+(\mu+L)^{2}}{4}$
$\displaystyle\frac{(1+\sqrt{m})^{2}}{(1-\sqrt{m})^{2}}$
$\displaystyle=\frac{4c^{2}+(\mu+L)^{2}}{4\mu L}$
$\displaystyle(1+\sqrt{m})\sqrt{4\mu L}$
$\displaystyle=(1-\sqrt{m})\sqrt{4c^{2}+(\mu+L)^{2}}$
$\displaystyle\sqrt{m}(\sqrt{4c^{2}+(\mu+L)^{2}}+\sqrt{4\mu L})$
$\displaystyle=\sqrt{4c^{2}+(\mu+L)^{2}}-\sqrt{4\mu L}$
$\displaystyle\sqrt{m}$
$\displaystyle=\frac{\sqrt{4c^{2}+(\mu+L)^{2}}-\sqrt{4\mu
L}}{\sqrt{4c^{2}+(\mu+L)^{2}}+\sqrt{4\mu L}}.$ (36)

Finally, to simplify (35) further, from (36), we have:

$\displaystyle 1-\sqrt{m}$ $\displaystyle=\frac{4\sqrt{\mu
L}}{\sqrt{4c^{2}+(\mu+L)^{2}}+\sqrt{4\mu L}}.$

Hence, from (35),

$\displaystyle h$ $\displaystyle=\frac{(\mu+L)(1-\sqrt{m})^{2}}{\mu
L}=\frac{\frac{16\mu L(\mu+L)}{(\sqrt{4c^{2}+(\mu+L)^{2}}+\sqrt{4\mu
L})^{2}}}{\mu L}=\frac{16(\mu+L)}{(\sqrt{4c^{2}+(\mu+L)^{2}}+\sqrt{4\mu
L})^{2}}.$ (37)

∎

###  A.7 Proof of Corollary 1

###### Proof.

This result can be obtained either by plugging in $c=\frac{L-\mu}{2}$ to (36)
and (37), or by using the same technique as the proof of Theorem 3. For
completeness, we prove by plugging in $c=\frac{L-\mu}{2}$ to (36) and (37)
here.

First, for $\sqrt{m}$, from (36), we have:

$\displaystyle\sqrt{m}$
$\displaystyle=\frac{\sqrt{4c^{2}+(\mu+L)^{2}}-\sqrt{4\mu
L}}{\sqrt{4c^{2}+(\mu+L)^{2}}+\sqrt{4\mu
L}}=\frac{\sqrt{4\left(\frac{L-\mu}{2}\right)^{2}+(\mu+L)^{2}}-\sqrt{4\mu
L}}{\sqrt{4\left(\frac{L-\mu}{2}\right)^{2}+(\mu+L)^{2}}+\sqrt{4\mu L}}$
$\displaystyle=\frac{\sqrt{2\mu^{2}+2L^{2}}-\sqrt{4\mu
L}}{\sqrt{2\mu^{2}+2L^{2}}+\sqrt{4\mu
L}}=\frac{\sqrt{\mu^{2}+L^{2}}-\sqrt{2\mu L}}{\sqrt{\mu^{2}+L^{2}}+\sqrt{2\mu
L}}.$

Now, for $h$, focusing on the denominator of (37):

$\displaystyle\left(\sqrt{4c^{2}+(\mu+L)^{2}}+\sqrt{4\mu L}\right)^{2}$
$\displaystyle=\left(\sqrt{4\left(\frac{L-\mu}{2}\right)^{2}+(\mu+L)^{2}}+\sqrt{4\mu
L}\right)^{2}$ $\displaystyle=\left(\sqrt{2\mu^{2}+2L^{2}}+\sqrt{2\mu
L}\right)^{2}$ $\displaystyle=2\mu^{2}+2L^{2}+4\mu L+4\sqrt{\mu
L(2\mu^{2}+2L^{2})}.$

Plugging this back into (37), we have:

$\displaystyle\frac{16(\mu+L)}{(\sqrt{4c^{2}+(\mu+L)^{2}}+\sqrt{4\mu L})^{2}}$
$\displaystyle=\frac{16(\mu+L)}{2\mu^{2}+2L^{2}+4\mu L+4\sqrt{\mu
L(2\mu^{2}+2L^{2})}}$ $\displaystyle=\frac{8(\mu+L)}{\mu^{2}+L^{2}+2\mu
L+2\sqrt{2\mu L(\mu^{2}+L^{2})}}$
$\displaystyle=\frac{8(\mu+L)}{(\sqrt{\mu^{2}+L^{2}}+\sqrt{2\mu L})^{2}}.$

∎

##  Appendix B Details on Remarks

###  B.1 Details on Remark 2

To compute the convergence rates of GD and EG from Theorem 4 applied to the
cross-shaped spectrum of game Jacobian in (2), we need to compute
$\min_{\lambda\in\Delta^{\star}}\mathfrak{R}(1/\lambda)$ and
$\min_{\lambda\in\Delta^{\star}}\mathfrak{R}(\lambda)$ for GD. For a complex
number $z=a+bi\in\mathbb{C},$ we can compute $\mathfrak{R}(1/z)$ as:

$\displaystyle\frac{1}{z}$
$\displaystyle=\frac{1}{a+bi}=\frac{a-bi}{a^{2}+b^{2}}=\frac{a}{a^{2}+b^{2}}-\frac{b}{a^{2}+b^{2}}i\implies\mathfrak{R}\left(\frac{1}{z}\right)=\frac{a}{a^{2}+b^{2}}.$

The four extreme points of the cross-shaped spectrum model in (2) with
$c=\frac{L-\mu}{2}$ are:

$\displaystyle\mu=\mu+0i,\quad
L=L+0i,\quad\text{and}\quad\frac{L-\mu}{2}\pm\frac{L-\mu}{2}i.$

Hence, $\mathfrak{R}(1/z)$ for each of the above points is:

$\displaystyle\mathfrak{R}\left(\frac{1}{\mu}\right)$
$\displaystyle=\frac{\mu}{\mu^{2}}=\frac{1}{\mu},$
$\displaystyle\mathfrak{R}\left(\frac{1}{L}\right)$
$\displaystyle=\frac{L}{L^{2}}=\frac{1}{L},\quad\text{and}$
$\displaystyle\mathfrak{R}\left(\frac{1}{\frac{L-\mu}{2}\pm\frac{L-\mu}{2}i}\right)$
$\displaystyle=\frac{\frac{L-\mu}{2}}{\left(\frac{L-\mu}{2}\right)^{2}+\left(\frac{L-\mu}{2}\right)^{2}}$
$\displaystyle=\frac{\frac{L-\mu}{2}}{\frac{2L^{2}+2\mu^{2}-4L\mu}{4}}=\frac{\frac{L-\mu}{2}}{\frac{(L-\mu)^{2}}{2}}=\frac{1}{L-\mu}.$

Therefore,
$\min_{\lambda\in\Delta^{\star}}\mathfrak{R}\left(\frac{1}{\lambda}\right)=\frac{1}{L}.$

For $\min_{\lambda\in\Delta^{\star}}\mathfrak{R}(\lambda)$, we have

$\displaystyle\mathfrak{R}(\mu)=\mu,\quad\mathfrak{R}(L)=L,\quad\text{and}\quad\mathfrak{R}\left(\frac{L-\mu}{2}\pm\frac{L-\mu}{2}i\right)=\frac{L-\mu}{2}.$

Therefore,

$\displaystyle\min_{\lambda\in\Delta^{\star}}\mathfrak{R}(\lambda)=\begin{cases}\mu\quad&\text{if}\quad
3\mu<L\\\ \frac{L-\mu}{2}&\text{otherwise}.\end{cases}$

As we consider the ill-conditioned regime where $\frac{\mu}{L}\rightarrow 0,$
we have $\min_{\lambda\in\Delta^{\star}}\mathfrak{R}(\lambda)=\mu.$ Thus, we
get

$\displaystyle\rho(\nabla F_{h}^{\text{GD}}(w^{\star}))^{2}\leqslant
1-\min_{\lambda\in\Delta^{\star}}\mathfrak{R}\left(\frac{1}{\lambda}\right)\min_{\lambda\in\Delta^{\star}}\mathfrak{R}(\lambda)=1-\frac{\mu}{L}.$

Similarly for EG, we need to compute
$\min_{\lambda\in\Delta^{\star}}\mathfrak{R}(\lambda)$, which was computed
above; additionally, we need to compute
$\min_{\lambda\in\Delta^{\star}}|\lambda|$,
$\min_{\lambda\in\Delta^{\star}}|\lambda|^{2},$ and
$\max_{\lambda\in\Delta^{\star}}|\lambda|^{2}$. For $z=a+bi\in\mathbb{C},$
$|z|=\sqrt{a^{2}+b^{2}}.$ Hence, we have

$\displaystyle|\mu+0i|=\mu,\quad|L+0i|=L,~{}~{}\text{and}~{}~{}\left|\frac{L-\mu}{2}\pm\frac{L-\mu}{2}i\right|=\frac{\sqrt{2}(L-\mu)}{2}.$

Thus, for $\min_{\lambda\in\Delta^{\star}}|\lambda|$, we have

$\displaystyle\min_{\lambda\in\Delta^{\star}}|\lambda|=\begin{cases}\mu\quad&\text{if}\quad\mu<\frac{\sqrt{2}}{2+\sqrt{2}}L\\\
\frac{\sqrt{2}(L-\mu)}{2}&\text{otherwise}.\end{cases}$

Similarly to the above argument, in the ill-conditioned regime, we have
$\min_{\lambda\in\Delta^{\star}}|\lambda|=\mu.$

For $\max_{\lambda\in\Delta^{\star}}|\lambda|$, we have:

$\displaystyle(2-\sqrt{2})L$ $\displaystyle>L>\mu>-\sqrt{2}\mu$ $\displaystyle
2L-\sqrt{2}L$ $\displaystyle>-\sqrt{2}\mu$ $\displaystyle 2L$
$\displaystyle>\sqrt{2}(L-\mu)$ $\displaystyle L$
$\displaystyle>\frac{\sqrt{2}(L-\mu)}{2}$
$\displaystyle\therefore\max_{\lambda\in\Delta^{\star}}|\lambda|$
$\displaystyle=L.$

Finally, for $\min_{\lambda\in\Delta^{\star}}|\lambda|^{2},$ we have
$|\lambda|^{2}\in\\{\mu^{2},L^{2},\frac{(L-\mu)^{2}}{2}\\}.$ Thus, we have

$\displaystyle\min_{\lambda\in\Delta^{\star}}|\lambda|^{2}=\begin{cases}\mu^{2}\quad&\text{if}\quad
L\geqslant(\sqrt{2}+1)\mu\\\
\frac{(L-\mu)^{2}}{2}&\text{otherwise}.\end{cases}$

Combining all three, we get:

$\displaystyle\rho_{\text{EG}}^{2}:=\rho(\nabla
F_{h}^{\text{EG}}(w^{\star}))\leqslant
1-\tfrac{1}{4}\left(\tfrac{\min_{\lambda\in\Delta^{\star}}\mathfrak{R}(\lambda)}{\max_{\lambda\in\Delta^{\star}}|\lambda|}+\tfrac{\min_{\lambda\in\Delta^{\star}}|\lambda|^{2}}{\max_{\lambda\in\Delta^{\star}}|\lambda|^{2}}\right)=\begin{cases}1-\tfrac{1}{4}\left(\tfrac{\mu}{L}+\tfrac{\mu^{2}}{16L^{2}}\right)&\text{if}~{}L\geqslant(\sqrt{2}+1)\mu,\\\
1-\tfrac{1}{4}\left(\tfrac{\mu}{L}+\tfrac{(L-\mu)^{2}}{16L^{2}}\right)&\text{otherwise}.\end{cases}$

###  B.2 Details on Remark 3

Per Theorem 5, the largest $\epsilon$ that permits acceleration for GDM is
$\epsilon=\sqrt{\mu L}$. Therefore, in the special case of (2) we consider,
i.e., when $c=\frac{L-\mu}{2},$ GDM _cannot_ achieve acceleration if
$\frac{L-\mu}{2}>\sqrt{\mu L}$. Hence, we have:

$\displaystyle\frac{L-\mu}{2}$ $\displaystyle>\sqrt{\mu L}$ $\displaystyle
L-\mu$ $\displaystyle>2\sqrt{\mu L}$ $\displaystyle L^{2}+\mu^{2}$
$\displaystyle>6\mu L>6\mu^{2}~{}~{}(\because L>\mu)$ $\displaystyle L$
$\displaystyle>\sqrt{5}\mu.$