Year: 2020

Abstract:
Min-max formulations have attracted great attention in the ML community due to the rise of deep generative models and adversarial methods, while understanding the dynamics of gradient algorithms for solving such formulations has remained a grand challenge. As a first step, we restrict to bilinear zero-sum games and give a systematic analysis of popular gradient updates, for both simultaneous and alternating versions. We provide exact conditions for their convergence and find the optimal parameter setup and convergence rates. In particular, our results offer formal evidence that alternating updates converge “better” than simultaneous ones.

heading: 1 INTRODUCTION
text: Min-max optimization has received significant attention recently due to the popularity of generative adversarial networks (GANs) (Goodfellow et al., 2014), adversarial training (Madry et al., 2018) and reinforcement learning (Du et al., 2017; Dai et al., 2018), just to name some examples. Formally, given a bivariate function f(x,y), we aim to find a saddle point (x∗,y∗) such that f(x∗,y) ≤ f(x∗,y∗) ≤ f(x,y∗), ∀x ∈ Rn, ∀y ∈ Rn. (1.1) Since the beginning of game theory, various algorithms have been proposed for finding saddle points (Arrow et al., 1958; Dem’yanov & Pevnyi, 1972; Gol’shtein, 1972; Korpelevich, 1976; Rockafellar, 1976; Bruck, 1977; Lions, 1978; Nemirovski & Yudin, 1983; Freund & Schapire, 1999). Due to its recent resurgence in ML, new algorithms specifically designed for training GANs were proposed (Daskalakis et al., 2018; Kingma & Ba, 2015; Gidel et al., 2019b; Mescheder et al., 2017). However, due to the inherent non-convexity in deep learning formulations, our current understanding of the convergence behaviour of new and classic gradient algorithms is still quite limited, and existing analysis mostly focused on bilinear games or strongly-convex-strongly-concave games (Tseng, 1995; Daskalakis et al., 2018; Gidel et al., 2019b; Liang & Stokes, 2019; Mokhtari et al., 2019b). Nonzero-sum bilinear games, on the other hand, are known to be PPAD-complete (Chen et al., 2009) (for finding approximate Nash equilibria, see e.g. Deligkas et al. (2017)).
In this work, we study bilinear zero-sum games as a first step towards understanding general min-max optimization, although our results apply to some simple GAN settings (Gidel et al., 2019a). It is well-known that certain gradient algorithms converge linearly on bilinear zero-sum games (Liang & Stokes, 2019; Mokhtari et al., 2019b; Rockafellar, 1976; Korpelevich, 1976). These iterative algorithms usually come with two versions: Jacobi style updates or Gauss–Seidel (GS) style. In a Jacobi style, we update the two sets of parameters (i.e., x and y) simultaneously whereas in a GS style we update them alternatingly (i.e., one after the other). Thus, Jacobi style updates are naturally amenable to parallelization while GS style updates have to be sequential, although the latter is usually found to converge faster (and more stable). In numerical linear algebra, the celebrated Stein–Rosenberg theorem (Stein & Rosenberg, 1948) formally proves that in solving certain linear systems, GS updates converge strictly faster than their Jacobi counterparts, and often with a larger set of convergent instances. However, this result does not readily apply to bilinear zero-sum games.
Our main goal here is to answer the following questions about solving bilinear zero-sum games:
• When exactly does a gradient-type algorithm converge?
Contributions We summarize our main results from §3 and §4 in Table 1 and 2 respectively, with supporting experiments given in §5. We use σ1 and σn to denote the largest and the smallest singular values of matrixE (see equation 2.1), and κ := σ1/σn denotes the condition number. The algorithms will be introduced in §2. Note that we generalize gradient-type algorithms but retain the same names. Table 1 shows that in most cases that we study, whenever Jacobi updates converge, the corresponding GS updates converge as well (usually with a faster rate), but the converse is not true (§3). This extends the well-known Stein–Rosenberg theorem to bilinear games. Furthermore, Table 2 tells us that by generalizing existing gradient algorithms, we can obtain faster convergence rates.

heading: 2 PRELIMINARIES
text: In the study of GAN training, bilinear games are often regarded as an important simple example for theoretically analyzing and understanding new algorithms and techniques (e.g. Daskalakis et al., 2018; Gidel et al., 2019a;b; Liang & Stokes, 2019). It captures the difficulty in GAN training and can represent some simple GAN formulations (Arjovsky et al., 2017; Daskalakis et al., 2018; Gidel et al., 2019a; Mescheder et al., 2018). Mathematically, bilinear zero-sum games can be formulated as the following min-max problem:
minx∈Rn maxy∈Rn x >Ey + b>x+ c>y. (2.1)
The set of all saddle points (see definition in eq. (1.1)) is:
{(x,y) |Ey + b = 0, E>x+ c = 0}. (2.2)
Throughout, for simplicity we assume E to be invertible, whereas the seemingly general case with non-invertible E is treated in Appendix G. The linear terms are not essential in our analysis and we take b = c = 0 throughout the paper1. In this case, the only saddle point is (0,0). For bilinear games, it is well-known that simultaneous gradient descent ascent does not converge (Nemirovski & Yudin, 1983) and other gradient-based algorithms tailored for min-max optimization have been proposed (Korpelevich, 1976; Daskalakis et al., 2018; Gidel et al., 2019a; Mescheder et al., 2017). These iterative algorithms all belong to the class of general linear dynamical systems (LDS, a.k.a.
1If they are not zero, one can translate x and y to cancel the linear terms, see e.g. Gidel et al. (2019b).
matrix iterative processes). Using state augmentation z(t) := (x(t),y(t)) we define a general k-step LDS as follows:
z(t) = ∑k i=1Aiz (t−i) + d, (2.3)
where the matrices Ai and vector d depend on the gradient algorithm (examples can be found in Appendix C.1). Define the characteristic polynomial, withA0 = −I:
p(λ) := det( ∑k i=0Aiλ k−i). (2.4)
The following well-known result decides when such a k-step LDS converges for any initialization:
Theorem 2.1 (e.g. Gohberg et al. (1982)). The LDS in eq. (2.3) converges for any initialization (z(0), . . . ,z(k−1)) iff the spectral radius r := max{|λ| : p(λ) = 0} < 1, in which case {z(t)} converges linearly with an (asymptotic) exponent r.
Therefore, understanding the bilinear game dynamics reduces to spectral analysis. The (sufficient and necessary) convergence condition reduces to that all roots of p(λ) lie in the (open) unit disk, which can be conveniently analyzed through the celebrated Schur’s theorem (Schur, 1917):
Theorem 2.2 (Schur (1917)). The roots of a real polynomial p(λ) = a0λn + a1λn−1 + · · ·+ an are within the (open) unit disk of the complex plane iff ∀k ∈ {1, 2, . . . , n}, det(PkP>k −Q>kQk) > 0, where Pk,Qk are k × k matrices defined as: [Pk]i,j = ai−j1i≥j , [Qk]i,j = an−i+j1i≤j .
In the theorem above, we denoted 1S as the indicator function of the event S, i.e. 1S = 1 if S holds and 1S = 0 otherwise. For a nice summary of related stability tests, see Mansour (2011). We therefore define Schur stable polynomials to be those polynomials whose roots all lie within the (open) unit disk of the complex plane. Schur’s theorem has the following corollary (proof included in Appendix B.2 for the sake of completeness):
Corollary 2.1 (e.g. Mansour (2011)). A real quadratic polynomial λ2 + aλ + b is Schur stable iff b < 1, |a| < 1 + b; A real cubic polynomial λ3 + aλ2 + bλ + c is Schur stable iff |c| < 1, |a+ c| < 1 + b, b− ac < 1− c2; A real quartic polynomial λ4 + aλ3 + bλ2 + cλ+ d is Schur stable iff |c− ad| < 1− d2, |a+ c| < b+ d+ 1, and b < (1 + d) + (c− ad)(a− c)/(d− 1)2.
Let us formally define Jacobi and GS updates: Jacobi updates take the form
x(t) = T1(x (t−1),y(t−1), . . . ,x(t−k),y(t−k)), y(t) = T2(x (t−1),y(t−1), . . . ,x(t−k),y(t−k)),
while Gauss–Seidel updates replace x(t−i) with the more recent x(t−i+1) in operator T2, where T1, T2 : Rnk × Rnk → Rn can be any update functions. For LDS updates in eq. (2.3) we find a nice relation between the characteristic polynomials of Jacobi and GS updates in Theorem 2.3 (proof in Appendix B.1), which turns out to greatly simplify our subsequent analyses:
Theorem 2.3 (Jacobi vs. Gauss–Seidel). Let p(λ, γ) = det( ∑k i=0(γLi +Ui)λ
k−i), whereAi = Li + Ui and Li is strictly lower block triangular. Then, the characteristic polynomial of Jacobi updates is p(λ, 1) while that of Gauss–Seidel updates is p(λ, λ).
Compared to the Jacobi update, in some sense the Gauss–Seidel update amounts to shifting the strictly lower block triangular matrices Li one step to the left, as p(λ, λ) can be rewritten as det (∑k i=0(Li+1 +Ui)λ k−i )
, with Lk+1 := 0. This observation will significantly simplify our comparison between Jacobi and Gauss–Seidel updates.
Next, we define some popular gradient algorithms for finding saddle points in the min-max problem
min x max y f(x,y). (2.5)
We present the algorithms for a general (bivariate) function f although our main results will specialize f to the bilinear case in eq. (2.1). Note that we introduced more “step sizes” for our refined analysis, as we find that the enlarged parameter space often contains choices for faster linear convergence (see §4). We only define the Jacobi updates, while the GS counterparts can be easily inferred. We always use α1 and α2 to define step sizes (or learning rates) which are positive.
Gradient descent (GD) The generalized GD update has the following form: x(t+1) = x(t) − α1∇xf(x(t),y(t)), y(t+1) = y(t) + α2∇yf(x(t),y(t)). (2.6)
When α1 = α2, the convergence of averaged iterates (a.k.a. Cesari convergence) for convex-concave games is analyzed in (Bruck, 1977; Nemirovski & Yudin, 1978; Nedić & Ozdaglar, 2009). Recent progress on interpreting GD with dynamical systems can be seen in, e.g., Mertikopoulos et al. (2018); Bailey et al. (2019); Bailey & Piliouras (2018).
Extra-gradient (EG) We study a generalized version of EG, defined as follows: x(t+1/2) = x(t) − γ1∇xf(x(t),y(t)), y(t+1/2) = y(t) + γ2∇yf(x(t),y(t)); (2.7) x(t+1) = x(t) − α1∇xf(x(t+1/2),y(t+1/2)), y(t+1) = y(t) + α2∇yf(x(t+1/2),y(t+1/2)). (2.8)
EG was first proposed in Korpelevich (1976) with the restriction α1 = α2 = γ1 = γ2, under which linear convergence was proved for bilinear games. Convergence of EG on convex-concave games was analyzed in Nemirovski (2004); Monteiro & Svaiter (2010), and Mertikopoulos et al. (2019) provides convergence guarantees for specific non-convex-non-concave problems. For bilinear games, a slightly more generalized version was proposed in Liang & Stokes (2019) where α1 = α2, γ1 = γ2, with linear convergence proved. For later convenience we define β1 = α2γ1 and β2 = α1γ2.
Optimistic gradient descent (OGD) We study a generalized version of OGD, defined as follows: x(t+1) = x(t) − α1∇xf(x(t),y(t)) + β1∇xf(x(t−1),y(t−1)), (2.9) y(t+1) = y(t) + α2∇yf(x(t),y(t))− β2∇yf(x(t−1),y(t−1)). (2.10)
The original version of OGD was given in Popov (1980) with α1 = α2 = 2β1 = 2β2 and rediscovered in the GAN literature (Daskalakis et al., 2018). Its linear convergence for bilinear games was proved in Liang & Stokes (2019). A slightly more generalized version with α1 = α2 and β1 = β2 was analyzed in Peng et al. (2019); Mokhtari et al. (2019b), again with linear convergence proved. The stochastic case was analyzed in Hsieh et al. (2019).
Momentum method Generalized heavy ball method was analyzed in Gidel et al. (2019b): x(t+1) = x(t) − α1∇xf(x(t),y(t)) + β1(x(t) − x(t−1)), (2.11) y(t+1) = y(t) + α2∇yf(x(t),y(t)) + β2(y(t) − y(t−1)). (2.12)
This is a modification of Polyak’s heavy ball (HB) (Polyak, 1964), which also motivated Nesterov’s accelerated gradient algorithm (NAG) (Nesterov, 1983). Note that for both x-update and the y-update, we add a scale multiple of the successive difference (e.g. proxy of the momentum). For this algorithm our result below improves those obtained in Gidel et al. (2019b), as will be discussed in §3.
EG and OGD as approximations of proximal point algorithm It has been observed recently in Mokhtari et al. (2019b) that for convex-concave games, EG (α1 = α2 = γ1 = γ2 = η) and OGD (α1/2 = α2/2 = β1 = β2 = η) can be treated as approximations of the proximal point algorithm (Martinet, 1970; Rockafellar, 1976) when η is small. With this result, one can show that EG and OGD converge to saddle points sublinearly for smooth convex-concave games (Mokhtari et al., 2019a). We give a brief introduction of the proximal point algorithm in Appendix A (including a linear convergence result for the slightly generalized version).
The above algorithms, when specialized to a bilinear function f (see eq. (2.1)), can be rewritten as a 1-step or 2-step LDS (see. eq. (2.3)). See Appendix C.1 for details.

heading: 3 EXACT CONDITIONS
text: With tools from §2, we formulate necessary and sufficient conditions under which a gradient-based algorithm converges for bilinear games. We sometimes use “J” as a shorthand for Jacobi style updates and “GS” for Gauss–Seidel style updates. For each algorithm, we first write down the characteristic polynomials (see derivation in Appendix C.1) for both Jacobi and GS updates, and present the exact conditions for convergence. Specifically, we show that in many cases the GS convergence regions strictly include the Jacobi convergence regions. The proofs for Theorem 3.1, 3.2, 3.3 and 3.4 can be found in Appendix C.2, C.3, C.4, and C.5, respectively.
GD The characteristic equations can be computed as:
J: (λ− 1)2 + α1α2σ2 = 0, GS: (λ− 1)2 + α1α2σ2λ = 0. (3.1)
Scaling symmetry From section 3 we obtain a scaling symmetry (α1, α2) → (tα1, α2/t), with t > 0. With this symmetry we can always fix α1 = α2 = α. This symmetry also holds for EG and momentum. For OGD, the scaling symmetry is slightly different with (α1, β1, α2, β2) → (tα1, tβ1, α2/t, β2/t), but we can still use this symmetry to fix α1 = α2 = α.
Theorem 3.1 (GD). Jacobi GD and Gauss–Seidel GD do not converge. However, Gauss–Seidel GD can have a limit cycle while Jacobi GD always diverges.
In the constrained case, Mertikopoulos et al. (2018) and Bailey & Piliouras (2018) show that FTRL, a more generalized algorithm of GD, does not converge for polymatrix games. When α1 = α2, the result of Gauss–Seidel GD has been shown in Bailey et al. (2019).
EG The characteristic equations can be computed as:
J: (λ− 1)2 + (β1 + β2)σ2(λ− 1) + (α1α2σ2 + β1β2σ4) = 0, (3.2) GS: (λ− 1)2 + (α1α2 + β1 + β2)σ2(λ− 1) + (α1α2σ2 + β1β2σ4) = 0. (3.3)
Theorem 3.2 (EG). For generalized EG with α1 = α2 = α and γi = βi/α, Jacobi and Gauss–Seidel updates achieve linear convergence iff for any singular value σ of E, we have:
J : |β1σ2 + β2σ2 − 2| < 1 + (1− β1σ2)(1− β2σ2) + α2σ2, (1− β1σ2)(1− β2σ2) + α2σ2 < 1, (3.4) GS : |(β1 + β2 + α2)σ2 − 2| < 1 + (1− β1σ2)(1− β2σ2), (1− β1σ2)(1− β2σ2) < 1. (3.5)
If β1 +β2 +α2 < 2/σ21 , the convergence region of GS updates strictly include that of Jacobi updates.
OGD The characteristic equations can be computed as:
J: λ2(λ− 1)2 + (λα1 − β1)(λα2 − β2)σ2 = 0, (3.6) GS: λ2(λ− 1)2 + (λα1 − β1)(λα2 − β2)λσ2 = 0. (3.7)
Theorem 3.3 (OGD). For generalized OGD with α1 = α2 = α, Jacobi and Gauss–Seidel updates achieve linear convergence iff for any singular value σ of E, we have:
J : { |β1β2σ2| < 1, (α− β1)(α− β2) > 0, 4 + (α+ β1)(α+ β2)σ2 > 0, α2 ( β21σ 2 + 1 ) ( β22σ 2 + 1 ) < (β1β2σ 2 + 1)(2α(β1 + β2) + β1β2(β1β2σ 2 − 3)); (3.8)
GS :
{ (α− β1)(α− β2) > 0, (α+ β1)(α+ β2)σ2 < 4, (αβ1σ 2 + 1)(αβ2σ 2 + 1) > (1 + β1β2σ 2)2. (3.9)
The convergence region of GS updates strictly include that of Jacobi updates.
Momentum The characteristic equations can be computed as:
J: (λ− 1)2(λ− β1)(λ− β2) + α1α2σ2λ2 = 0, (3.10) GS: (λ− 1)2(λ− β1)(λ− β2) + α1α2σ2λ3 = 0. (3.11)
Theorem 3.4 (momentum). For the generalized momentum method with α1 = α2 = α, the Jacobi updates never converge, while the GS updates converge iff for any singular value σ of E, we have:
|β1β2| < 1, | − α2σ2 + β1 + β2 + 2| < β1β2 + 3, 4(β1 + 1)(β2 + 1) > α2σ2, α2σ2β1β2 < (1− β1β2)(2β1β2 − β1 − β2). (3.12)
This condition implies that at least one of β1, β2 is negative.
Prior to our work, only sufficient conditions for linear convergence were given for the usual EG and OGD; see §2 above. For the momentum method, our result improves upon Gidel et al. (2019b) where they only considered specific cases of parameters. For example, they only considered β1 = β2 ≥ −1/16 for Jacobi momentum (but with explicit rate of divergence), and β1 = −1/2, β2 = 0 for GS momentum (with convergence rate). Our Theorem 3.4 gives a more complete picture and formally justifies the necessity of negative momentum.
In the theorems above, we used the term “convergence region” to denote a subset of the parameter space (with parameters α, β or γ) where the algorithm converges. Our result shares similarity with the celebrated Stein–Rosenberg theorem (Stein & Rosenberg, 1948), which only applies to solving linear systems with non-negative matrices (if one were to apply it to our case, the matrix S in eq. (F.1) in Appendix F needs to have non-zero diagonal entries, which is not possible). In this sense, our results extend the Stein–Rosenberg theorem to cover nontrivial bilinear games.

heading: 4 OPTIMAL EXPONENTS OF LINEAR CONVERGENCE
text: In this section we study the optimal convergence rates of EG and OGD. We define the exponent of linear convergence as r = limt→∞ ||z(t)||/||z(t−1)|| which is the same as the spectral radius. For ease of presentation we fix α1 = α2 = α > 0 (using scaling symmetry) and we use r∗ to denote the optimal exponent of linear convergence (achieved by tuning the parameters α, β, γ). Our results show that by generalizing gradient algorithms one can obtain better convergence rates.
Theorem 4.1 (EG optimal). Both Jacobi and GS EG achieve the optimal exponent of linear convergence r∗ = (κ2 − 1)/(κ2 + 1) at α→ 0 and β1 = β2 = 2/(σ21 + σ2n). As κ→∞, r∗ → 1− 2/κ2.
Note that we defined βi = γiα in Section 2. In other words, we are taking very large extra-gradient steps (γi →∞) and very small gradient steps (α→ 0). Theorem 4.2 (OGD optimal). For Jacobi OGD with β1 = β2 = β, to achieve the optimal exponent of linear convergence, we must have α ≤ 2β. For the original OGD with α = 2β, the optimal exponent of linear convergence r∗ satisfies
r2∗ = 1
2 +
1
4 √
2σ21
√ (σ21 − σ2n)(5σ21 − σ2n + √ (σ21 − σ2n)(9σ21 − σ2n)), at (4.1)
β∗ = 1
4 √ 2
√ 3σ41 − (σ21 − σ2n)3/2 √ 9σ21 − σ2n + 6σ21σ2n − σ4n σ41σ 2 n . (4.2)
If κ → ∞, r∗ ∼ 1 − 1/(6κ2). For GS OGD with β2 = 0, the optimal exponent of convergence is r∗ = √ (κ2 − 1)/(κ2 + 1), at α = √ 2/σ1 and β1 = √ 2σ1/(σ 2 1 + σ 2 n). If κ→∞, r∗ ∼ 1− 1/κ2.
Remark The original OGD (Popov, 1980; Daskalakis et al., 2018) with α = 2β may not always be optimal. For example, take one-dimensional bilinear game and σ = 1, and denote the spectral radius given α, β as r(α, β). If we fix α = 1/2, by numerically solving section 3 we have
r(1/2, 1/4) ≈ 0.966, r(1/2, 1/3) ≈ 0.956, (4.3)
i.e, α = 1/2, β = 1/3 is a better choice than α = 2β = 1/2.
Numerical method We provide a numerical method for finding the optimal exponent of linear convergence, by realizing that the unit disk in Theorem 2.2 is not special. Let us call a polynomial to be r-Schur stable if all of its roots lie within an (open) disk of radius r in the complex plane. We can scale the polynomial with the following lemma:
Lemma 4.1. A polynomial p(λ) is r-Schur stable iff p(rλ) is Schur stable.
With the lemma above, one can rescale the Schur conditions and find the convergence region where the exponent of linear convergence is at most r (r < 1). A simple binary search would allow one to find a better and better convergence region. See details in Appendix D.3.

heading: 5 EXPERIMENTS
text: Bilinear game We run experiments on a simple bilinear game and choose the optimal parameters as suggested in Theorem 4.1 and 4.2. The results are shown in the left panel of Figure 1, which confirms the predicted linear rates.
Density plots We show the density plots (heat maps) of the spectral radii in Figure 2. We make plots for EG, OGD and momentum with both Jacobi and GS updates. These plots are made when β1 = β2 = β and they agree with our theorems in §3.
Wasserstein GAN As in Daskalakis et al. (2018), we consider a WGAN (Arjovsky et al., 2017) that learns the mean of a Gaussian:
minφmaxθ f(φ,θ) := Ex∼N (v,σ2I)[s(θ >x)]− Ez∼N (0,σ2I)[s(θ>(z + φ))], (5.1)
where s(x) is the sigmoid function. It can be shown that near the saddle point (θ∗,φ∗) = (0,v) the min-max optimization can be treated as a bilinear game (Appendix E.1). With GS updates, we find that Adam diverges, SGD goes around a limit cycle, and EG converges, as shown in the middle panel of Figure 1. We can see that Adam does not behave well even in this simple task of learning a single two-dimensional Gaussian with GAN.
Our next experiment shows that generalized algorithms may have an advantage over traditional ones. Inspired by Theorem 4.1, we compare the convergence of two EGs with the same parameter β = αγ, and find that with scaling, EG has better convergence, as shown in the right panel of Figure 1. Finally,
we compare Jacobi updates with GS updates. In Figure 3, we can see that GS updates converge even if the corresponding Jacobi updates do not.
Mixtures of Gaussians (GMMs) Our last experiment is on learning GMMs with a vanilla GAN (Goodfellow et al., 2014) that does not directly fall into our analysis. We choose a 3-hidden layer ReLU network for both the generator and the discriminator, and each hidden layer has 256 units. We find that for GD and OGD, Jacobi style updates converge more slowly than GS updates, and whenever Jacobi updates converge, the corresponding GS updates converges as well. These comparisons can be found in Figure 4 and 5, which implies the possibility of extending our results to non-bilinear games. Interestingly, we observe that even Jacobi GD converges on this example. We provide additional comparison between the Jacobi and GS updates of Adam (Kingma & Ba, 2015) in Appendix E.2.

heading: 6 CONCLUSIONS
text: In this work we focus on the convergence behaviour of gradient-based algorithms for solving bilinear games. By drawing a connection to discrete linear dynamical systems (§2) and using Schur’s theorem, we provide necessary and sufficient conditions for a variety of gradient algorithms, for both simultaneous (Jacobi) and alternating (Gauss–Seidel) updates. Our results show that Gauss–Seidel updates converge more easily than Jacobi updates. Furthermore, we find the optimal exponents of linear convergence for EG and OGD, and provide a numerical method for searching that exponent. We performed a number of experiments to validate our theoretical findings and suggest further analysis.
There are many future directions to explore. For example, our preliminary experiments on GANs suggest that similar (local) results might be obtained for more general games. Indeed, the local convergence behaviour of min-max nonlinear optimization can be studied through analyzing the spectrum of the Jacobian matrix of the update operator (see, e.g., Nagarajan & Kolter (2017); Gidel et al. (2019b)). We believe our framework that draws the connection to linear discrete dynamic systems and Schur’s theorem is a powerful machinery that can be applied in such problems and beyond. It would be interesting to generalize our results to the constrained case (even for bilinear games), as studied in Daskalakis & Panageas (2019); Carmon et al. (2019). Extending our results to account for stochastic noise (as empirically tested in our experiments) is another interesting direction, with results in Gidel et al. (2019a); Hsieh et al. (2019).

heading: ACKNOWLEDGEMENTS
text: We would like to thank Argyrios Deligkas, Sarath Pattathil and Georgios Piliouras for pointing out several related references. GZ is supported by David R. Cheriton Scholarship. We gratefully acknowledge funding support from NSERC and the Waterloo-Huawei Joint Innovation Lab.

heading: A PROXIMAL POINT (PP) ALGORITHM
text: PP was originally proposed by Martinet (1970) with α1 = α2 and then carefully studied by Rockafellar (1976). The linear convergence for bilinear games was also proved in the same reference. Note that we do not consider Gauss–Seidel PP since we do not get a meaningful solution after a shift of steps2.
x(t+1) = x(t) − α1∇xf(x(t+1),y(t+1)), y(t+1) = y(t) + α2∇yf(x(t+1),y(t+1)), (A.1)
where x(t+1) and y(t+1) are given implicitly by solving the equations above. For bilinear games, one can derive that:
z(t+1) =
[ I α1E
−α2E> I
]−1 z(t). (A.2)
We can compute the exact form of the inverse matrix, but perhaps an easier way is just to compute the spectrum of the original matrix (the same as Jacobi GD except that we flip the signs of αi) and perform λ→ 1/λ. Using the fact that the eigenvalues of a matrix are reciprocals of the eigenvalues of its inverse, the characteristic equation is:
(1/λ− 1)2 + α1α2σ2 = 0. (A.3)
With the scaling symmetry (α1, α2) → (tα1, α2/t), we can take α1 = α2 = α > 0. With the notations in Corollary 2.1, we have a = −2/(1 + α2σ2) and b = 1/(1 + α2σ2), and it is easy to check |a| < 1 + b and b < 1 are always satisfied, which means linear convergence is always guaranteed. Hence, we have the following theorem: Theorem A.1. For bilinear games, the proximal point algorithm always converges linearly.
Although the proximal point algorithm behaves well, it is rarely used in practice since it is an implicit method, i.e., one needs to solve (x(t+1),y(t+1)) from equation A.1.

heading: B PROOFS IN SECTION 2

heading: B.1 PROOF OF THEOREM 2.3
text: In this section we apply Theorem 2.1 to prove Theorem 2.3, an interesting connection between Jacobi and Gauss–Seidel updates:
Theorem 2.3 (Jacobi vs. Gauss–Seidel). Let p(λ, γ) = det( ∑k i=0(γLi +Ui)λ
k−i), whereAi = Li + Ui and Li is strictly lower block triangular. Then, the characteristic polynomial of Jacobi updates is p(λ, 1) while that of Gauss–Seidel updates is p(λ, λ).
Let us first consider the block linear iterative process in the sense of Jacobi (i.e., all blocks are updated simultaneously):
z(t) = z (t) 1 ... z (t) b  = k∑ i=1 Ai z (t−i) 1 ... z (t−i) b  = k∑ i=1  l−1∑ j=1 Ai,jz (t−i) j + b∑ j=l Ai,jz (t−i) j + d, (B.1) whereAi,j is the j-th column block ofAi. For each matrixAi, we decompose it into the sum
Ai = Li +Ui, (B.2)
where Li is the strictly lower block triangular part and Ui is the upper (including diagonal) block triangular part. Theorem 2.1 indicates that the convergence behaviour of equation B.1 is governed by the largest modulus of the roots of the characteristic polynomial:
det ( −λkI +
k∑ i=1 Aiλ k−i
) = det ( −λkI +
k∑ i=1 (Li +Ui)λ k−i
) . (B.3)
2If one uses inverse operators this is in principle doable.
Alternatively, we can also consider the updates in the sense of Gauss–Seidel (i.e., blocks are updated sequentially):
z (t) l = k∑ i=1  l−1∑ j=1 Ai,jz (t−i+1) j + b∑ j=l Ai,jz (t−i) j  l + dl, l = 1, . . . , b. (B.4)
We can rewrite the Gauss–Seidel update elegantly3 as:
(I −L1)z(t) = k∑ i=1 (Li+1 +Ui)z (t−i) + d, (B.5)
i.e.,
z(t) = k∑ i=1 (I −L1)−1(Li+1 +Ui)z(t−i) + (I −L1)−1d, (B.6)
where Lk+1 := 0. Applying Theorem 2.1 again we know the convergence behaviour of the Gauss– Seidel update is governed by the largest modulus of roots of the characteristic polynomial:
det ( −λkI +
k∑ i=1
(I −L1)−1(Li+1 +Ui)λk−i )
(B.7)
= det ( (I −L1)−1 ( − λkI + λkL1 +
k∑ i=1
(Li+1 +Ui)λ k−i ))
(B.8)
= det(I −L1)−1 · det ( k∑ i=0 (λLi +Ui)λ k−i ) (B.9)
Note thatA0 = −I and the factor det(I −L1)−1 can be discarded since multiplying a characteristic polynomial by a non-zero constant factor does not change its roots.

heading: B.2 PROOF OF COROLLARY 2.1
text: Corollary 2.1 (e.g. Mansour (2011)). A real quadratic polynomial λ2 + aλ + b is Schur stable iff b < 1, |a| < 1 + b; A real cubic polynomial λ3 + aλ2 + bλ + c is Schur stable iff |c| < 1, |a+ c| < 1 + b, b− ac < 1− c2; A real quartic polynomial λ4 + aλ3 + bλ2 + cλ+ d is Schur stable iff |c− ad| < 1− d2, |a+ c| < b+ d+ 1, and b < (1 + d) + (c− ad)(a− c)/(d− 1)2.
Proof. It suffices to prove the result for quartic polynomials. We write down the matrices:
P1 = [1], Q1 = [d], (B.10)
P2 = [ 1 0 a 1 ] , Q2 = [ d c 0 d ] , (B.11)
P3 = [ 1 0 0 a 1 0 b a 1 ] ,Q3 = [ d c b 0 d c 0 0 d ] , (B.12)
P4 = 1 0 0 0a 1 0 0b a 1 0 c b a 0  ,Q4 = d c b a0 d c b0 0 d c 0 0 0 d  . (B.13) We require det(PkP>k − Q>kQk) =: δk > 0, for k = 1, 2, 3, 4. If k = 1, we have 1 − d2 > 0, namely, |d| < 1. δ2 > 0 reduces to (c − ad)2 < (1 − d2)2 and thus |c − ad| < 1 − d2 due to the first condition. δ4 > 0 simplifies to:
−((a+ c)2 − (b+ d+ 1)2)((b− d− 1)(d− 1)2 − (a− c)(c− ad))2 < 0, (B.14) 3This is well-known when k = 1, see e.g. Saad (2003).
which yields |a+ c| < |b+ d+ 1|. Finally, δ3 > 0 reduces to:
((b− d− 1)(d− 1)2 − (a− c)(c− ad))((d2 − 1)(b+ d+ 1) + (c− ad)(a+ c)) > 0. (B.15)
Denote p(λ) := λ4 + aλ3 + bλ2 + cλ + d, we must have p(1) > 0 and p(−1) > 0, as otherwise there is a real root λ0 with |λ0| ≥ 1. Hence we obtain b + d + 1 > |a + c| > 0. Also, from |c− ad| < 1− d2, we know that:
|c− ad| · |a+ c| < |b+ d+ 1|(1− d2) = (b+ d+ 1)(1− d2). (B.16)
So, the second factor in B.15 is negative and the positivity of the first factor reduces to:
b < (1 + d) + (c− ad)(a− c)
(d− 1)2 . (B.17)
To obtain the Schur condition for cubic polynomials, we take d = 0, and the quartic Schur condition becomes:
|c| < 1, |a+ c| < b+ 1, b− ac < 1− c2. (B.18)
To obtain the Schur condition for quadratic polynomials, we take c = 0 in the above and write:
b < 1, |a| < 1 + b. (B.19)
The proof is now complete.

heading: C PROOFS IN SECTION 3
text: Some of the following proofs in Appendix C.4 and C.5 rely on Mathematica code (mostly with the built-in function Reduce) but in principle the code can be verified manually using cylindrical algebraic decomposition.4

heading: C.1 DERIVATION OF CHARACTERISTIC POLYNOMIALS
text: In this appendix, we derive the exact forms of LDSs (eq. (2.3)) and the characteristic polynomials for all gradient-based methods introduced in §2, with eq. (2.4). The following lemma is well-known and easy to verify using Schur’s complement:
Lemma C.1. GivenM ∈ R2n×2n,A ∈ Rn×n and
M = [ A B C D ] . (C.1)
If C andD commute, then detM = det(AD −BC).
Gradient descent From equation 2.6 the update equation of Jacobi GD can be derived as:
z(t+1) =
[ I −α1E
α2E > I
] z(t), (C.2)
and with Lemma C.1, we compute the characteristic polynomial as in eq. (2.4):
det [ (λ− 1)I α1E −α2E> (λ− 1)I ] = det[(λ− 1)2I + α1α2EE>], (C.3)
With spectral decomposition we obtain equation 3.1. Taking α2 → λα2 and with Theorem 2.3 we obtain the corresponding GS updates. Therefore, the characteristic polynomials for GD are:
J: (λ− 1)2 + α1α2σ2 = 0, GS: (λ− 1)2 + α1α2σ2λ = 0. (C.4) 4See the online Mathematica documentation.
Extra-gradient From eq. (2.7) and eq. (2.8), the update of Jacobi EG is:
z(t+1) = [ I − β2EE> −α1E α2E > I − β1E>E ] z(t), (C.5)
the characteristic polynomial is:
det
[ (λ− 1)I + β2EE> α1E
−α2E> (λ− 1)I + β1E>E
] . (C.6)
Since we assumed α2 > 0, we can left multiply the second row by β2E/α2 and add it to the first row. Hence, we obtain:
det [ (λ− 1)I α1E + (λ− 1)β2E/α2 + β1β2EE>E/α2 −α2E> (λ− 1)I + β1E>E ] . (C.7)
With Lemma C.1 the equation above becomes:
det[(λ− 1)2I + (β1 + β2)E>E(λ− 1) + (α1α2E>E + β1β2E>EE>E)], (C.8) which simplifies to equation 3.2 with spectral decomposition. Note that to obtain the GS polynomial, we simply take α2 → λα2 in the Jacobi polynomial as shown in Theorem 2.3. For the ease of reading we copy the characteristic equations for generalized EG:
J: (λ− 1)2 + (β1 + β2)σ2(λ− 1) + (α1α2σ2 + β1β2σ4) = 0, (C.9) GS: (λ− 1)2 + (α1α2 + β1 + β2)σ2(λ− 1) + (α1α2σ2 + β1β2σ4) = 0. (C.10)
Optimistic gradient descent We can compute the LDS for OGD with eq. (2.9) and eq. (2.10):
z(t+2) =
[ I −α1E
α2E > I
] z(t+1) + [ 0 β1E
−β2E> 0
] z(t), (C.11)
With eq. (2.4), the characteristic polynomial for Jacobi OGD is
det
[ (λ2 − λ)I (λα1 − β1)E
(−λα2 + β2)E> (λ2 − λ)I
] . (C.12)
Taking the determinant and with Lemma C.1 we obtain equation 3.6. The characteristic polynomial for GS updates in equation 3.7 can be subsequently derived with Theorem 2.3, by taking (α2, β2)→ (λα2, λβ2). For the ease of reading we copy the characteristic polynomials from the main text as:
J: λ2(λ− 1)2 + (λα1 − β1)(λα2 − β2)σ2 = 0, (C.13) GS: λ2(λ− 1)2 + (λα1 − β1)(λα2 − β2)λσ2 = 0. (C.14)
Momentum method With eq. (2.11) and eq. (2.12), the LDS for the momentum method is:
z(t+2) = [ (1 + β1)I −α1E α2E > (1 + β2)I ] z(t+1) + [ −β1I 0 0 −β2I ] z(t), (C.15)
From eq. (2.4), the characteristic polynomial for Jacobi momentum is
det
[ (λ2 − λ(1 + β1) + β1)I λα1E
−λα2E> (λ2 − λ(1 + β2) + β2)I
] . (C.16)
Taking the determinant and with Lemma C.1 we obtain equation 3.10, while equation 3.11 can be derived with Theorem 2.3, by taking α2 → λα2. For the ease of reading we copy the characteristic polynomials from the main text as:
J: (λ− 1)2(λ− β1)(λ− β2) + α1α2σ2λ2 = 0, (C.17) GS: (λ− 1)2(λ− β1)(λ− β2) + α1α2σ2λ3 = 0. (C.18)

heading: C.2 PROOF OF THEOREM 3.1: SCHUR CONDITIONS OF GD
text: Theorem 3.1 (GD). Jacobi GD and Gauss–Seidel GD do not converge. However, Gauss–Seidel GD can have a limit cycle while Jacobi GD always diverges.
Proof. With the notations in Corollary 2.1, for Jacobi GD, b = 1 + α2σ2 > 1. For Gauss–Seidel GD, b = 1. The Schur conditions are violated.

heading: C.3 PROOF OF THEOREM 3.2: SCHUR CONDITIONS OF EG
text: Theorem 3.2 (EG). For generalized EG with α1 = α2 = α and γi = βi/α, Jacobi and Gauss–Seidel updates achieve linear convergence iff for any singular value σ of E, we have:
J : |β1σ2 + β2σ2 − 2| < 1 + (1− β1σ2)(1− β2σ2) + α2σ2, (1− β1σ2)(1− β2σ2) + α2σ2 < 1, (3.4) GS : |(β1 + β2 + α2)σ2 − 2| < 1 + (1− β1σ2)(1− β2σ2), (1− β1σ2)(1− β2σ2) < 1. (3.5)
If β1 +β2 +α2 < 2/σ21 , the convergence region of GS updates strictly include that of Jacobi updates.
Both characteristic polynomials can be written as a quadratic polynomial λ2 + aλ+ b, where:
J: a = (β1 + β2)σ2 − 2, b = (1− β1σ2)(1− β2σ2) + α2σ2, (C.19) GS: a = (β1 + β2 + α2)σ2 − 2, b = (1− β1σ2)(1− β2σ2). (C.20)
Compared to Jacobi EG, the only difference between Gauss–Seidel and Jacobi updates is that the α2σ2 in b is now in a, which agrees with Theorem 2.3. Using Corollary 2.1, we can derive the Schur conditions equation 3.4 and equation 3.5.
More can be said if β1 +β2 is small. For instance, if β1 +β2 +α2 < 2/σ21 , then equation 3.4 implies equation 3.5. In this case, the first conditions of equation 3.4 and equation 3.5 are equivalent, while the second condition of equation 3.4 strictly implies that of equation 3.5. Hence, the Schur region of Gauss–Seidel updates includes that of Jacobi updates. The same holds true if β1 + β2 < 43σ21 .
More precisely, to show that the GS convergence region strictly contains that of the Jacobi convergence region, simply take β1 = β2 = β. The Schur condition for Jacobi EG and Gauss–Seidel EG are separately:
J: α2σ2 + (βσ2 − 1)2 < 1, (C.21) GS: 0 < βσ2 < 2 and |ασ| < 2− βσ2. (C.22)
It can be shown that if β = α2/3 and α→ 0, equation C.21 is always violated whereas equation C.22 is always satisfied.
Conversely, we give an example when Jacobi EG converges while GS EG does not. Let β1σ2 = β2σ 2 ≡ 32 , then Jacobi EG converges iff α 2σ2 < 34 while GS EG converges iff α 2σ2 < 14 .

heading: C.4 PROOF OF THEOREM 3.3: SCHUR CONDITIONS OF OGD
text: In this subsection, we fill in the details of the proof of Theorem 3.3, by first deriving the Schur conditions of OGD, and then studying the relation between Jacobi OGD and GS OGD. Theorem 3.3 (OGD). For generalized OGD with α1 = α2 = α, Jacobi and Gauss–Seidel updates achieve linear convergence iff for any singular value σ of E, we have:
J : { |β1β2σ2| < 1, (α− β1)(α− β2) > 0, 4 + (α+ β1)(α+ β2)σ2 > 0, α2 ( β21σ 2 + 1 ) ( β22σ 2 + 1 ) < (β1β2σ 2 + 1)(2α(β1 + β2) + β1β2(β1β2σ 2 − 3)); (3.8)
GS :
{ (α− β1)(α− β2) > 0, (α+ β1)(α+ β2)σ2 < 4, (αβ1σ 2 + 1)(αβ2σ 2 + 1) > (1 + β1β2σ 2)2. (3.9)
The convergence region of GS updates strictly include that of Jacobi updates.
The Jacobi characteristic polynomial is now quartic in the form λ4 + aλ3 + bλ2 + cλ+ d, with
a = −2, b = α2σ2 + 1, c = −α(β1 + β2)σ2, d = β1β2σ2. (C.23) Comparably, the GS polynomial equation 3.7 can be reduced to a cubic one λ3 + aλ2 + bλ+ c with
a = −2 + α2σ2, b = −α(β1 + β2)σ2 + 1, c = β1β2σ2. (C.24)
First we derive the Schur conditions equation 3.8 and equation 3.9. Note that other than Corollary 2.1, an equivalent Schur condition can be read from Cheng & Chiou (2007, Theorem 1) as:
Theorem C.1 (Cheng & Chiou (2007)). A real quartic polynomial λ4 + aλ3 + bλ2 + cλ + d is Schur stable iff:
|d| < 1, |a| < d+ 3, |a+ c| < b+ d+ 1, (1− d)2b+ c2 − a(1 + d)c− (1 + d)(1− d)2 + a2d < 0. (C.25)
With equation C.23 and Theorem C.1, it is straightforward to derive equation 3.8. With equation C.24 and Corollary 2.1, we can derive equation 3.9 without much effort.
Now, let us study the relation between the convergence region of Jacobi OGD and GS OGD, as given in equation 3.8 and equation 3.9. Namely, we want to prove the last sentence of Theorem 3.3. The outline of our proof is as follows. We first show that each region of (α, β1, β2) described in equation 3.8 (the Jacobi region) is contained in the region described in equation 3.9 (the GS region). Since we are only studying one singular value, we slightly abuse the notations and rewrite βiσ as βi (i = 1, 2) and ασ as α. From equation 3.6 and equation 3.7, β1 and β2 can switch. WLOG, we assume β1 ≥ β2. There are four cases to consider:
• β1 ≥ β2 > 0. The third Jacobi condition in equation 3.8 now is redundant, and we have α > β1 or α < β2 for both methods. Solving the quadratic feasibility condition for α gives:
0 < β2 < 1, β2 ≤ β1 < β2 +
√ 4 + 5β22
2(1 + β22) , β1 < α <
u+ √ u2 + tv
t , (C.26)
where u = (β1β2 + 1)(β1 + β2), v = β1β2(β1β2 + 1)(β1β2 − 3), t = (β21 + 1)(β22 + 1). On the other hand, assume α > β1, the first and third GS conditions are automatic. Solving the second gives:
0 < β2 < 1, β2 ≤ β1 < −β2 +
√ 8 + β22 2 , β1 < α < − 1 2 (β1+β2)+ 1 2 √ (β1 − β2)2 + 16.
(C.27) Define f(β2) := −β2 + √ 8 + β22/2 and g(β2) := (β2 + √ 4 + 5β22)/(2(1 + β 2 2)), and one can show that
f(β2) ≥ g(β2). (C.28)
Furthermore, it can also be shown that given 0 < β2 < 1 and β2 ≤ β1 < g(β2), we have (u+ √ u2 + 4v)/t < −(β1 + β2)/2 + (1/2) √ (β1 − β2)2 + 16. (C.29)
• β1 ≥ β2 = 0. The Schur condition for Jacobi and Gauss–Seidel updates reduces to:
Jacobi: 0 < β1 < 1, β1 < α < 2β1
1 + β21 , (C.30)
GS: 0 < β1 < √ 2, β1 < α < −β1 +
√ 16 + β21
2 . (C.31)
One can show that given β1 ∈ (0, 1), we have 2β1/(1 + β21) < (−β1 + √
16 + β21)/2. • β1 ≥ 0 > β2. Reducing the first, second and fourth conditions of equation 3.8 yields:
β2 < 0, 0 < β1 < β2 +
√ 4 + 5β22
2(1 + β22) , β1 < α <
u+ √ u2 + tv
t . (C.32)
This region contains the Jacobi region. It can be similarly proved that even within this larger region, GS Schur condition equation 3.9 is always satisfied. • β2 ≤ β1 < 0. We have u < 0, tv < 0 and thus α < (u + √ u2 + tv)/t < 0. This
contradicts our assumption that α > 0.
Combining the four cases above, we know that the Jacobi region is contained in the GS region.
To show the strict inclusion, take β1 = β2 = α/5 and α → 0. One can show that as long as α is small enough, all the Jacobi regions do not contain this point, each of which is described with a
singular value in equation 3.8. However, all the GS regions described in equation 3.9 contain this point.
The proof above is still missing some details. We provide the proofs of equation C.26, equation C.28, equation C.29 and equation C.32 in the sub-sub-sections below, with the help of Mathematica, although one can also verify these claims manually. Moreover, a one line proof of the inclusion can be given with Mathematica code, as shown in Section C.4.5.

heading: C.4.1 PROOF OF EQUATION C.26
text: The fourth condition of equation 3.8 can be rewritten as:
α2t− 2uα− v < 0, (C.33)
where u = (β1β2 + 1)(β1 + β2), v = β1β2(β1β2 + 1)(β1β2 − 3), t = (β21 + 1)(β22 + 1). The discriminant is 4(u2 + tv) = (1 − β1β2)2(1 + β1β2)(β21 + β22 + β21β22 − β1β2) ≥ 0. Since if β1β2 < 0,
β21 + β 2 2 + β 2 1β 2 2 − β1β2 = β21 + β22 + β1β2(β1β2 − 1) > 0,
If β1β2 ≥ 0, β21 + β 2 2 + β 2 1β 2 2 − β1β2 = (β1 − β2)2 + β1β2(1 + β1β2) ≥ 0, where we used |β1β2| < 1 in both cases. So, equation C.33 becomes:
u− √ u2 + tv
t < α <
u+ √ u2 + tv
t . (C.34)
Combining with α > β1 or α < β2 obtained from the second condition, we have:
u− √ u2 + tv
t < α < β2 or β1 < α <
u+ √ u2 + tv
t . (C.35)
The first case is not possible, with the following code:
u = (b1 b2 + 1) (b1 + b2); v = b1 b2 (b1 b2 + 1) (b1 b2 - 3); t = (b1^2 + 1) (b2^2 + 1); Reduce[b2 t > u - Sqrt[u^2 + t v] && b1 >= b2 > 0 && Abs[b1 b2] < 1],
and we have:
False.
Therefore, the only possible case is β1 < α < (u+ √ u2 + tv)/t. Where the feasibility region can be solved with:
Reduce[b1 t < u + Sqrt[u^2+t v]&&b1>=b2>0&&Abs[b1 b2] < 1].
What we get is:
0<b2<1 && b2<=b1<b2/(2 (1+b2^2))+1/2 Sqrt[(4+5 b2^2)/(1+b2^2)^2].
Therefore, we have proved equation C.26.

heading: C.4.2 PROOF OF EQUATION C.28
text: With
Reduce[-(b2/2) + Sqrt[8 + b2^2]/2 >= (b2 + Sqrt[4 + 5 b2^2])/(2 (1 + b2^2)) && 0 < b2 < 1],
we can remove the first constraint and get:
0 < b2 < 1.

heading: C.4.3 PROOF OF EQUATION C.29
text: Given
Reduce[-1/2 (b1 + b2) + 1/2 Sqrt[(b1 - b2)^2 + 16] > (u + Sqrt[u^2 + t v])/t &&
0 < b2 < 1 && b2 <= b1 < (b2 + Sqrt[4 + 5 b2^2])/(2 (1 + b2^2)), {b2, b1}],
we can remove the first constraint and get:
0 < b2 < 1 && b2 <= b1 < b2/(2 (1 + b2^2)) + 1/2 Sqrt[(4 + 5 b2^2)/(1 + b2^2)^2].

heading: C.4.4 PROOF OF EQUATION C.32
text: The second Jacobi condition simplifies to α > β1 and the fourth simplifies to equation C.34. Combining with the first Jacobi condition:
Reduce[Abs[b1 b2] < 1 && a > b1 && (u - Sqrt[u^2 + t v])/t < a < (u + Sqrt[u^2 + t v])/t && b1 >= 0 && b2 < 0, {b2, b1, a} ] // Simplify,
we have:
b2 < 0 && b1 > 0 && b2/(1 + b2^2) + Sqrt[(4 + 5 b2^2)/(1 + b2^2)^2] > 2 b1 && b1 < a < (b1 + b2 + b1^2 b2 + b1 b2^2)/((1 + b1^2) (1 + b2^2)) +
Sqrt[((-1 + b1 b2)^2 (b1^2 + b2^2 + b1 b2 (-1 + b2^2) + b1^3 (b2 + b2^3)))/((1 + b1^2)^2 (1 + b2^2)^2)].
This can be further simplified to achieve equation C.32.

heading: C.4.5 ONE LINE PROOF
text: In fact, there is another very simple proof:
Reduce[ForAll[{b1, b2, a}, (a - b1) (a - b2) > 0 && (a + b1) (a + b2) > -4 && Abs[b1 b2] < 1 && a^2 (b1^2 + 1) (b2^2 + 1) < (b1 b2 + 1) (2 a (b1 + b2) + b1 b2 (b1 b2 - 3)), (a - b1) (a - b2) > 0 && (a + b1) (a + b2) < 4 && (a b1 + 1) (a b2 + 1) > (1 + b1 b2)^2], {b2, b1, a}] True.
However, this proof does not tell us much information about the range of our variables.

heading: C.5 PROOF OF THEOREM 3.4: SCHUR CONDITIONS OF MOMENTUM
text: Theorem 3.4 (momentum). For the generalized momentum method with α1 = α2 = α, the Jacobi updates never converge, while the GS updates converge iff for any singular value σ of E, we have:
|β1β2| < 1, | − α2σ2 + β1 + β2 + 2| < β1β2 + 3, 4(β1 + 1)(β2 + 1) > α2σ2, α2σ2β1β2 < (1− β1β2)(2β1β2 − β1 − β2). (3.12)
This condition implies that at least one of β1, β2 is negative.

heading: C.5.1 SCHUR CONDITIONS OF JACOBI AND GS UPDATES
text: Jacobi condition We first rename ασ as al and β1, β2 as b1, b2. With Theorem C.1:
{Abs[d] < 1, Abs[a] < d + 3, a + b + c + d + 1 > 0, -a + b - c + d + 1 > 0, (1 - d)^2 b - (c - a d) (a - c) - (1 + d) (1 - d)^2 < 0} /. {a -> -2 - b1 - b2, b -> al^2 + 1 + 2 (b1 + b2) + b1 b2, c -> -b1 - b2 - 2 b1 b2, d -> b1 b2} // FullSimplify.
We obtain:
{Abs[b1 b2] < 1, Abs[2 + b1 + b2] < 3 + b1 b2, al^2 > 0, al^2 + 4 (1 + b1) (1 + b2) > 0, al^2 (-1 + b1 b2)^2 < 0}.
The last condition is never satisfied and thus Jacobi momentum never converges.
Gauss–Seidel condition With Theorem C.1, we compute:
{Abs[d] < 1, Abs[a] < d + 3, a + b + c + d + 1 > 0, -a + b - c + d + 1 > 0, (1 - d)^2 b + c^2 - a (1 + d) c - (1 + d) (1 - d)^2 + a^2 d < 0} /. {a -> al^2 - 2 - b1 - b2, b -> 1 + 2 (b1 + b2) + b1 b2, c -> -b1 - b2 - 2 b1 b2, d -> b1 b2} // FullSimplify.
The result is:
{Abs[b1 b2] < 1, Abs[2 - al^2 + b1 + b2] < 3 + b1 b2, al^2 > 0, 4 (1 + b1) (1 + b2) > al^2, al^2 (b1 + b2 + (-2 + al^2 - b1) b1 b2 + b1 (-1 + 2 b1) b2^2) < 0},
which can be further simplified to equation 3.12.

heading: C.5.2 NEGATIVE MOMENTUM
text: With Theorem 3.4, we can actually show that in general at least one of β1 and β2 must be negative. There are three cases to consider, and in each case we simplify equation 3.12:
1. β1β2 = 0. WLOG, let β2 = 0, and we obtain
−1 < β1 < 0 and α2σ2 < 4(1 + β1). (C.36)
2. β1β2 > 0. We have
−1 < β1 < 0, −1 < β2 < 0 , α2σ2 < 4(1 + β1)(1 + β2). (C.37)
3. β1β2 < 0. WLOG, we assume β1 ≥ β2. We obtain: −1 < β2 < 0, 0 < β1 < min { − 1 3β2 , ∣∣∣− β2 1 + 2β2 ∣∣∣} . (C.38) The constraints for α are α > 0 and:
max
{ (1− β1β2)(2β1β2 − β1 − β2)
β1β2 , 0
} < α2σ2 < 4(1 + β1)(1 + β2). (C.39)
These conditions can be further simplified by analyzing all singular values. They only depend on σ1 and σn, the largest and the smallest singular values. Now, let us derive equation C.37, equation C.38 and equation C.39 more carefully. Note that we use a for ασ.

heading: C.5.3 PROOF OF EQUATION C.37
text: Reduce[Abs[b1 b2] < 1 && Abs[-a^2 + b1 + b2 + 2] < b1 b2 + 3 && 4 (b1 + 1) (b2 + 1) > a^2 && a^2 b1 b2 < (1 - b1 b2) (2 b1 b2 - b1 - b2) && b1 b2 > 0 && a > 0, {b2, b1, a}]
-1 < b2 < 0 && -1 < b1 < 0 && 0 < a < Sqrt[4 + 4 b1 + 4 b2 + 4 b1 b2]

heading: C.5.4 PROOF OF EQUATIONS C.38 AND C.39
text: Reduce[Abs[b1 b2] < 1 && Abs[-a^2 + b1 + b2 + 2] < b1 b2 + 3 && 4 (b1 + 1) (b2 + 1) > a^2 && a^2 b1 b2 < (1 - b1 b2) (2 b1 b2 - b1 - b2) && b1 b2 < 0 && b1 >= b2 && a > 0, {b2, b1, a}]
(-1 < b2 <= -(1/3) && ((0 < b1 <= b2/(-1 + 2 b2) && 0 < a < Sqrt[4 + 4 b1 + 4 b2 + 4 b1 b2]) || (b2/(-1 + 2 b2) < b1 < -(1/(3 b2)) &&
Sqrt[(-b1 - b2 + 2 b1 b2 + b1^2 b2 + b1 b2^2 - 2 b1^2 b2^2)/( b1 b2)] < a < Sqrt[4 + 4 b1 + 4 b2 + 4 b1 b2]))) || (-(1/3) <
b2 < 0 && ((0 < b1 <= b2/(-1 + 2 b2) && 0 < a < Sqrt[4 + 4 b1 + 4 b2 + 4 b1 b2]) || (b2/(-1 + 2 b2) < b1 < -(b2/(1 + 2 b2)) &&
Sqrt[(-b1 - b2 + 2 b1 b2 + b1^2 b2 + b1 b2^2 - 2 b1^2 b2^2)/( b1 b2)] < a < Sqrt[4 + 4 b1 + 4 b2 + 4 b1 b2])))
Some further simplication yields equation C.38 and equation C.39.

heading: D PROOFS IN SECTION 4
text: For bilinear games and gradient-based methods, a Schur condition defines the region of convergence in the parameter space, as we have seen in Section 3. However, it is unknown which setting of parameters has the best convergence rate in a Schur stable region. We explore this problem now. Due to Theorem 3.1, we do not need to study GD. The remaining cases are EG, OGD and GS momentum (Jacobi momentum does not converge due to Theorem 3.4). Analytically (Section D.1 and D.2), we study the optimal linear rates for EG and special cases of generalized OGD (Jacobi OGD with β1 = β2 and Gauss–Seidel OGD with β2 = 0). The special cases include the original form of OGD. We also provide details for the numerical method described at the end of Section 4.
The optimal spectral radius is obtained by solving another min-max optimization problem:
min θ max σ∈Sv(E) r(θ, σ), (D.1)
where θ denotes the collection of all hyper-parameters, and r(θ, σ) is defined as the spectral radius function that relies on the choice of parameters and the singular value σ. We also use Sv(E) to denote the set of singular values of E.
In general, the function r(θ, σ) is non-convex and thus difficult to analyze. However, in the special case of quadratic characteristic polynomials, it is possible to solve equation D.1. This is how we will analyze EG and special cases of OGD, as r(θ, σ) can be expressed using root functions of quadratic polynomials. For cubic and quartic polynomials, it is in principle also doable as we have analytic formulas for the roots. However, these formulas are extremely complicated and difficult to optimize and we leave it for future work. For EG and OGD, we will show that the optimal linear rates depend only on the conditional number κ := σ1/σn.
For simplicity, we always fix α1 = α2 = α > 0 using the scaling symmetry studied in Section 3.

heading: D.1 PROOF OF THEOREM 4.1: OPTIMAL CONVERGENCE RATE OF EG
text: Theorem 4.1 (EG optimal). Both Jacobi and GS EG achieve the optimal exponent of linear convergence r∗ = (κ2 − 1)/(κ2 + 1) at α→ 0 and β1 = β2 = 2/(σ21 + σ2n). As κ→∞, r∗ → 1− 2/κ2.

heading: D.1.1 JACOBI EG
text: For Jacobi updates, if β1 = β2 = β, by solving the roots of equation 3.2, the min-max problem is:
min α,β max σ∈Sv(E)
√ α2σ2 + (1− βσ2)2. (D.2)
If σ1 = σn = σ, we can simply take α→ 0 and β = 1/σ2 to obtain a super-linear convergence rate. Otherwise, let us assume σ1 > σn. We obtain a lower bound by taking α → 0 and equation D.2 reduces to:
min β max σ∈Sv(E)
|1− βσ2|. (D.3)
The optimal solution is given at 1− βσ2n = βσ21 − 1, yielding β = 2/(σ21 + σ2n). The optimal radius is thus (σ21 − σ2n)/(σ21 + σ2n) since the lower bound equation D.3 can be achieved by taking α→ 0. From general β1, β2, it can be verified that the optimal radius is achieved at β1 = β2 and the problem reduces to the previous case. The optimization problem is:
min α,β1,β2 max σ∈Sv(E) r(α, β1, β2, σ), (D.4)
where
r(α, β1, β2, σ) = {√ (1− β1σ2)(1− β2σ2) + α2σ2 4α2 > (β1 − β2)2σ2, |1− 12 (β1 + β2)σ 2|+ 12 √ (β1 − β2)2σ4 − 4α2σ2 4α2 ≤ (β1 − β2)2σ2.
In the first case, a lower bound is obtained at α2 = (β1 − β2)2σ2/4 and thus the objective only depends on β1 + β2. In the second case, the lower bound is obtained at α → 0 and β1 → β2. Therefore, the function is optimized at β1 = β2 and α→ 0. Our analysis above does not mean that α → 0 and β1 = β2 = 2/(σ21 + σ2n) is the only optimal choice. For example, when σ1 = σn = 1, we can take β1 = 1 + α and β2 = 1 − α to obtain a super-linear convergence rate.

heading: D.1.2 GAUSS–SEIDEL EG
text: For Gauss–Seidel updates and β1 = β2 = β, we do the following optimization:
min α,β max σ∈Sv(E) r(α, β, σ), (D.5)
where by solving equation 3.3:
r(α, β, σ) = { 1− βσ2 α2σ2 < 4(1− βσ2), α2
2 σ 2 − (1− βσ2) +
√ α2σ2(α2σ2 − 4(1− βσ2))/2 α2σ2 ≥ 4(1− βσ2).
r(σ, β, σ2) is quasi-convex in σ2, so we just need to minimize over α, β at both end points. Hence, equation D.5 reduces to:
min α,β
max{r(α, β, σ1), r(α, β, σn)}.
By arguing over three cases: α2 + 4β < 4/σ21 , α 2 + 4β > 4/σ2n and 4/σ 2 1 ≤ α2 + 4β ≤ 4/σ2n, we find that the minimum (κ2 − 1)/(κ2 + 1) can be achieved at α→ 0 and β = 2/(σ21 + σ2n), the same as Jacobi EG. This is because α→ 0 decouples x and y and it does not matter whether the update is Jacobi or GS.
For general β1, β2, it can be verified that the optimal radius is achieved at β1 = β2. We do the following transformation: βi → ξi − α2/2, so that the characteristic polynomial becomes:
(λ− 1)2 + (ξ1 + ξ2)σ2(λ− 1) + α2σ2 + (ξ1 − α2/2)(ξ2 − α2/2)σ4 = 0. (D.6) Denote ξ1 + ξ2 = φ, and (ξ1 − α2/2)(ξ2 − α2/2) = ν, we have:
λ2 − (2− σ2φ)λ+ 1− σ2φ+ σ4v + σ2α2 = 0. (D.7) The discriminant is ∆ := σ2(σ2(φ2 − 4ν)− 4α2). We discuss two cases:
1. φ2 − 4ν < 0. We are minimizing:
min α,u,v
√ 1 + (α2 − φ)σ21 + σ41ν ∨ √ 1 + (α2 − φ)σ2n + σ4nν,
with a ∨ b := max{a, b} a shorthand. A minimizer is at α → 0 and ν → φ2/4 (since φ2 < 4ν), where β1 = β2 = 2/(σ21 + σ 2 n) and α→ 0.
2. φ2 − 4ν ≥ 0. A lower bound is: min u |1− φσ21/2| ∨ |1− φσ2n/2|,
which is obtained iff 4α2 ∼ (φ2 − 4ν)t for all σ2. This is only possible if α → 0 and φ2 → 4ν, which yields β1 = β2 = 2/(σ21 + σ2n).
From what has been discussed, the optimal radius is (κ2 − 1)/(κ2 + 1) which can be achieved at β1 = β2 = 2/(σ 2 1 + σ 2 n) and α → 0. Again, this might not be the only choice. For instance, take σ1 = σ 2 n = 1, from equation 3.3, a super-linear convergence rate can be achieved at β1 = 1 and β2 = 1− α2.

heading: D.2 PROOF OF THEOREM 4.2: OPTIMAL CONVERGENCE RATE OF OGD
text: Theorem 4.2 (OGD optimal). For Jacobi OGD with β1 = β2 = β, to achieve the optimal linear rate, we must have α ≤ 2β. For the original OGD with α = 2β, the optimal linear rate r∗ satisfies
r2∗ = 1
2 +
1
4 √
2σ21
√ (σ21 − σ2n)(5σ21 − σ2n + √ (σ21 − σ2n)(9σ21 − σ2n)), (D.8)
at
β∗ = 1
4 √ 2
√ 3σ41 − (σ21 − σ2n)3/2 √ 9σ21 − σ2n + 6σ21σ2n − σ4n σ41σ 2 n . (D.9)
If κ → ∞, r∗ ∼ 1 − 1/(6κ2). For Gauss–Seidel OGD with β2 = 0, the optimal linear rate is r∗ = √ (κ2 − 1)/(κ2 + 1), at α = √ 2/σ1 and β1 = √ 2σ1/(σ 2 1 + σ 2 n). If κ→∞, r∗ ∼ 1− 1/κ2.
For OGD, the characteristic polynomials equation 3.6 and equation 3.7 are quartic and cubic separately, and thus optimizing the spectral radii for generalized OGD is difficult. However, we can study two special cases: for Jacobi OGD, we take β1 = β2; for Gauss–Seidel OGD, we take β2 = 0. In both cases, the spectral radius functions can be obtained by solving quadratic polynomials.

heading: D.2.1 JACOBI OGD
text: We assume β1 = β2 = β in this subsection. The characteristic polynomial for Jacobi OGD equation 3.6 can be written as: λ2(λ− 1)2 + (λα− β)2σ2 = 0. (D.10) Factorizing it gives two equations which are conjugate to each other:
λ(λ− 1)± i(λα− β)σ = 0. (D.11) The roots of one equation are the conjugates of the other equation. WLOG, we solve λ(λ − 1) + i(λα− β)σ = 0 which gives (1/2)(u± v), where
u = 1− iασ, v = √ 1− α2σ2 − 2i(α− 2β)σ. (D.12)
Denote ∆1 = 1− α2σ2 and ∆2 = 2(α− 2β)σ. If α ≥ 2β, v can be expressed as:
v = 1√ 2
(√√ ∆21 + ∆ 2 2 + ∆1 − i √√ ∆21 + ∆ 2 2 −∆1 ) =:
1√ 2 (a− ib), (D.13)
therefore, the spectral radius r(α, β, σ) satisfies:
r(α, β, σ)2 = 1
4
( (1 + a/ √ 2)2 + (ασ + b/ √ 2)2 ) = 1
4 (1 + α2σ2 +
√ ∆21 + ∆ 2 2 + √ 2(bσα+ a)),
(D.14) and the minimum is achieved at α = 2β. From now on, we assume α ≤ 2β, and thus v = a + ib. We write: r(α, β, σ)2 = 1
4 max{
( (1 + a/ √ 2)2 + (ασ − b/ √ 2)2 ) , ( (1− a/ √ 2)2 + (ασ + b/ √ 2)2 ) },
= 1
4 (1 + α2σ2 +
√ ∆21 + ∆ 2 2 + √ 2|bσα− a|).
=
{ 1 4 (1 + α 2σ2 + √ ∆21 + ∆ 2 2 − √
2(bσα− a)) 0 < ασ ≤ 1, 1 4 (1 + α 2σ2 + √ ∆21 + ∆ 2 2 + √ 2(bσα− a)) ασ > 1. (D.15)
This is a non-convex and non-differentiable function, which is extremely difficult to optimize. At α = 2β, in this case, a = √ 1− 4β2σ2sign(1− 4β2σ2) and b = √
4β2σ2 − 1sign(4β2σ2 − 1). The sign function sign(x) is defined to be 1 if x > 0 and 0 otherwise. The function we are optimizing is a quasi-convex function:
r(β, σ)2 =
{ 1 2 (1 + √ 1− 4β2σ2) 4β2σ2 ≤ 1,
2β2σ2 + βσ √ 4β2σ2 − 1 4β2σ2 > 1. (D.16)
We are maximizing over σ and minimizing over β. There are three cases:
• 4β2σ21 ≤ 1. At 4β2σ21 = 1, the optimal radius is:
r2∗ = 1
2
( 1 + √ 1− 1
κ2
) .
• 4β2σ2n ≥ 1. At 4β2σ2n = 1, the optimal radius satisfies:
r2∗ = κ2 2 + κ 2
√ κ2 − 1.
• 4β2σ2n ≤ 1 and 4β2σ21 ≥ 1. The optimal β is achieved at:
1
2
( 1 + √ 1− 4β2σ2n ) = 2β2σ21 + βσ1 √ 4β2σ21 − 1.
The solution is unique since the left is decreasing and the right is increasing. The optimal β is:
β∗ = 1
4 √ 2
√ 3σ41 − (σ21 − σ2n)3/2 √ 9σ21 − σ2n + 6σ21σ2n − σ4n σ41σ 2 n . (D.17)
The optimal radius satisfies:
r2∗ = 1
2 +
1
4 √
2σ21
√ (σ21 − σ2n)(5σ21 − σ2n + √ (σ21 − σ2n)(9σ21 − σ2n)). (D.18)
This is the optimal solution among the three cases. If σ2n/σ 2 1 is small enough we have r2 ∼ 1− 1/(3κ2).

heading: D.2.2 GAUSS–SEIDEL OGD
text: In this subsection, we study Gauss–Seidel OGD and fix β2 = 0. The characteristic polynomial equation 3.7 now reduces to a quadratic polynomial:
λ2 + (α2σ2 − 2)λ+ 1− αβ1σ2 = 0.
For convenience, we reparametrize β1 → β/α. So, the quadratic polynomial becomes:
λ2 + (α2σ2 − 2)λ+ 1− βσ2 = 0.
We are doing a min-max optimization minα,β maxσ r(α, β, σ), where r(α, β, σ) is:
r(α, β, σ) =
{√ 1− βσ2 α4σ2 < 4(α2 − β)
1 2 |α 2σ2 − 2|+ 12 √ α4σ4 − 4(α2 − β)σ2 α4σ2 ≥ 4(α2 − β).
(D.19)
There are three cases to consider:
• α4σ21 ≤ 4(α2 − β). We are minimizing 1− βσ2n over α and β. Optimizing over β1 gives β = α2−α4σ21/4. Then we minimize over α and obtain α2 = 2/σ21 . The optimal β = 1/σ21 and the optimal radius is √ 1− 1/κ2.
• α4σ2n > 4(α2 − β). Fixing α, the optimal β = α2 − α4σ2n/4, and we are solving
min α max
{ 1
2 |α2σ21 − 2|+
1 2 α2 √ σ21(σ 2 1 − σ2n), 1 2 |α2σ2n − 2| } .
We need to discuss three cases: α2σ2n > 2, α 2σ21 < 2 and 2/σ 2 1 < α 2 < 2/σ2n. In the first case, the optimal radius is κ2 − 1 + κ √
(κ2 − 1). In the second case, α2 → 2/σ21 and the optimal radius is √ 1− 1/κ2. In the third case, the
optimal radius is also √
1− 1/κ2 minimized at α2 → 2/σ21 . • α4σ21 > 4(α2 − β) and α4σ2n < 4(α2 − β). In this case, we have α2σ21 < 4. Otherwise, r(α, β, σ1) > 1. We are minimizing over:
max{ √ 1− βσ2n, 1
2 |α2σ21 − 2|+
1
2
√ α4σ41 − 4α2σ21 + 4βσ21}.
The minimum over α is achieved at α2σ21 = 2, and β = 2/(σ 2 1 +σ 2 n), this gives α = √ 2/σ1
and β1 = √ 2σ1/(σ 2 1 + σ 2 n). The optimal radius is r∗ = √ (κ2 − 1)/(κ2 + 1).
Out of the three cases, the optimal radius is obtained in the third case, where r ∼ 1− 1/κ2. This is better than Jacobi OGD, but still worse than the optimal EG.

heading: D.3 NUMERICAL METHOD
text: We first prove Lemma 4.1: Lemma 4.1. A polynomial p(λ) is r-Schur stable iff p(rλ) is Schur stable. Proof. Denote p(λ) = ∏n i=1(λ− λi). We have p(rλ) ∝ ∏n i=1(λ− λi/r), and:
∀i ∈ [n], |λi| < r ⇐⇒ ∀i ∈ [n], |λi/r| < 1. (D.20)
With Lemma 4.1 and Corollary 2.1, we have the following corollary: Corollary D.1. A real quadratic polynomial λ2 + aλ+ b is r-Schur stable iff b < r2, |a| < r+ b/r; A real cubic polynomial λ3 + aλ2 + bλ + c is r-Schur stable iff |c| < r3, |ar2 + c| < r3 + br, br4 − acr2 < r6 − c2; A real quartic polynomial λ4 + aλ3 + bλ2 + cλ + d is r-Schur stable iff |cr5 − adr3| < r8 − d2, |ar2 + c| < br + d/r + r3, and
b < r2 + dr−2 + r2 (cr2 − ad)(ar2 − c)
(d− r4)2 .
Proof. In Corollary 2.1, rescale the coefficients according to Lemma 4.1.
We can use the corollaries above to find the regions where r-Schur stability is possible, i.e., a linear rate of exponent r. A simple algorithm might be to start from r0 = 1, find the region S0. Then recursively take rt+1 = srt and find the Schur stable region St+1 inside St. If the region is empty then stop the search and return St. s can be taken to be, say, 0.99. Formally, this algorithm can be described as follows in Algorithm 1:
r0 = 1, t = 0, s = 0.99; Find the r0-Schur region S0; while St is not empty do
rt+1 = srt; Find the rt+1-Schur region St+1; t = t+ 1;
end Algorithm 1: Numerical method for finding the optimal convergence rate
In this algorithm, Corollary D.1 can be applied to obtain any r-Schur region.

heading: E SUPPLEMENTARY MATERIAL FOR SECTIONS 5 AND 6
text: We provide supplementary material for Sections 5 and 6. We first prove that when learning the mean of a Gaussian, WGAN is locally a bilinear game in Appendix E.1. For mixtures of Gaussians, we provide supplementary experiments about Adam in Appendix E.2. This result implies that in some cases, Jacobi updates are better than GS updates. We further verify this claim in Appendix E.3 by showing an example of OGD on bilinear games. Optimizing the spectral radius given a certain singular value is possible numerically, as in Appendix E.4.

heading: E.1 WASSERSTEIN GAN
text: Inspired by Daskalakis et al. (2018), we consider the following WGAN (Arjovsky et al., 2017):
f(φ,θ) = min φ max θ Ex∼N (v,σ2I)[s(θ
>x)]− Ez∼N (0,σ2I)[s(θ>(z + φ))], (E.1)
with s(x) := 1/(1 + e−x) the sigmoid function. We study the local behavior near the saddle point (v,0), which depends on the Hessian:[
∇2φφ ∇2φθ ∇2θφ ∇2θθ
] = [ −Eφ[s′′(θ>z)θθ>] −Eφ[s′′(θ>z)θz> + s′(θ>z)I]
(∇2φθ)> Ev[s′′(θ>x)xx>]− Eφ[s′′(θ>z)zz>]
] ,
with Ev a shorthand for Ex∼N (v,σ2I) and Eφ for Ez∼N (φ,σ2I). At the saddle point, the Hessian is simplified as:
[ ∇2φφ ∇2φθ ∇2θφ ∇2θθ ] = [ 0 −s′(0)I −s′(0)I 0 ] = [ 0 −I/4 −I/4 0 ] .
Therefore, this WGAN is locally a bilinear game.

heading: E.2 MIXTURES OF GAUSSIANS WITH ADAM
text: Given the same parameter settings as in Section 5, we train the vanilla GAN using Adam, with the step size α = 0.0002, and β1 = 0.9, β2 = 0.999. As shown in Figure 6, Jacobi updates converge faster than the corresponding GS updates.

heading: E.3 JACOBI UPDATES MAY CONVERGE FASTER THAN GS UPDATES
text: Take α = 0.9625, β1 = β2 = β = 0.5722, and σ = 1, the Jacobi and GS OGD radii are separately 0.790283 and 0.816572 (by solving equation 3.6 and equation 3.7), which means that Jacobi OGD has better performance for this setting of parameters. A more intuitive picture is given as Figure 7, where we take β1 = β2 = β.

heading: E.4 SINGLE SINGULAR VALUE
text: We minimize r(θ, σ) for a given singular value numerically. WLOG, we take σ = 1, since we can rescale parameters to obtain other values of σ. We implement grid search for all the parameters within the range [−2, 2] and step size 0.05. For the step size α, we take it to be positive. We use {a, b, s} as a shorthand for {a, a+ s, a+ 2s, . . . , b}.
• We first numerically solve the characteristic polynomial for Jacobi OGD equation 3.6, fixing α1 = α2 = α with scaling symmetry. With α ∈ {0, 2, 0.05}, βi ∈ {−2, 2, 0.05}, the best parameter setting is α = 0.7, β1 = 0.1 and β2 = 0.6. β1 and β2 can be switched. The optimal radius is 0.6.
• We also numerically solve the characteristic polynomial for Gauss–Seidel OGD equation 3.7, fixing α1 = α2 = α with scaling symmetry. With α ∈ {0, 2, 0.05}, βi ∈ {−2, 2, 0.05}, the best parameter setting is α = 1.4, β1 = 0.7 and β2 = 0. β1 and β2 can be switched. The optimal rate is 1/(5 √ 2). This rate can be further improved to be zero where α = √ 2,
β1 = 1/ √ 2 and β2 = 0.
• Finally, we numerically solve the polynomial for Gauss–Seidel momentum equation 3.11, with the same grid. The optimal parameter choice is α = 1.8, β1 = −0.1 and β2 = −0.05. β1 and β2 can be switched. The optimal rate is 0.5.

heading: F SPLITTING METHOD
text: In this appendix, we interpret the gradient-based algorithms (except PP) we have studied in this paper as splitting methods (Saad, 2003), for both Jacobi and Gauss–Seidel updates. By doing this, one can understand our algorithms better in the context of numerical linear algebra and compare our results in Section 3 with the Stein–Rosenberg theorem.

heading: F.1 JACOBI UPDATES
text: From equation 2.2, finding a saddle point is equivalent to solving:
Sz := [ 0 E −E> 0 ] [ x y ] = [ −b c ] =: d. (F.1)
Now, we try to understand the Jacobi algorithms using splitting method. For GD and EG, the method splits S intoM −N and solve
zt+1 = M −1Nzt +M −1d. (F.2) For GD, we can obtain that:
M =
[ α−11 I 0
0 α−12 I
] , N = [ α−11 I −E E> α−12 I ] . (F.3)
For EG, we need to compute an inverse:
M−1 = [ α1I −β1E β2E > α2I ] , N = M − S. (F.4)
Given det(α1α2I + β1β2EE>) 6= 0, the inverse always exists. The splitting method can also work for second-step methods, such as OGD and momentum. We split S = M −N − P and solve:
zt+1 = M −1Nzt +M −1Pzt−1 +M −1d. (F.5)
For OGD, we have:
M =
[ I
α1−β1 0
0 Iα2−β2
] , N =  Iα1−β1 − α1Eα1−β1 α2E >
α2−β2 I α2−β2
 , P = [ 0 β1Eα1−β1− β2E>α2−β2 0 ] . (F.6)
For the momentum method, we can write:
M =
[ α−11 I 0
0 α−12 I
] , N = [ 1+β1 α1 I −E
E> 1+β2α2 I
] , P = [ − β1α1 I 0
0 − β2α2 I
] . (F.7)

heading: F.2 GAUSS–SEIDEL UPDATES
text: Now, we try to understand the GS algorithms using splitting method. For GD and EG, the method splits S intoM −N and solve
zt+1 = M −1Nzt +M −1d. (F.8) For GD, we can obtain that:
M = [ α−11 I 0 −E> α−12 I ] , N = [ α−11 I −E 0 α−12 I ] . (F.9)
For EG, we need to compute an inverse:
M−1 =
[ α1I −β1E
(β2 + α1α2)E > α2(I − β1E>E)
] , N = M − S. (F.10)
The splitting method can also work for second-step methods, such as OGD and momentum. We split S = M −N − P and solve:
zt+1 = M −1Nzt +M −1Pzt−1 +M −1d. (F.11)
For OGD, we obtain:
M =  Iα1−β1 0 − α2E >
α2−β2 I α2−β2
 , N =  I α1−β1 − α1E α1−β1
− β2E >
α2−β2 I α2−β2
 , P = 0 β1Eα1−β1
0 0  . (F.12) For the momentum method, we can write:
M = [ α−11 I 0 −E> α−12 I ] , N = [ 1+β1 α1 I −E 0 1+β2α2 I ] , P = [ − β1α1 I 0 0 − β2α2 I ] . (F.13)

heading: G SINGULAR BILINEAR GAMES
text: In this paper we considered the bilinear game when E is a non-singular square matrix for simplicity. Now let us study the general case where E ∈ Rm×n. As stated in Section 2, saddle points exist iff
b ∈ R(E), c ∈ R(E>). (G.1)
Assume b = Eb′, c = E>c′. One can shift the origin of x and y: x→ x− b′, y → y − c′, such that the linear terms cancel out. Therefore, the min-max optimization problem becomes:
min x∈Rm max y∈Rn
x>Ey. (G.2)
The set of saddle points is:
{(x,y)|y ∈ N (E),x ∈ N (E>)}. (G.3)
For all the first-order algorithms we study in this paper, x(t) ∈ x(0)+R(E) and y(t) ∈ y(0)+R(E>). Since for any matrix X ∈ Rp×q, R(X) ⊕ N (X>) = Rp, if the algorithm converges to a saddle point, then this saddle point is uniquely defined by the initialization:
x∗ = P⊥Ex (0), y∗ = P⊥E>y (0), (G.4)
where
P⊥X := I −X†X, (G.5)
is the orthogonal projection operator onto the null space ofX , andX† denotes the Moore–Penrose pseudoinverse. Therefore, the convergence to the saddle point is described by the distances of x(t) and y(t) to the null spaces N (E>) and N (E). We consider the following measure:
∆2t = ||E†Ey(t)||2 + ||EE†x(t)||2, (G.6)
as the Euclidean distance of z(t) = (x(t),y(t)) to the space of saddle points N (E>) × N (E). Consider the singular value decomposition of E:
E = U [ Σr 0 0 0 ] V >, (G.7)
with Σr ∈ Rr×r diagonal and non-singular. Define:
v(t) = V >y(t), u(t) = U>x(t), (G.8)
and equation G.6 becomes:
∆2t = ||v(t)r ||2 + ||u(t)r ||2, (G.9)
with vr denoting the sub-vector with the first r elements of v. Hence, the convergence of the bilinear game with a singular matrix E reduces to the convergence of the bilinear game with a non-singular matrix Σr, and all our previous analysis still holds.

Review summaries:
This paper studies the convergence of multiple methods (Gradient, extragradient, optimistic and momentum) on a bilinear minmax game. More precisely, this paper uses spectral condition to study the difference between simultaneous (Jacobi) and alternating (Gau\ss-Seidel) updates. The analysis is based on Schur theorem and give necessary and sufficient condition for convergence. 

The paper presents exact conditions for the convergence of several gradient based methods for solving bilinear games. In particular, the methods under study are Gradient Descent(GD), Extragradient (EG), Optimizatic Gradient descent (OGD) and Momentum methods. For these methods, the authors provide convergence rates (with optimal parameter setup) for both alternating (Gauss-Seidel) and simultaneous (Jacobi) updates.