GDA-AM: ON THE EFFECTIVENESS OF SOLVING MIN- IMAX OPTIMIZATION VIA
Year: 2021

Abstract:
Many modern machine learning algorithms such as generative adversarial networks (GANs) and adversarial training can be formulated as minimax optimization. Gradient descent ascent (GDA) is the most commonly used algorithm due to its simplicity. However, GDA can converge to non-optimal minimax points. We propose a new minimax optimization framework, GDA-AM, that views the GDA dynamics as a fixed-point iteration and solves it using Anderson Mixing to converge to the local minimax. It addresses the diverging issue of simultaneous GDA and accelerates the convergence of alternating GDA. We show theoretically that the algorithm can achieve global convergence for bilinear problems under mild conditions. We also empirically show that GDA-AM solves a variety of minimax problems and improves GAN training on several datasets.

heading: 1 INTRODUCTION
text: Minimax optimization has received a surge of interest due to its wide range of applications in modern machine learning, such as generative adversarial networks (GAN), adversarial training and multiagent reinforcement learning (Goodfellow et al., 2014; Madry et al., 2018; Li et al., 2019). Formally, given a bivariate function f(x,y), the objective is to find a stable solution where the players cannot improve their objective, i.e., to find the Nash equilibrium of the underlying game (von Neumann & Morgenstern, 1944):
arg min x∈X arg max y∈Y f(x,y). (1)
It is commonplace to use simple algorithms such as gradient descent ascent (GDA) to solve such problems, where both players take a gradient update simultaneously or alternatively. Despite its simplicity, GDA is known to suffer from a generic issue for minimax optimization: it may cycle around a stable point, exhibit divergent behavior, or converge very slowly since it requires very small learning rates (Gidel et al., 2019a; Mertikopoulos et al., 2019). Given the widespread usage of gradient-based methods for solving machine learning problems, first-order optimization algorithms to solve minimax problems have gained considerable popularity in the last few years. Algorithms such as optimistic Gradient Descent Ascent (OG) (Daskalakis et al., 2018; Mertikopoulos et al., 2019) and extra-gradient (EG) (Gidel et al., 2019a) can alleviate the issue of GDA for some problems. Yet, it has been shown that these methods can still diverge or cycle around a stable point (Adolphs et al.; Mazumdar et al., 2019; Parker-Holder et al., 2020). For example, these algorithms even fail to find a local minimax (the set of local minimax is a superset of local Nash (Jin et al., 2020; Wang et al., 2020)) as shown in Figure 1. This leads to the following question: Can we design better algorithms for minimax problems? We answer this in the affirmative, by introducing GDA-AM. We cast the GDA dynamics as a fixed-point iteration problem and compute the iterates effectively using an advanced nonlinear extrapolation method. We show that indeed our algorithm has theoretical and empirical guarantees across a broad range of minimax problems, including GANs.
Our contributions: In this paper, we propose a different approach to solve minimax optimization. Our starting point is to cast the GDA dynamics as a fixed-point iteration. We then highlight that the fixed-point iteration can be solved effectively by using advanced non-linear extrapolation methods such as Anderson Mixing (Anderson, 1965), which we name as GDA-AM. Although first mentioned in Azizian et al. (2020), to our best knowledge, this is still the first work to investigate and improve the GDA dynamics by tapping into advanced fixed-point algorithms.
We demonstrate that GDA dynamics can benefit from Anderson Mixing. In particular, we study bilinear games and give a systematic analysis of GDA-AM for both simultaneous and alternating versions of GDA. We theoretically show that GDA-AM can achieve global convergence guarantees under mild conditions.
We complement our theoretical results with numerical simulations across a variety of minimax problems. We show that for some convex-concave and non-convex-concave functions, GDA-AM can converge to the optimal point with little hyper-parameter tuning whereas existing first-order methods are prone to divergence and cycling behaviors.
We also provide empirical results for GAN training across two different datasets, CIFAR10 and CelebA. Given the limited computational overhead of our method, the results suggest that an extrapolation add-on to GDA can lead to significant performance gains. Moreover, the convergence behavior across a variety of problems and the ease-of-use demonstrate the potential of GDA-AM to become the minimax optimization workhorse.

heading: 2 PRELIMINARIES AND BACKGROUND

heading: 2.1 MINIMAX OPTIMIZATION
text: Definition 1. Point (x∗,y∗) is a local Nash equilibrium of f if there exists δ > 0 such that for any (x,y) satisfying ‖x− x∗‖ ≤ δ and ‖y − y∗‖ ≤ δ we have: f (x∗,y) ≤ f (x∗,y∗) ≤ f (x,y∗) .
To find the Nash equilibria, common algorithms including GDA, EG and OG, can be formulated as follows. For the two variants of GDA, simultaneous GDA (SimGDA) and alternating GDA (AltGDA), the updates have the following forms:
Simultaneous : xt+1 = xt − η∇xf(xt,yt), yt+1 = yt + η∇yf(xt,yt) Alternating : xt+1 = xt − η∇xf(xt,yt), yt+1 = yt + η∇yf(xt+1,yt).
(2)
The EG update has the following form:
xt+ 12 = xt − η∇xf(xt,yt), yt+ 12 = yt + η∇yf(xt,yt) xt+1 = xt − η∇xf(xt+ 12 ,yt+ 12 ), yt+1 = yt + η∇yf(xt+ 12 ,yt+ 12 ).
(3)
The OG update has the following form:
xt+1 = xt−η∇xf(xt,yt)+ η
2 ∇xf(xt−1,yt−1),yt+1 = yt+η∇yf(xt,yt)−
η 2 ∇yf(xt−1,yt−1).
(4)

heading: 2.2 FIXED-POINT ITERATION AND ANDERSON MIXING (AM)
text: Definition 2. w? is a fixed point of the mapping g if w? = g (w?) .
Consider the simple fixed-point iteration wt+1 = g(wt) which produces a sequence of iterates {w0, w1, · · · , wN}. In most cases, this converges to the fixed-point, w∗ = g(w∗). Take gradient descent as an example, it can be viewed as iteratively applying the operation: wt+1 = g (wt) , wt−αt∇f (wt) ,where the limit is the fixed-point w? = g (w?) (i.e.∇f (wt = 0) . SimGDA updates can be defined as the repeated application of a nonlinear operator:
wt+1 = G (sim) η (wt) , wt − ηV (wt) with w = [ x y ] , V (w) = [ ∇xf(x,y) −∇yf(x,y) ] Similarly, we can write AltGDA updates as wt+1 = G (alt) η (wt). An issue with fixed-point iteration is that it does not always converge, and even in the cases where it does converge, it might do so very slowly. GDA is one example that it could result in the possibility of the operator converging to a limit cycle instead of a single point for the GDA dynamic. A way of dealing with these problems is to use acceleration methods, which can potentially speed up the convergence process and in some cases even decrease the likelihood for divergence.
There are many different acceleration methods, but we will put our focus on an algorithm which we refer to as Anderson Mixing (or Anderson Acceleration). In short, Anderson Mixing (AM) shares the same idea as Nesterov’s acceleration. Given a fixed-point iteration wt = g (wt−1), Anderson Mixing argues that a good approximation to the final solution w∗ can be obtained as a linear combination of the previous p iterates wt+1 = ∑p i=0 βig (wt−pt+i). Since obtaining the proper coefficients βi is a nonlinear procedure, Anderson Mixing is also known as a nonlinear extrapolation method. The general form of Anderson Mixing is shown in Algorithm 1. For efficiency, we prefer a ‘restarted’ version with a small table size p that cleans up the table F every p iterations because it avoids solving a linear system of increasing size.
Algorithm 1: Anderson Mixing Prototype (truncated version) Input: Initial point w0, Anderson restart dimension p, fixed-point mapping g : Rn → Rn. Output: wt+1 for t = 0, 1, . . . do
Set pt = min{t, p}. Set Ft = [ft−pt , . . . , ft], where fi = g(wi)− wi for each i ∈ [t− pt . . t]. Determine weights β = (β0, . . . , βpt) T that solves minβ ‖Ftβ‖2 , s. t. ∑pt i=0 βi = 1.
Set wt+1 = ∑pt i=0 βig (wt−pt+i).
end

heading: 2.3 AM AND GENERALIZED MINIMAL RESIDUAL (GMRES)
text: Developed by Saad & Schultz (1986), Generalized Minimal Residual method (GMRES) is a Krylov subspace method for solving linear system equations. The method approximates the solution by the vector in a Krylov subspace with minimal residual, which is described below. Definition 3. Assume we have the linear system of equations Ax = b with A ∈ Rn×n,b ∈ Rn and an initial guess x0. Then we denote the initial residual by r0 = b−Ax0 and define the tth Krylov subspace as Kt = span{r0,Ar0, · · · ,At−1r0}.
The tth iterate xt of GMRES minimizes the norm of the residual rt = b −Axt in Kt, that is, xt solves
min xt∈x0+Kt ‖b−Axt‖2 .
The following formulation is equivalent to GMRES minimization problem and more convenient for implementation. It computes x̂t such that
x̂t = arg min x̂t∈Kt ‖b−A (x0 + x̂t)‖2 = arg min x̂t∈Kt ‖r0 −Ax̂t‖2 .
Using a larger Krylov dimension will improve the convergence of the method, but will require more memory. For this reason, a smaller Krylov subspace dimension t and ‘restarted’ versions of the method are used in practice Saad (2003).
The convergence of GMRES can be studied through the magnitude of the residual polynomial.
Theorem 2.1 (Lemma 6.31 of Saad (2003)). Let x̂t be the approximate solution obtained at the t-th iteration of GMRES being applied to solve Ax = b, and denote the residual as rt = b−Ax̂t. Then, rt is of the form
rt = ft(A)r0, (5)
where
‖rt‖2 = ‖ft(A)r0‖2 = min ft∈Pt ‖ft(A)r0‖2, (6)
where Pp is the family of polynomials with degree p such that fp(0) = 1,∀fp ∈ Pp, which are usually called residual polynomials.
Although GMRES is applied to a system of linear equations not a fixed-point problem, there is a strong connection between Anderson Mixing and GMRES. In AM we are looking for a fixed-point x such that Gx− b− x = 0 and by rearranging this equation we get
b + (G− I)x = 0⇔ (I−G)x = b.
Theorem 2.2 shows that if GMRES is applied to the system (I −G)x = b and AM is applied to g(x) = Gx + b with the same initial guess and I−G is non-singular, then these are equivalent in the sense that the iterates of each algorithm can be obtained directly from the iterates of the other algorithm.
Theorem 2.2 (Equivalence between AM with restart and GMRES (Walker & Ni, 2011a)). Consider the fixed point iteration x = g(x) where g(x) = Gx + b for G ∈ Rn×n and b ∈ Rn. If I −G is non-singular, Algorithm 1 produces exactly the same iterates as GMRES being applied to solve (I−G)x = b when both algorithms start with the same initial guess.
Theorem 2.2 can also be generalized to the restart version of AM an GMRES as well.
3 GDA-AM : GDA WITH ANDERSON MIXING
We propose a novel minimax optimizer, called GDA-AM, that is inspired by recent advances in parameter (or weight) averaging (Wu et al., 2020; Yazici et al., 2019). We argue that a nonlinear adaptive average (combination) is a more appropriate choice for minimax optimization.

heading: 3.1 GDA WITH NAÏVE ANDERSON MIXING
text: We propose to exploit the dynamic information present in the GDA iterates to “smartly” combine the past iterates. This is in contrast to the classical averaging methods (moving averaging and exponential moving averaging) (Yang et al., 2019) that “blindly” combine past iterates. A naı̈ve adoption of Anderson Mixing using the past p GDA iterates for both simGDA and altGDA has the following form:
Anderson mixing : xt+1 = p∑ i=0 βixt−p+i,yt+1 = p∑ i=0 βiyt−p+i.
Since Zhang et al. (2021); Gidel et al. (2019b) show the AltGDA is superior to SimGDA in many aspects, we briefly summarized both Simultaneous and Alternating GDA-AM in Algorithms 2 and 3 with the truncated Anderson Mixing Algorithm 1 using a table size p.
Algorithm 2: Simultaneous GDA-AM Input: x0,y0, stepsize η, Anderson table size p Output: xt,yt Set w0 = [x0,y0], sx = length(x0) for t = 0, 1, . . . do
xt,yt = wt[0 : sx− 1],wt[sx : end] xt+1 = xt − η∇xf(xt,yt) yt+1 = yt − η∇yf(xt,yt)
wt+1 = [ xt+1 yt+1 ] Use Anderson Mixing with table size p
to extrapolate wt+1 end xt,yt = wt+1[0 : sx− 1],wt+1[sx : end] return xt,yt
Algorithm 3: Alternating GDA-AM Input: x0,y0, stepsize η, Anderson table size p Output: xt,yt Set w0 = [x0,y0], sx = length(x0) for t = 0, 1, . . . do
xt,yt = wt[0 : sx− 1],wt[sx : end] xt+1 = xt − η∇xf(xt,yt) yt+1 = yt − η∇yf(xt+1,yt)
wt+1 = [ xt+1 yt+1 ] Use Anderson Mixing with table size p to extrapolate wt+1 end xt,yt = wt+1[0 : sx− 1],wt+1[sx : end] return xt,yt

heading: 3.2 IMPLEMENTATION DETAILS
text: It is important to note that the Anderson Mixing form shown in Algorithm 1 is for illustrative purpose and not computationally efficient. For example, only one column of Ft needs to be updated at each iteration. In addition, the solution of the least-square problem in Algorithm 1 can also be solved by a quick QR update scheme which costs (2n+ 1)p2 (Walker & Ni, 2011a). Thus, from Algorithms 2 and 3, we can see that the major cost of GDA-AM arises from solving the additional linear least squares problem compared to regular GDA at each iteration. Additional implementation details are provided in the Appendix.
4 CONVERGENCE RESULTS FOR GDA-AM
In this section, we show that both simultaneous and alternating version GDA-AM converge to the equilibrium for bilinear problems. First, we do not require the learning rate to be sufficiently small. Second, we explicitly provide a linear convergence rate that is faster than EG and OG.

heading: 4.1 BILINEAR GAMES
text: Bilinear games are often regarded as an important simple example for theoretically analyzing and understanding new algorithms and techniques for solving general minimax problems (Gidel et al., 2019a; Mertikopoulos et al., 2019; Schaefer & Anandkumar, 2019). In this section, we analyze the convergence property of simultaneous GDA-AM and alternating GDA-AM schemes on the following zero-sum bilinear games:
min x∈Rn max y∈Rn
f(x,y) = xTAy + bTx + cTy, A is full rank. (7)
The Nash equilibrium to the above problem is given by (x∗,y∗) = (−A−T c,−A−1b). We also investigate bilinear-quadratic games from a spectrum analysis perspective. In addition, we show that analysis based on the numerical range (Bollapragada et al., 2018) can be also extended to such games, although it can not help derive a convergent bound for equation 7. Detailed discussion can be found in Appendix B.1 and B.4.1.
4.2 SIMULTANEOUS GDA-AM
Suppose x0 and y0 are the initial guesses for x∗ and y∗, respectively. Then each iteration of simultaneous GDA can be written in the following matrix form:[
xt+1 yt+1
] = [ I −ηA
ηAT I ] ︸ ︷︷ ︸
G(Sim)
[ xt yt ] ︸︷︷︸ w
(Sim) t
−η [ b c ] ︸︷︷︸ b(Sim) . (8)
It has been shown that the iteration in equation 8 often cycles and fails to converge for the bilinear problem due to the poor spectrum/numerical range of the fixed point operator G(Sim) (Gidel et al., 2019a; Azizian et al., 2020; Mokhtari et al., 2020a). Next we show that the convergence can be improved with Algorithm 2. Theorem 4.1. [Global convergence for simultaneous GDA-AM on bilinear problem] Denote the distance between the stationary point w∗ and current iterate w(k+1)p of Algorithm 2 with table size p as N(k+1)p = ‖w∗ −w(k+1)p‖. Then we have the following bound for Nt
N2(k+1)p ≤ ρ(A)N 2 kp (9)
where ρ(A) = ( 1 Tp(1+ 2 κ(ATA)−1 ) )2. Here, Tp is the Chebyshev polynomial of first kind of degree p and 1 Tp(1+ 2 κ(ATA)−1 ) < 1 since 1 + 2 κ(ATA)−1 > 1.
It is worthy emphasizing that the convergence rate of Algorithm 2 is independent of learning rate η while the convergence results of other methods like EG and OG depend on the learning rate. Remark 4.1.1. Both EG and OG have the following form of convergence rate (Mokhtari et al., 2020a) for bilinear problem
N2t+1 ≤ (1− c
κ(ATA) )N2t ,
where c is a positive constant independent of the problem parameters.
4.3 ALTERNATING GDA-AM
The underlying fixed point iteration in Algorithm 3 can be written in the following matrix form:[ xt+1 yt+1 ] = [ I −ηA ηAT I− η2ATA ] ︸ ︷︷ ︸
G(Alt)
[ xt yt ] ︸︷︷︸ w
(Alt) t
−η [ b c ] ︸︷︷︸ b(Alt) .
According to the equivalence between truncated Anderson acceleration and GMRES with restart, we can analyze the convergence of Algorithm 3 through the convergence analysis of applying GMRES to solve linear systems associated with G = I−G(Alt):
G =
[ 0 ηA
−ηAT η2ATA
] .
Theorem 4.2. [Global convergence for alternating GDA-AM on bilinear problem] Denote the distance between the stationary point w∗ and current iterate w(k+1)p of Algorithm 3 with table size p as N(k+1)p = ‖w∗ − w(k+1)p‖. Assume A is normalized such that its largest singular value is equal to 1. Then when the learning rate η is less than 2, we have the following bound for Nt
N2(k+1)p ≤ √ 1 + 2η
2− η ( r c )pN2kp
where c and r are the center and radius of a diskD(c, r) which includes all the eigenvalues of G. Especially, rc < 1.
Figure 3: An illustration of the spectrum of G (red) and the closing circle (blue) in Theorem 4.2.
Theorem 4.2 shows that when p > log
√ 2−η 2+η
log rc , alternating GDA-AM will converge globally.

heading: 5 RELATED WORK
text: There is a rich literature on different strategies to alleviate the issue of minimax optimization. A useful add-on technique, Momentum, has been shown to be effective for bilinear games and strongly-convexstrongly-concave settings (Zhang & Wang, 2021; Gidel et al., 2019b; Azizian et al., 2020). Several second-order methods (Adolphs et al.; Mescheder et al., 2017; Mazumdar et al., 2019; Parker-Holder et al., 2020) show that their stable fixed points are exactly either Nash equilibria or local minimax by incorporating second-order information. However, such methods are computationally expensive and thus unsuitable for large applications such as image generation. Focusing on variants of GDA, EG and OG are two widely studied algorithms on improving the GDA dynamics. EG proposed to apply extra-gradient to overcome the cycling behaviour of GDA. OG, originally proposed in Popov (1980) and rediscovered in Daskalakis et al. (2018); Mertikopoulos et al. (2019), is more efficient by storing and re-using the extrapolated gradient for the extrapolation step. Without projection, OG is equivalent to extrapolation from past. Mokhtari et al. (2020b) shows that both of these algorithms can be interpreted as approximations of the classical proximal point method and did a unified analysis for bilinear games. These approaches mentioned the GDA dynamics can be viewed as a fixed-point iteration, but none of them further provides a solution to improve it. In this work, we fill this gap by proposing the application of the extrapolation method directly on the entire GDA dynamics. Unlike OG, EG and their variants (Hsieh et al., 2019; Lei et al., 2021; Thekumparampil et al., 2019; Yang et al., 2019), which regard minimax problems as variational inequality problems (Bruck, 1977; Nemirovski, 2004), our work is from a new perspective and thus orthogonal to these previous approaches.
In addition, several recent works consider nonconvex-concave minimax problems. Zhang et al. (2020) introduced a “smoothing” scheme combined with GDA to stabilize the dynamic of GDA. Luo et al. (2020) proposed a method called Stochastic Recursive gradiEnt Descent Ascent (SREDA) for stochastic nonconvex-strongly-concave minimax problems, by estimating gradients recursively and reducing its variance. Lin et al. (2020) showed that the two-timescale GDA can find a stationary point of nonconvex-concave minimax problems effectively. Ostrovskii et al. (2021) proposed a variant of Nesterov’s accelerated algorithm to find -first-order Nash equilibrium that is a stronger criterion than the commonly used proximal gradient norm. Nouiehed et al. (2019) proposed a iterative method that finds -first-order Nash equilibrium in O( −2) iterations under Polyak-Lojasiewicz (PL) condition. Focusing on nonconvex minimax problems, they studied an interesting and difficult problem. Since our work cast insight on the effectiveness of solving minimax optimization via Anderson Mixing, we expect the extension of this algorithm to general nonconvex problems can be further investigated in the future.

heading: 6 EXPERIMENTS
text: In this section, we conduct experiments to see whether GDA-AM improves GDA for minimax optimization from simple to practical problems. We first investigate performance of GDA-AM on bilinear games. In addition, we evaluate the efficacy of our approach on GANs.

heading: 6.1 BILINEAR PROBLEMS
text: In this section, we answer following questions: Q1: How is GDA-AM perform in terms of iteration number and running time? Q2: How is the scalability of GDA-AM ? Q3: How is the performance of GDA-AM using different table size p? Q4: Does GDA-AM converge for large step size η?
We compare the performance with SimGDA, AltGDA, EG, and OG, and EG with Negative Momentum(Azizian et al. (2020)) on bilinear minimax games shown in equation 7 without any constraint.
A,b, c, and initial points are generated using normally distributed random number. We set the maximum iteration number as 1× 106, stopping criteria 1× 10−5 and depict convergence by use of the norm of distance to optima, which is defined as ‖w∗ −wt‖. Similar to Azizian et al. (2020); Wei et al. (2021a), the step size is set as 1 after rescaling A to have 2-norm 1. We present results of different settings in Figures 4, 5, and 6.
We first generate different problem size (n = 100, 1000, 5000) and present results of convergence in terms of iteration number in Figure 4. It can be observed that GDA-AM converges in much fewer iterations for different problem sizes. Note that EG, EG-NM, and OG converge in the end but requires many iterations, thus we plot only a portion for illustrative purposes. Figure 5 depicts the convergence for all methods in terms of time. It can be observed that the running time of GDA-AM is faster than EG. Although slower than OG, we can observe GDA-AM converges in much less time for all problems. Figure 4 and Figure 5 answer Q1 and Q2; although there is additional computation for GDA-AM , it does not hinder the benefits of adopting Anderson Mixing. Even for a large problem size, GDA-AM still converges in much less time than the baselines.
Next, we run GDA-AM using different table size p and show the results in Figure 6a and Figure 6b. Figure 6a indicates an increasing of table size results in faster convergence in terms of iteration number, which also verifies our claim in Theorem 4.1. However, we also observe an increased running time when using a larger table size in Figure 6b. Further, we can see that p = 50 converges in a comparable time and iterations to p = 100. Similar results are found in repeated experiments as well. As a result, our answer to Q3 is that although a larger p means less iterations, a medium p is sufficient and a small p still outperforms the baselines. The optimal choice of p is related to the condition number and step size, which is another interesting topic in the Anderson Mixing community.
Next, we answer Q4 on convergence under different step sizes. Although GDA-AM usually converges with suitable step size, our theorem suggests it requires a larger table size when combined with a extremely aggressive step size. Figure 6c shows the convergence under such circumstance. We can observe that although a very large step size goes the wrong way in the beginning, Anderson Mixing can still make it back on track except when η > 1. It answers the question and confirms our claim that GDA-AM can achieve global convergence for bilinear problems for a large step size η > 0.

heading: 6.2 GAN EXPERIMENTS: IMAGE GENERATION
text: We apply our method to the CIFAR10 dataset (Krizhevsky, 2009) and use the ResNet architecture with WGAN-GP (Gulrajani et al., 2017) and SNGAN (Miyato et al., 2018) objective. We also compared
the performance of GDA-AM using cropped CelebA (64×64) (Liu et al., 2015) on WGAN-GP. We compare with Adam and extra-gradient with Adam (EG) as it offers significant improvement over OG. Models are evaluated using the inception score (IS) (Salimans et al., 2016) and FID (Heusel et al., 2017) computed on 50,000 samples. For fair comparison, we fixed the same hyperparamters of Adam for all methods after an extensive search. Experiments were run with 5 random seeds. We show results in Table 1. Table 1 reports the best IS and FID (averaged over 5 runs) achieved on these datasets by each method. We see that GDA-AM yields improvements over the baselines in terms of generation quality.

heading: 7 CONCLUSION
text: We prove the convergence property of GDA-AM and obtain a faster convergence rate than EG and OG on the bilinear problem. Empirically, we verify our claim for such a problem and show the efficacy of GDA-AM in a deep learning setting as well. We believe our work is different from previous approaches and takes an important step towards understanding and improving minimax optimization by exploiting the GDA dynamic and reforming it with numerical techniques.

heading: A ANDERSON MIXING IMPLEMENTATION DETAILS
text: In this section, we discuss the efficient implementation of Anderson Mixing. We start with generic Anderson Mixing prototype (Algorithm 4) and then present the idea of Quick QR-update Anderson Mixing implementation as described in Walker & Ni (2011b), which is commonly used in practice. For each iteration t ≥ 0 , AM prototype solves a least squares problem with a normalization constraint. The intuition is to minimize the norm of the weighted residuals of the previous m iterates.
Algorithm 4: Anderson Mixing Prototype (truncated version) Input: Initial point w0, Anderson restart dimension p, fixed-point mapping g : Rn → Rn. Output: wt+1 for t = 0, 1, . . . do
Set pt = min{t, p}. Set Ft = [ft−pt , . . . , ft], where fi = g(wi)− wi for each i ∈ [t− pt . . t]. Determine weights β = (β0, . . . , βpt) T that solves minβ ‖Ftβ‖2 , s. t. ∑pt i=0 βi = 1.
Set wt+1 = ∑pt i=0 βig (wt−pt+i).
end
The constrained linear least-squares problem in Algorithm AA can be solved in a number of ways. Our preference is to recast it in an unconstrained form suggested in Fang & Saad (2009); Walker & Ni (2011b) that is straightforward to solve and convenient for implementing efficient updating of QR.
Define fi = g(wi) − wi, 4fi = fi+1 − fi for each i and set Ft = [ft−pt , . . . , ft], Ft = [4ft−pt , . . . ,4ft]. Then solving the least-squares problem (minβ ‖Ftβ‖2 , s. t. ∑pt i=0 βi = 1) is equivalent to min
γ=(γ0,...,γpt−1) T ‖ft −Ftγ‖2 (10)
where α and γ are related by α0 = γ0, αi = γi − γi−1 for 1 ≤ i ≤ pt − 1, and αpt = 1− γpt−1. Now the inner minimization subproblem can be efficiently solved as an unconstrained least squares problem by a simple variable elimination. This unconstrained least-squares problem leads to a modified form of Anderson Mixing
wt+1 = g (wt)− pt−1∑ i=0 γ (t) i [g (wt−pt+i+1)− g (wt−pt+i)] = g (wt)− Gtγ (t)
where Gt = [4gt−pt , . . . ,4gt−1] with4gi = g(wi+1)− g(wi) for each i. To obtain γ(t) = ( γ
(t) 0 , . . . , γ (t) pt−1
)T by solving equation 10 efficiently, we show how the successive
least-squares problems can be solved efficiently by updating the factors in the QR decomposition Ft = QtRt as the algorithm proceeds. We assume a think QR decomposition, for which the solution of the least-squares problem is obtained by solving the pt × pt linear system Rγ = Q′ ∗ ft. Each Ft is n× pt and is obtained from Ft−1 by adding a column on the right and, if the resulting number of columns is greater than p, also cleaning up (re-initialize) the table. That is,we never need to delete the left column because cleaning up the table stands for a restarted version of AM. As a result, we only need to handle two cases; 1 the table is empty(cleaned). 2 the table is not full. When the table is empty, we initialize F1 = Q1R1 with Q1 = 4f0/ ‖4f0‖2 and R = ‖4f0‖2. If the table size is smaller than p, we add a column on the right of Ft−1. Have Ft−1 = QR, we update Q and R so that Ft = [Ft−1,∆ft−1] = QR. It is a single modified Gram–Schmidt sweep that is described as follows: Algorithm 5: QR-updating procedures for i = 1, . . . , pt−1 do
Set R(i, pt) = Q(:, i)′ ∗ 4ft−1. Update4ft−1 ← ∆ft−1 −R (i, pt) ∗Q(:, i)
end Set Q (:, pt) = 4ft−1/ ‖4ft−1‖2 and R (pt, pt) = ‖∆ft−1‖2
Note that we do not explicitly conduct QR decomposition in each iteration, instead we update the factors (O(p2n)) and then solve a linear system using back substitution which has a complexity of
O(p2). Based on this complexity analysis, we can find Anderson Mixing with QR-updating scheme has limited computational overhead than GDA (or OG). This explains why GDA-AM is faster than EG but slower than OG in terms of running time of each iteration.

heading: B THEORETICAL RESULTS

heading: B.1 DIFFICULTY OF ANALYSIS ON GDA WITH ANDERSON MIXING
text: In the analysis, we study the inherent structures of the dynamics of the fixed point iteration and provide the convergence analysis for both simultaneous and alternating schemes. We want to emphasize that the direct application of existing convergence results of GMRES can not lead to convergent results. A recent paper Bollapragada et al. (2018) study the convergence acceleration schemes for multi-step optimization algorithms using Regularized Nonlinear Acceleration. We also want to point out that a naı̈ve application of Crouzeix’s bound to the minimax optimization problem can not be used to derive the convergent result. Theorem B.1 (Fischer & Freund (1991)). Let n ≥ 5 be an integer, r > 1, and c ∈ R. Consider the following constrained polynomial minmax problem
min p∈Pn:p(c)=1 max z∈Er
|p(z)| (11)
where
Er :=
{ z ∈ C ∣∣∣ |z − 1|+ |z + 1| ≤ r + 1 r } (12)
and c ∈ C \ Er. Then this problem can be solved uniquely by
tn(z; c) := Tn(z)
Tn(c) , (13)
where
Tn(z) = 1
2
( vn + 1
vn
) , z = 1
2
( v + 1
v
) (14)
if (a) |c| ≥ 12 ( r √ 2 + r− √ 2 ) or
(b) |c| ≥ (1/2ar) ( 2a2r − 1 + √ 2a4r − a2r + 1 ) , where ar := 12 ( r + 1r ) .
This is because the point 0 where all the residual polynomials take the fixed value of 1 is included in the numerical range of the iteration matrix, which violates the assumption of Theorem B.1. As a result, it can not be used to prove that the residual norm is decreasing based on this approach. Instead, we show that although the coefficient matrix is non-normal, it is diagonalizable. We then give the convergence results based on the eigenvalues instead of the numerical range. More specifically, Anderson mixing is equivalent to GMRES being applied to solve the following linear system:
(I−G(Alt))w = b(Alt), with w0 = w(Alt)0 . (15) Writing this linear system in the block form:[
0 ηA −ηAT η2ATA
] w = b(Alt). (16)
The residual norm bound for GMRES reads:
‖rt‖2 = min p∈P1t ‖p(I−G(Alt))r0‖2. (17)
Notice that the matrix (I −G(Alt)) is non-normal. If we apply Crouzeix’s bound in Crouzeix & Palencia (2017) to our problem as Bollapragada et al. (2018) did, then we have the following bound
‖rt‖2 ‖r0‖2 ≤ min p∈P1t ‖p(I−G(Alt))‖ ≤ (1 + √ 2) min p∈P1t sup z∈W (I−G(Alt)) ‖p(z)‖ (18)
where W (I −G(Alt)) = {z∗(I −G(Alt))z,∀z ∈ C2n \ {0}, ‖z‖ = 1} is the numerical range for I−G(Alt). In order to simplify the upper bound in the previous theorem, we study the numerical
range of I −G(Alt) similar to Bollapragada et al. (2018). Writing z = [ z1 z2 ] and computing the
numerical range of I−G(Alt) explicitly yields:
[z∗1, z ∗ 2]
[ 0 ηA
−ηAT η2ATA ] [ z1 z2 ] = η2z∗2A TAz2 + ηz ∗ 1Az2 − ηz∗2AT z1. (19)
For a general matrixA, there is no special structure about the numerical range of I−G(Alt). However, when A is symmetric, we can decompose A as A = ∑n i=1 λiviv T i where {λi}ni=1 are eigenvalues
of A in decreasing order and {vi}ni=1 are associated eigenvectors, and write ATA = ∑n i=1 λ 2 iviv T i . Then we can compute the numerical range of G(Alt) as follows: n∑ i [z∗1, z ∗ 2] [ 0 ηλiviv T i −ηλivivTi η2λ2ivivTi ] [ z1 z2 ] = n∑ i [z∗1vi, z ∗ 2vi] [ 0 ηλi −ηλi η2λ2i ] . [ vTi z1 vTi z2 ] (20)
Following the techniques proposed in Bollapragada et al. (2018) to analyze the numerical range of general 2× 2 matrices, we can show that the numerical range of I−G(Alt) is equal to the convex hull of the union of the numerical range of
Gi = [ 0 ηλi −ηλi η2λ2i ] , i = 1, . . . , n. (21)
And the boundary of numerical range of Gi is an ellipse whose axes are the line segments joining the points x to y and w to z, respectively, with
x = 0, y = η2λ2i , , w = η2λ2i 2 − √ −1η|λi|, z = η2λ2i 2 + √ −1η|λi|. (22)
Thus, the numerical range of I−G(Alt) can be spanned by convex hull of the union of the numerical range of a set of 2-by-2 matrices and the numerical range of each such a 2-by-2 matrix is an ellipse. We can compute the center o and focal distance d of the ellipse generated by numerical range of I−G(Alt) explicitly. Then a linear transformation enables us to use Theorem B.1 to show that the near-best polynomial for the minimax problem on the numerical range of I −G(Alt) is given by tn(z; c) := Tn( z−o d )
Tn( c−o d ) if 0 is excluded from the numerical range of I−G(Alt). However, according to equation 22 the numerical range includes the point 0 where the residual polynomial takes value 1, thus the analysis based on numerical range can not help derive the convergent result as the upper bound is not guaranteed to be less than 1.

heading: B.2 PROOFS OF THEOREM
text: We first provide proof of Theorem 4.1. Theorem B.2 (Global convergence for simultaneous GDA-AM on bilinear problem). Denote the distance between the stationary point w∗ and current iterate w(k+1)p of Algorithm 2 with Anderson restart dimension p as N(k+1)p = dist(w∗,w(k+1)p). Then we have the following bound for Nt Algorithm 2 is unconditionally convergent
N(k+1)p ≤ 1
Tp(1 + 2 κ(ATA)−1 ) Nkp (23)
where Tp is the Chebyshev polynomial of first kind of degree p and 1Tp(1+ 2 κ(ATA)−1 ) < 1 since 1 + 2 κ(ATA)−1 > 1.
Proof of Theorem 4.1. Note that I −G(Sim) is a normal matrix which will be denoted as G for notational simplicity. Thus it admits the following eigendecomposition:
G = UΛUT , UUT = I, Λ = diag(λ1, . . . , λ2n). (24)
Based on the equivalence between GMRES and Anderson Mixing, we know that the convergence rate of simultaneous GDA-AM can be estimated by the spectrum of G. Especially, it holds that
r(k+1)p = Ufp(Λ)U T rkp. fp ∈ Pp (25)
where Pp is the family of residual polynomials with degree p such that fp(0) = 1,∀fp ∈ Pp. According to Lemma 2.1, we have the following estimation
‖r(k+1)p‖2 = min fp∈Pp ‖fp(G)rkp‖2 ≤ min fp∈Pp max i |fp(λi)|‖rkp‖2. (26)
Due to the block structure of G, the eigenvalues of G can be computed explicitly as
± ησi √ −1, i = 1, . . . , n, (27)
where σi is the ith largest singular value of matrix A. This shows that the eigenvalues of G are n pairs of purely imaginary numbers excluding 0 since A has full rank.
Since the eigenvalues of G are distributed in two intervals excluding the origin
I = [−ησmax(A) √ −1,−ησmin(A) √ −1] ∪ [ησmin(A) √ −1, ησmax(A) √ −1],
it can be shown that the following p-th degree polynomial with value 1 at the origin that has the minimal maximum deviation from 0 on I is given by:
fp(z) = Tl(q(
√ −1z))
Tl(q(0)) , q(
√ −1z) = 1− 2(
√ −1z − ησmin)( √ −1z + ησmin)
(ησmax(A))2 − (ησmin(A))2 (28)
where l = [p2 ] and Tl is the Chebyshev polynomial of first kind of degree l. The function q( √ −1z) maps I to [−1, 1]. Thus the numerator of the polynomial fp is bounded by 1 on I. The size of denominator can be determined by the method discussed in Chapter 3 of Greenbaum (1997). Assume q(0) = 12 (y + y −1), then Tl(q(0)) = 12 (y l + y−l). Then y can be determined by solving
q(0) = (ησmax(A))
2 + (ησmin(A)) 2
(ησmax(A))2 − (ησmin(A))2 . (29)
The solutions to this equation are
y1 = ησmax(A) + ησmin(A)
ησmax(A)− ησmin(A) or y2 = ησmax(A)− ησmin(A) ησmax(A) + ησmin(A) . (30)
Then plugging the value of q(0) into the polynomial fp yields
‖r(k+1)p‖ ‖rkp‖
≤ 2 (√η2σ2max(A)−√η2σ2min(A)√
η2σ2max(A) + √ η2σ2min(A) )l = 2 (σmax(A)− σmin(A) σmax(A) + σmin(A) )l = 2 (κ(A)− 1 κ(A) + 1
)l (31) Note that Nt and rt is related through G(wt −w∗) = rt. Therefore,
N(k+1)p = ‖w(k+1)p −w∗‖2 = ‖G−1r(k+1)p‖2 = min fp∈Pp ‖G−1fp(G)G(wkp −w∗)‖2
≤ min fp∈Pp max i |fp(λi)|‖wkp −w∗‖2 ≤ 2
( 1− 2
κ(A) + 1
) p 2
Nkp. (32)
Actually a tighter bound can be proved after noting that the problem is essentially equivalent to polynomial minmax problem on the interval:
Ĩ = [η2σ2min(A), η 2σ2max(A)],
Then it is well known that,
N(k+1)p ≤ min fp∈Pp max λi∈[η2σ2min(A), η2σ2max(A)]
|fp(λi)|‖wkp −w∗‖2 ≤ 1
Tp(1 + 2 σ2min σ2max−σ2min ) Nkp
≤ 1 Tp(1 + 2 κ(ATA)−1 ) Nkp
(33)
where Tp Chebyshev polynomial of degree p of the first kind and 1Tp(1+ 2 κ(ATA)−1 ) < 1. Explicitly,
Tp(1 + 2
κ(ATA)− 1 ) =
1
2
[( 1 +
2
κ(ATA)− 1 +
√ (1 +
2
κ(ATA)− 1 )2 − 1 )p + ( 1 + 2
κ(ATA)− 1 +
√ (1 +
2
κ(ATA)− 1 )2 − 1
)−p]
Next, we give the proof of Theorem 4.2. Theorem B.3 (Global convergence for alternating GDA-AM on bilinear problem). Denote the distance between the stationary point w∗ and current iterate w(k+1)p of Algorithm 3 with Anderson restart dimension p as N(k+1)p = dist(w∗,w(k+1)p). Assume A is normalized such that its largest singular value is equal to 1. Then when the learning rate η is less than 2, we have the following bound for Nt
N2(k+1)p ≤ √ 1 + 2η
2− η ( r c )pN2kp
where c and r are the center and radius of a disk D(c, r) which includes all the eigenvalues of G in equation 4.3. Especially, rc < 1.
Proof. Since the residual rp of AA at p-th iteration has the form of
rp = (I− p∑ i=1 Gi)r0,
and AA minimizes the residual, we have
‖r(k+1)p‖22 ≤ min β ‖rkp − βGirkp‖22 ≤ min fp∈Pp ‖fp(G)rkp‖22,
where Pp is the family of polynomials with degree p such that fp(0) = 1,∀fp ∈ Pp . It’s easy to see that G is unitarily similar to a block diagonal matrix Λ with 2× 2 blocks as follows:[
0 ησi −ησi (ησi)2
] ∀ i ∈ [n].
Thus the eigenvalues of G can be easily identified as
λ±i = (ησi(ησi ±
√ (ησi)2 − 4))
2 , i ∈ [n].
where σ1 ≥ σ2 ≥ · · · ≥ σn are the singular values of A. Furthermore, the eigenvector and eigenvalue associated with each 2× 2 diagonal block are[
0 ησi −ησi (ησi)2 ] [ 1 λ±i ησi ] = λ±i [ 1 λ±i ησi ] Thus G is diagonalizable and denote the matrix with the columns of eigenvectors of G by X. The real part of the eigenvalues of G are at least
R(λ±i) ≥ (ησi)
2
2 , i ∈ [n]. (34) And since |ησi| ≥ | √
(ησi)2 − 4)|, all the eigenvalues will be included in a disk D(c, r) which is included in the right half plane. Moreover, both c and r being greater than zero indicates that rc < 1. Start from the following inequality:
N(k+1)p = ‖w(k+1)p −w∗‖2 = ‖G−1r(k+1)p‖2 ≤ min fp∈Pp ‖G−1fp(G)rkp‖2
= min fp∈Pp ‖G−1fp(G)G(wkp −w∗)‖2 = min fp∈Pp ‖G−1p (G)(wkp −w∗)‖2
= min fp∈Pp
‖fp(G)(wkp −w∗)‖2.
We will use the eigendeomposition of G and the special polynomial ( c−tc ) p to derive the inequality in Theorem 3. Now we know rc < 1. If we choose gp(t) = ( c−t c ) p, we can obtain
min fp∈Pp
‖fp(G)(wkp −w∗)‖2 ≤ ‖gp(G)(wkp −w∗)‖2
which implies
min fp∈Pp
‖fp(G)(wkp −w∗)‖2 ≤ ‖gp(XΛX−1)‖‖(wkp −w∗)‖2
Since G is diagonalizable (which has been shown above), we assume the eigendecomposition of G is G = XΛX−1. Then
min fp∈Pp ‖gp(G)(wkp−w∗)‖2 ≤ ‖X‖‖X−1‖ max {λi}2ni=1
‖gp(Λ)‖‖(wkp−w∗)‖2 ≤ κG( r
c )p‖wkp−w∗‖2
where κG is the condition number of X . The last inequality comes from Lemma 6.26 and Proposition 6.32 in Saad (2003).. Since G and Λ are unitarily similar, κG is equal to the condition number of the eigenvector matrix of Λ. The eigenvector matrix of Λ is a block diagonal matrix with the ith block
as [
1 1 λ+i ησi λ−i ησi
] . Thus the singluar values of the eigenvector matrix of Λ is equal to the union of the
singular values of these 2-by-2 blocks. Under the assumption that the largest singular value of A are equal to 1 and the learning rate is less than 2, it is easy to find the singular values of the eigenvector matrix of Λ are √ 2± ησi. Thus, κG = √
2+ησmax√ 2−ησmax
= √
2+η√ 2−η = √ 1 + 2η2−η .

heading: B.3 DISCUSSION OF OBTAINED RATES
text: We would like to first explain on why taking Chebyshev polynomial of degree p at the point 1 + 2κ−1 . We evaluate the Chebyshev polynomial at this specific point because the reciprocal of this value gives the minimal value of infinite norm of the all polynomials of degree p defined on the interval Ĩ = [η2σ2min(A), η
2σ2max(A)] based on Theorem 6.25 (page 209) (Saad, 2003). In other words, taking the function value at this point leads to the tight bound.
When comparing between existing bounds, we would like to point our our derived bounds are hard to compare directly. Alternatively, we can derive another bound for comparison with existing bounds for simultaneous GDA-AM. If we use the inequality that Tp(t) ≥ 12 ((t+ √ t2 − 1)p), we can obtain
the bound ρ(A) = 4( √ κ(ATA)−1√ κ(ATA)+1 )2 = 4(1−O( 1√ κ(ATA) )), which is in a form that is comparable with EG and can compete with EG + positive momentum. The numerical experiments in figure 2b numerically verify that our bound is smaller than EG. We wanted to numerically compare our rate with EG with positive momentum. However the bound of EG with positive momentum is asymptotic. Moreover, it does not specify the constants so we can not numerically compare them. We do provide empirical comparison between GDA-AM and EG with positive momentum for bilinear problems in Appendix C.1. It shows GDA-AM outperforms EG with positive momentum. Regarding alternating GDA-AM , we would like to note that the bound in Theorem 4.2 depends on the eigenvalue distribution of the matrix G. Condition number is not directly related to the distribution of eigenvalues of a nonsymmetric matrix G. Thus, the condition number is not a precise metric to characterize the convergence. If these eigenvalues are clustered, then our bound can be small. On the other hand, if these eigenvalues are evenly distributed in the complex plane, then the bound can very close to 1.
More importantly, we would like to stress several technical contributions.
1. Our obtained Theorem 4.1 and 4.2 provide nonasymptotic guarantees, while most other work are asymptotic. For example, EG with positive momentum can achieve a asymptotic rate of 1−O(1/ √ κ) under strong assumptions (Azizian et al., 2020).
2. Our contribution is not just about fix the convergence issue of GDA by applying Anderson Mixing; another contribution is that we arrive at a convergent and tight bound on the original work and not just adopting existing analyses. We developed Theorem 4.1 and 4.2 from a new perspective because applying existing theoretical results fail to give us neither convergent nor tight bounds.
3. Theorem 4.1 and 4.2 only requires mild conditions and reflects how the table size p controls the convergence rate. Theorem 4.1 is independent of the learning rate η. However, the convergence results of other methods like EG and OG depend on the learning rate, which may yield less than desirable results for ill-specified learning rates.

heading: B.4 CONVEX-CONCAVE AND GENERAL CASE
text: Given the widespread usage of minimax problems in applications of machine learning, it is natural to ask about its properties when being applied to general nonconvex-nonconcave settings. If f is a nonconvex-nonconcave function, the problem of finding global Nash equilibrium is NP-hard in general. Recently, Jin et al. (2020) show that local or global Nash equilibrium may not exist in nonconvex-nonconcave settings and propose a new notation local minimax as defined below:
Definition 4. A point (x?,y?) is said to be a local minimax point of f , if there exists δ0 > 0 and a function h satisfying h(δ)→ 0 as δ → 0, such that for any δ ∈ (0, δ0], and any (x,y) satisfying ‖x− x?‖ ≤ δ and ‖y − y?‖ ≤ δ, we have
f (x?,y) ≤ f (x?,y?) ≤ max y′:‖y′−y?‖≤h(δ) f (x,y′) .
Jin et al. (2020) also establishes the following first- and second-order conditions to characterize local minimax:
Proposition 1 (First-order Condition). Any local minimax point (x∗,y∗) satisfies∇f(x∗,y∗) = 0. Proposition 2 (Second-order Necessary Condition). Any local minimax point (x∗,y∗) satisfies ∇yyf(x∗,y∗) 4 0 and∇xxf(x∗,y∗)−∇xyf(x∗,y∗)(∇yyf(x∗,y∗))−1∇yxf(x∗,y∗) < 0. Proposition 3 (Second-order Sufficient Condition). Any stationary point (x∗,y∗) satisfies ∇yyf(x∗,y∗) ≺ 0 and ∇xxf(x∗,y∗)−∇xyf(x∗,y∗)(∇yyf(x∗,y∗))−1∇yxf(x∗,y∗) 0 is a local minimax point.
Given the second-order conditions of local minimax, it turns out that above question is extremely challenging—GDA-AM is a first-order method. But we can prove the following result for GDA-AM:
Theorem B.4 (Local minimax as subset of limiting points of GDA-AM). Consider a general objective function f(x,y). The set of limiting points of GDA-AM for minimax problem
min x∈Rn max y∈Rn f(x,y)
includes the local minimax points of this function.
The definition of local minimax is stronger than that of first order point. The convergence analysis for complexity of finding stationary point is included in the next section. The proof of Theorem B.4 needs the result from the following theorem.
Theorem B.5 (Calvetti et al. (2002)). Let δ satisfy 0 < δ ≤ δ0 for some constant δ0 > 0 (refer to Calvetti et al. (2002) for details), and let bδ ∈ X satisfy ∥∥b− bδ∥∥ ≤ δ. Let k ≤ ` and let xδk denote the kth iterate determined by the GMRES method applied to equation Ax = bδ, with initial guess xδ0 = 0. Similarly, let xk denote the kth iterate determined by the GMRES method applied to equation Ax = b with initial guess x0 = 0. Then, there are constants σk independent of δ, such that∥∥xk − xδk∥∥ ≤ σkδ, 1 ≤ k ≤ ` Then, we give the proof of Theorem B.4.
Proof of Theorem B.4. For notational simplicity, we will denote ∇xxf(x∗,y∗), ∇xyf(x∗,y∗) and ∇yyf(x∗,y∗) by Hx∗x∗ , Hx∗y∗ and Hy∗y∗ , respectively. Simultaneous GDA can be written as
wt+1 = [ xt+1 yt+1 ] = [ xt − η∇xf(xt,yt) yt + η∇yf(xt,yt) ] .
Since the function is differentiable, Taylor expansion holds for ∇xf(xt,yt) and ∇yf(xt,yt) at a local minimx point w∗ = (x∗,y∗),
∇xf(xt,yt) = ∇xf(x∗,y∗) + Hx∗x∗(xt − x∗) + Hx∗y∗(yt − y∗) + o(‖wt −w∗‖2) ∇yf(xt,yt) = ∇yf(x∗,y∗) + Hy∗y∗(yt − y∗) + Hy∗x∗(xt − x∗) + o(‖wt −w∗‖2).
Use the fact that∇f(x∗,y∗) = 0 to simplify the above equations and obtain
∇xf(xt,yt) = Hx∗x∗(xt − x∗) + Hx∗y∗(yt − y∗) + o(‖wt −w∗‖2) ∇yf(xt,yt) = Hy∗y∗(yt − y∗) + Hy∗x∗(xt − x∗) + o(‖wt −w∗‖2).
Inserting the above formulas into the iteration scheme, it yields
wt+1 = [ xt+1 yt+1 ] = [ I− ηHx∗x∗ −ηHx∗y∗ ηHy∗x∗ I + ηHy∗y∗ ] [ xt yt ] + [ ηHx∗x∗x ∗ + ηHx∗y∗y ∗ + −ηHy∗y∗y∗ − ηHx∗y∗x∗ + ] where denotes the higher order error o(‖wt −w∗‖2). According to Theorem 2.2, we know that simultaneous GDA-AM is equivalent to applying GMRES to solve the following linear system
(I− [ (1− α)I− ηHx∗x∗ −ηHx∗y∗
ηHy∗x∗ (1− α)I + ηHy∗y∗
] )w = [ αI + ηHx∗x∗ ηHx∗y∗ −ηHy∗x∗ αI− ηHy∗y∗ ] w = b+
where b = [ ηHx∗x∗x ∗ + ηHx∗y∗y ∗ −ηHy∗y∗y∗ − ηHx∗y∗x∗ ] . We now know that GDA-AM is equivalent to GMRES being applied to solve the following linear system[ αI + ηHx∗x∗ ηHx∗y∗ −ηHy∗x∗ αI− ηHy∗y∗ ] w̃ = b
The symmetric part of the coefficient matrix of the above linear system is[ αI + ηHx∗x∗ 0
0 αI− ηHy∗y∗
] .
According to Proposition 2, αI − ηHy∗y∗ is positive definite since Hy∗y∗ 4 0. If Hx∗x∗ is positive semidefinite, then αI + ηHx∗x∗ is positive definite and we’re done. Otherwise, assume λmin(Hx∗x∗) < 0. Then for fixed α, when η < − αλmin(Hx∗x∗ ) , αI+ηHx∗x∗ will be positive definite. Then according to Theorem 2.2, we know GDA-AM indeed converges. Let’s create a new companion linear system as follows [
αI + ηHx∗x∗ ηHx∗y∗ −ηHy∗x∗ αI− ηHy∗y∗
] ŵ = b + αw∗
Note that ŵ = w∗ and GMRES on this companion linear system is convergent under suitable choice of learning rate η. Let the iterates of GMRES for w̃, ŵ,w be denoted by w̃t, ŵt,wt. Then ‖w̃t − ŵt‖ ≤ ‖w̃t −wt‖+ ‖ŵt −wt‖. According to Theorem B.5, we also have ‖w̃t −wt‖ ≤ σk , 1 ≤ k ≤ t. Further more, again according to Theorem B.5, we know ‖ŵt−ŵt‖ ≤ σk(αw∗+ ). Starting from an initial point very close to w∗ and let t → ∞ and α, → 0, ŵt will converge to w∗ = (x∗,y∗), which means the local minimax w∗ = (x∗,y∗) is a limiting point of GDA-RAM.
Theorem B.6. For strongly-convex-strongly-concave function f(x,y), GDA-AM will converge to the Nash equilibrium of this function.
Proof: Since strongly-convex-strongly-concave function f(x,y) has unique Nash equilibrium which is also the unique minimax point, this minimax point must be the limiting point of GDA-AM according to Theorem B.4.

heading: B.4.1 BILINEAR-QUADRATIC GAMES
text: Moreover, we can further show that the GDA-AM converges on bilinear-quadratic games. Consider a quadratic problem as follows,
min x∈Rn max y∈Rn
f(x,y) = xTAy + xTBx− yTCy + bTx + cTy, (35)
where A is full rank, B and C are both positive definite.
Theorem B.7. [Global convergence for simultaneous GDA-AM on bilinear-quadratic problem] Let r
(Sim) t be the residual of Algorithm 2 being applied to problem equation 35. For some constant ρ < 1,
‖r(Sim)t ‖2 ≤ ( 1− (λmin(J T + J))2
4λmax (JTJ) )t/2 ︸ ︷︷ ︸
ρt/2
‖r0‖2, (36)
where J = [ ηB ηA −ηAT ηC ] and λmin and λmax denote the smallest and largest eigenvalue, respec-
tively.
The convergence property of GMRES has been studied in the next theorem. We use this theorem to show the convergence rate of GDA-AM for bilinear-quadratic games. Theorem B.8 (Elman (1982)). Consider solving a linear system Ex = b using GMRES. Let rt = b−Ext be the residual at tth iteration. If the Hermitian part of E is positive definite, then for some positive constant ρ < 1, it holds that
‖rt‖2 ≤ ( 1− (λmin(E H + E))2
4λmax (EHE) )t/2 ︸ ︷︷ ︸
ρt/2
‖r0‖2. (37)
Proof of Theorem B.7. Applying simultaneous GDA-AM to solve the above problem is equivalent to applying Anderson Mixing on the following fixed point iteration:[
xt+1 yt+1
] = [ I− ηB −ηA ηAT I − ηC ] ︸ ︷︷ ︸
G(Quad−sim)
[ xt yt ] ︸︷︷︸
w (Quad−sim) t
+ [ −ηb −ηc ] ︸ ︷︷ ︸
b (Quad−sim)
. (38)
We know that we need to study the convergence properties of GMRES for solving the following linear system [
ηB ηA −ηAT ηC
] w = b. (39)
For notational simplicity, the superscripts has been dropped. Denote the coefficient matrix[ ηB ηA −ηAT ηC ] by J. The symmetric part of J is
J + JT
2 =
[ η 2 (B + B
T ) 0 0 η2 (C + C T ) ] which is positive definite. Then immediately by Theorem B.8, the following convergence rate holds For some constant 0 < ρ < 1,
‖rt‖2 = min p∈P1t ‖p(J)r0‖2 ≤
( 1− (λmin(J + J T ) 2
(4λmax (JTJ)) )t/2 ︸ ︷︷ ︸
ρt/2
‖r0‖2
= ρt/2‖r0‖2
(40)
Note that the convergence of GDA-AM for bilinear-quadratic games can also be analyzed by numerical range as shown in (Bollapragada et al., 2018). Although we previously show that analysis based on the numerical range can not help us derive a convergent bound for bilinear games, we show analysis in Bollapragada et al. (2018) can be extended to bilinear-quadratic games. When B and C are positive definite, 1 is outside of the numerical range of matrix G(Quad−sim) as shown in 7a. When B or C is not positive definite, 1 can be included in the numerical range of matrix G(Quad−sim) as shown in 7b. That is saying analysis based on the numerical range (Crouzeix & Palencia, 2017; Bollapragada et al., 2018) to the bilinear-quadratic problem can lead to a convergent result when B and C are positive definite. And analysis based on the numerical range can not help us derive convergent results when B or C is not positive definite.
-40 -30 -20 -10 0 10
-6
-4
-2
0
2
4
6
1
(a) Positive definite B and C -4 -2 0 2 4 6
-10
-5
0
5
10
1
(b) Random generated B and C
Figure 7: Numerical range of fixed-point operator (Simultaneous GDA-AM ) G =[ I− ηB −ηA
ηAT I − ηC
]
for bilinear-quadratic games.

heading: B.5 STOCHASTIC CONVEX-NONCONVACE CASE
text: In this section, we study the convergence of GDA-AM for convex-noncovace problem in the stochastic setting with the same assumptions in Wei et al. (2021b); Xu et al. (2021). The recent work Wei et al. (2021b) proves the convergence of the stochastic gradient descent with Anderson Mixing for min optimization. The convergence of GDA-AM for minimax optimization builds on top of it with several modifications. The minimax problem is equivalent to minimizing a function Φ(·) = maxy∈Y f(·,y) (Lin et al., 2020). And we are interested in complexity of a pair of -stationary point (x, y) instead of analysis of a point x.
Definition 5. (Lin et al., 2020) A pair of points (x,y) is an -stationary point ( ≥ 0) of a differentiable function Φ if
‖∇xf(x,y)‖ ≤ ‖PY (y + (1/`)∇yf(x,y))− y‖ ≤ /`
Assumption 1. f : Rd 7→ R is continuously differentiable. f(x) ≥ f low > −∞ for any x ∈ Rd. ∇f is globally L-Lipschitz continuous; namely ‖∇f(x)−∇f(y)‖2 ≤ L‖x− y‖2 for any x, y ∈ Rd.
Assumption 2. For any iteration k, the stochastic gradient ∇fξk (xk) satisfies Eξk [∇fξk (xk)] = ∇f (xk) ,Eξk [ ‖∇fξk (xk)−∇f (xk)‖ 2 2 ] ≤ σ2, where σ > 0, and ξk, k = 0, 1, . . ., are indepen-
dent samples that are independent of {xi}k
Theorem B.9. For a general convex-nonconcave function f , suppose that Assumptions 1 and 2 hold. Batch size nt = n for t = 0, . . . , N − 1. C > 0 is a constant. βt = µ4L(1+C−1) .δt ≥
Cβ−2t , 0 ≤ αt ≤ min { 1, β 1 2 t } and αt is chosen to make sure the positive definiteness of Ht. Let
R be a random variable following PR(t) def = Prob{R = t} = 1/N , and N̄ be the total number of stochastic GDA-AM calls needed to calculate stochastic gradients ∇̃fSt (wt) in our algorithm. To ensure E [∥∥∥∇̃f (wR)∥∥∥ 2 ] ≤ , total number of stochastic GDA-AM calls needed to calculate stochastic gradients ∇̃fSt (wt) is O( −4).
Recall that we can recast GDA scheme as the following fixed point iteration.
wt+1 = G (sim) η (wt) , wt + ηV (wt) with w = [ x y ] , V (w) = [ −∇xf(x,y) ∇yf(x,y) ] Ignoring the stepsize η and let Wt and Rt record the first and second order diffrence of recent m iterates:
Wt = [∆wt−m,∆wt−m+1, · · · ,∆wt−1] ,Rt = [∆Vt−m,∆Vt−m+1, · · · ,∆Vt−1] Similarly as Wei et al. (2021b),the Anderson mixing can be decoupled into
w̄t+1 = wt −WtΓt, (Projection step) w̄t+1 = wt + βtV̄t, (Mixing step)
where βt is the mixing parameter, and V̄t = Vt −WtΓt and Γt is solved by Γt = arg min
Γ∈Rm ‖Vt −RtΓ‖2 + δt‖Γ‖2
We want to argue that similar arguments in Wei et al. (2021b) can be applied to the problem here. To see why Anderson mixing works for minimax optimization, we assume function f is smooth. Then the hessian matrix for G(sim)η is
H = ( −∇2xxf −∇2xyf ∇2yxf ∇2yyf ) Notice that in a small neighborhood of wt+1, we have
Rt = −HWt = ( ∇2xxf ∇2xyf −∇2yxf −∇2yyf ) Wt
Thus ‖Vt − RtΓ‖2 ≈ ‖Vt + HWtΓ‖2, which is equivalent to solving for a vector pt such that Hpk = Vt. This is exactly the second order method for the fixed point iteration problem. Also at each step the AM is minimizing the residual, the reason that AM is equivalent to GMRES for linear problem is that this quadratic approximation is exact. Finally, we rewrite AM as the quasi-newton framework as Wei et al. (2021b) did. wt+1 = wt +HtVt where
min Ht ‖Ht − βtI‖F subject to HtRt = −Xt
Finally, with damping parameter, Anderson mixing has the following form
Wt+1 = Wt + βtVt − αt (Wt + βtRt) Γt (41)
we can also apply the very similar arguments to prove key results in lemma 1, lemma 2 in Wei et al. (2021b). There is also a key difference with Wei et al. (2021b). Here we are considering minimax optimization problem. Thus our gradient is actually V (w) = [ −∇xf(x,y) ∇yf(x,y) ] rather than
∇f(w) = [ ∇xf(x,y) ∇yf(x,y) ] This will introduce some difficulty to the dynamics of the fixed pointe iteration. However, noticing that ‖V ‖ = ‖∇f(w)‖ and
f (wt+1) ≤ f (wt) +∇f (wt)T (wt+1 −wt) + L
2 ‖wt+1 −wt‖22
≤ f (wt) + ∇̃f (wt)T (wt+1 −wt) + L
2 ‖wt+1 −wt‖22
= f (wt) + ∇̃f (wt)THtVt + L
2 ‖HtVt‖22
(42)
where
∇̃f (wt) = [ −∇xf(x,y) ∇yf(x,y) ] (43)
we call this the ascent-descent gradient (ADG) which is the gradient for minimax optimization problem
min x∈Rd max y∈Rd f(x,y).
To see why ∇f (wt)T (wt+1 −wt) ≤ ∇̃f (wt)T (wt+1 −wt), we consider their difference (∇̃ − ∇)f (wt)T (wt+1 −wt) = −2∇xf(xt,yt)T (xt+1 − xt).
For fixed yt, f(xt,yt+1) has the Talyor expansion:
f(xt+1,yt) = f(xt,yt)+∇xf(xt,yt)T (xt+1−xt)+(xt+1−xt)T∇xxf(xt+θ(xt+1−xt),yt)(xt+1−xt) Assuming f is convex w.r.t x and apply safeguard to ensure f(xt+1,yt) ≤ f(xt,yt) can guarantee (∇̃ − ∇)f (wt)T (wt+1 −wt) ≥ 0. Now applying lemmas in Wei et al. (2021b), we can derive the convergence of our method for general convex-nonconcave function similarly.

heading: C ADDITIONAL EXPERIMENTS

heading: C.1 COMPARISON WITH EG WITH POSITIVE MOMENTUM
text: In this section, we include additional comparison between GDA-AM and EG with positive momentum. GDA-AM has two big theoretical advantages over EG with positive momentum. First, convergence of GDA-AM does not require strong assumptions on choices of hyperparamters. Second, 4.1 and 4.2 provide nonasymptotic guarantees while convergence of EG with positive mometum is asymptotic. Experimental results are shown in 8. It indicates GDA-AM outperforms EG with positive momentum. Finding a good choice of the inner and outer step size of EG and momentum term is hard. For EG with positive momentum, we set the step size of extrapolation step as 1, the step size of update as 0.5, and the positive momentum term as 0.3 after grid search as shown in 8b and 8c. On the other hand, GDA-AM converges fast for different step size without hyperparameter tuning.

heading: C.2 ROBUST NEURAL NETWORK TRAINING
text: In this section, we test the effectiveness of GDA-AM by training a robust neural network on MNIST data set against adversarial attacks (Madry et al., 2019; Goodfellow et al., 2015; Kurakin et al., 2017) . The optimization formulation is
min w N∑ i=1 max δi, s.t. |δi|∞≤ε ` (f (xi + δi; w) , yi) (44)
where w is the parameter of the neural network, the pair (xi, yi) denotes the i-th data point, and δi is the perturbation added to data point i. The accuracy of our formulation against popular attacks, FGSM (Goodfellow et al., 2015) and PGD (Kurakin et al., 2017), are summarized in Table 2.. Since solving such problem is computationally challenging, Nouiehed et al. (2019) proposed an approximation of the above optimization problem with a new objective function as the following nonconvex-concave problem:
min w N∑ i=1 max t∈T 9∑ j=0 tj` ( f ( xKij ; w ) , yi ) , T =
{ (t1, · · · , tm) |
m∑ i=1 ti = 1, ti ≥ 0
} (45)
where K is a parameter in the approximation, and xKij is an approximated attack on sample xi by changing the output of the network to label j. We use the public available implementation (Nouiehed et al., 2019) 1. We apply our algorithm on top of (Nouiehed et al., 2019) and compare our results (p = 50) with (Madry et al., 2019; Zhang et al., 2019; 2020; Nouiehed et al., 2019). Results are summarized in table 2. We can observe that GDA-AM leads to a comparable or slightly better performance to the other methods. In addition, GDA-AM does not exhibit a significant drop in accuracy when is larger and this suggests the learned model is more robust.
C.3 1D MINIMAX FUNCTIONS
We begin with investigating the empirical performance of GDA-AM for 6 non-trivial 1d bivariate functions. We set initial points as (3, 3) and m as 20 or 5 for all functions. We use optimal learning rates for all methods on each problem. Results are shown in Figure 9, 10, 11, 12, 13 and 14. We observe GDA-AM consistently outperforms all baselines and improves convergence. It is worthwhile to mention that the difference between GDA-AM and traditional averaging is twofold. First, traditional averaging does not involve an adaptive averaging scheme and thus blindly converge to (0, 0) for all 1d bivariate functions. In contrast, GDA-AM obtains optimal weights by solving a small linear system on past iterates. Using different weights for each iteration, GDA-AM is able to minimize the residual of past iterates and thus find the solution of a fixed-point iteration. More importantly, averaging does not change the GDA dynamic because averaging generates a new sequence of parameters based on GDA iterates. This means averaging is independent with base training algorithm (GDA here). However, GDA-AM changes the dynamic directly by overwriting the latest iterate. It means Anderson Mixing interacts with GDA, which is another major difference from averaging.

heading: C.4 DENSITY ESTIMATION
text: To test our proposed method, we evaluate our method on two low-dimension density estimation problems, mixture of 25 Gaussians and Swiss roll. For both generator and discriminator, we use fully connected neural networks with 3 hidden layers and 128 hidden units in each layer. Except for the output layer of discriminator that uses a sigmoid activation, we use tanh-activation for all other layers. We run Adam and GDA-AM for 50000 steps. The learning rate is set as 2×10−4 and β1 = 0, β2 = 0.9 after an extensive grid search, which is close to the maximal possible stepsize under which the methods rarely diverge. Figure 15 and 16 show the output after {1K, 10K, 30K, 50K} iterations. It can be seen that our method converges faster to the target distribution offers a improvement over Adam. In
1https://github.com/optimization-for-data-driven-science/Robust-NN-Training
addition, we can observe that the generated samples using our method gather around the circle and are less connected with other circles.
C.5 IMAGE GENERATION
In this section, we provide additional experimental results that are not given in Section 6. Figure 18a and 18b show the Inception Score for CIFAR10 using WGAN-GP and SNGAN. It can be observed that our method consistently performs better than Adam and EG during training. Further, on CIFAR-10 using WGAN-GP and SNGAN, GDA-AM is slightly slower than Adam (about 110-115% computational time), but significantly faster than EG (about 65-75% computational time).

heading: C.6 DETAILS ON THE EXPERIMENTS
text: For our experiments, we used the PyTorch 2 deep learning framework. Experiments were run one NVIDIA V100 GPU. The residual network architecture for generator and discriminator are summarized in Table 3 and 4. We use a WGAN-GP loss, with gradient penalty λ = 10. When using the gradient penalty (WGAN-GP), we remove the batch normalization layers in the discriminator. When using SNGAN, we replace the batch normalization layers with spectral normalization. Hyperparamters of Adam are selected after grid search. We use a learning rate of 2× 10−4 and batch size of 64. For table size of GDA-AM , we set it as 120 for CIFAR10 and 150 for CelebA. We set β1 = 0.0 and β2 = 0.9 as we find it gives us better models than default settings.
2https://pytorch.org/
Input: x ∈ R3×32×32 Linear 128→ 128× 4× 4
ResBlock 128→ 128 ResBlock 128→ 128 ResBlock 128→ 128
Linear 128→ 1
Input: x ∈ R3×64×64 Linear 128→ 128× 4× 4
ResBlock 128→ 128 ResBlock 128→ 256 ResBlock 256→ 512
Linear 512→ 1

Review summaries:
The authors propose to use the Anderson Acceleration on min max problem. They show two theoretical results: convergence rate on bilinear problems when using simultaneous gradient decent-ascent, and convergence rate on bilinear problems when using alternating gradient decent-ascent. Finally, they present numerical results in favor of their approach.

This paper focuses on solving minimax optimization using Anderson mixing. Anderson mixing is a framework to accelerate fixpoint iteration. The authors tried to use  Anderson mixing for the minimax optimization problem. In their framework, they use a weighted average of the solution of last p iterates and optimized over the weights to come up with a new point.  They show that for the bilinear case, this algorithm converges to the optimal stationary point. Finally, the authors simulated their method on synthetic examples as well as image generation on CIFAR10. 

Minimax optimization is a central problem in machine learning. The gradient descent ascent method is the most commonly used algorithm to solve this problem. This paper views the gradient descent ascent method as a fixed-point iteration and solves it using Anderson Mixing to converge to the local minimax. They show theoretically that the algorithm can achieve global convergence for bilinear problems under mild conditions. Some numerical experiments have been conducted. 

This paper proposes a new optimization method for min-max optimization (That could actually be generalized to any variational inequality problem) based on Anderson Mixing. They prove the convergence of their algorithm in the case of bilinear min-max games for both alternated and simultaneous updates operator. They eventually try their algorithm on toy bilnear problems and GANs.