PROVABLY CONVERGENT QUASISTATIC DYNAMICS FOR MEAN-FIELD TWO-PLAYER ZERO-SUM GAMES
Year: 2022

Abstract:
In this paper, we study the problem of finding mixed Nash equilibrium for meanfield two-player zero-sum games. Solving this problem requires optimizing over two probability distributions. We consider a quasistatic Wasserstein gradient flow dynamics in which one probability distribution follows the Wasserstein gradient flow, while the other one is always at the equilibrium. Theoretical analysis are conducted on this dynamics, showing its convergence to the mixed Nash equilibrium under mild conditions. Inspired by the continuous dynamics of probability distributions, we derive a quasistatic Langevin gradient descent method with inner-outer iterations, and test the method on different problems, including training mixture of GANs.

heading: 1 INTRODUCTION
text: Finding Nash equilibrium has seen many important applications in machine learning, such as generative adversarial networks (GANs) (Goodfellow et al., 2014a) and reinforcement learning (Busoniu et al., 2008). In these problems, pure Nash equilibria are usually search for a function f(x, y). Yet, the problems arising from machine learning are usually nonconvex in x and nonconcave in y, in which case pure Nash equilibrium may not exist. And even if it exists, there is no guarantee for any optimization algorithm to find it efficiently. This difficulty is reflected in practice, that compared with simple minimization, machine learning applications involving Nash equilibria usually have more complicated behaviors and more subtle dependence on hyper-parameters. For example, stable and efficient training of GANs requires a number of carefully designed tricks (Gao et al., 2018).
On the other hand, the mixed Nash equilibrium (MNE) is known to exist in much more general settings, e.g. when the strategy spaces are compact and the payoff function is continuous (Glicksberg, 1952). In the mixed Nash equilibrium problem, instead of taking “pure strategies” x and y, two ”mixed strategies” for x and y, in the form of probability distributions, are considered, resulting in the following functional, ∫
f(x, y)p(x)q(y)dxdy,
where p and q are density functions of probability distributions of x and y, respectively. Efforts are invested to develop theoretically endorsed algorithms that can efficiently find MNE for high dimensional problems, with applications on the training of mixture of GANs. In Hsieh et al. (2019), a mirror-descent algorithm is proposed and its convergence is proven. In Domingo-Enrich et al. (2020), theoretical analysis and empirical experiments are conducted for a gradient descent-ascent flow under a Wasserstein-Fisher-Rao metric and its particle discretization.
In this paper, we also consider the mixed Nash equilibrium problem, and propose a simple QuasiStatic Wasserstein Gradient Flow (QSWGF) for solving the problem. In our dynamics, we treat q as a component with much faster speed than p, hence is always at equilibrium as p moves. With entropy regularization for both p and q (without requirement on the strength of the regularizations), we prove that the QSWGF converges to the unique mixed Nash equilibrium from any initialization (under mild conditions). Furthermore, we show there is a simple way to discretize the QSWGF, regardless of the complexity of the Wasserstein gradient flow of p induced by the fact that q is always
at equilibrium. Concretely, a partition function related with p appears in the QSWGF dynamics, and we find an efficient way to approximate the partition function. By discretizing the QSWGF, we derive a particle dynamics with an inner-outer structure, named the QuasiStatic Langevin Gradient Descent algorithm (QSLGD). In QSLGD, after each iteration of the outer problem (for the x particles), the inner loop conducts sufficient iterations to bring the y particles to equilibrium. Numerical experiments show the effectiveness of QSLGD on synthetic examples and training mixture of GANs. Our method outperforms the vanilla Langevin gradient descent-ascent method when the entropy regularization is weak.
As a summary, our two major contributions are:
1. We propose the quasistatic Wasserstein gradient flow dynamics for mixed Nash equilibrium problems, and show its convergence to the unique Nash equilibrium under weak assumptions. Our result neither requires the entropy regularization to be sufficiently strong, nor assumes the dynamics to converge a priori.
2. We derive a simple while practical quasistatic Langevin gradient descent algorithm by discretizing the quasistatic Wasserstein gradient flow, by finding an efficient way to approximate the partition function appearing in the dynamics of p. The proposed algorithm is applied on several problems including training mixtures of GANs.

heading: 2 RELATED WORK
text: The mixed Nash equilibrium problem has a long history, with the proof of its existence dates back to Morgenstern & Von Neumann (1953). It draws new attention in recent years, especially in the machine learning community, due to the development of GANs (Goodfellow et al., 2014a) and adversarial training (Goodfellow et al., 2014b). Training mixture of GANs is already discussed in paper (Goodfellow et al., 2014a). Some numerical experiments were conducted in (Arora et al., 2017). In Grnarova et al. (2017), the authors proposed an online learning approach for training mixture of GANs, and proved its effectiveness for semi-shallow GANs (GANs whose discriminator is a shallow neural network). Yet, rigorous theoretical treatment to an algorithm started from (Hsieh et al., 2019), in which a mirror descent method was studied and proven to converge. The implementation of the mirror descent method involves big computational cost that asks for heuristics to alleviate. Later, (Domingo-Enrich et al., 2020) studied more efficient algorithms under a mixture of Wasserstein and Fisher-Rao metrics. Theoretically, the time average of the dynamics’ trajectories is shown to converge to the mixed Nash equilibrium. As a comparison, in this work we show the global convergence of the quasistatic Wasserstein gradient flow without the need of taking time average. Meanwhile, the Wasserstein nature of our dynamics makes it easy to implement as well.
The Wasserstein gradient flow in the density space has been explored in previous works. For example, (Wang & Li, 2019) studied the Nesterov’s accelerated gradient flows for probability distributions under the Wasserstein metric, and (Arbel et al., 2019) studied practical implementations of the natural gradient method for the Wasserstein metric. Both works focus on minimization problems instead of min-max problems considered in this work. A more related work is Lin et al. (2021b), where a natural gradient based algorithm is proposed for training GANs. Yet, the method still optimizes one generator and one discriminator, searching for pure Nash equilibrium. Another work that derives algorithms for GANs from a Wasserstein perspective is (Lin et al., 2021a).
Another volume of works that studies the Wasserstein gradient flow in the machine learning context is the mean-field analysis of neural networks. This line of works started from two-layer neural networks (Mei et al., 2018; Rotskoff & Vanden-Eijnden, 2018; Chizat & Bach, 2018; Sirignano & Spiliopoulos, 2020), to deep fully-connected networks (Araújo et al., 2019; Sirignano & Spiliopoulos, 2021; Nguyen, 2019; Wojtowytsch et al., 2020), and residual networks (Lu et al., 2020; E et al., 2020). The mean-field formulations treat parameters as probability distributions, and the training dynamics are usually the gradient flow under Wasserstein metric. Attempts to prove convergence of the dynamics to global minima are made (Mei et al., 2018; Chizat & Bach, 2018; Rotskoff et al., 2019), though in the case without entropy regularization a convergence assumption should usually be made a priori.

heading: 3 THE QUASISTATIC DYNAMICS
text: We consider the entropy regularized mixed Nash equilibrium problem, which in our case is equivalent with solving the following minimax problem:
min p∈P(Ω) max q∈P(Ω) ∫ Ω×Ω K(x, y)p(x)q(y)dxdy + β−1 ∫ Ω p log pdx− β−1 ∫ Ω q log qdy. (1)
In (1), Ω is a compact Riemannian manifold without boundary, and P(Ω) is the set of probability distributions on Ω. Since Ω is compact, any probability distribution in P(Ω) naturally has finite moments. Let E(p, q) = ∫ Ω×ΩK(x, y)p(dx)q(dy), and S(p) = ∫ Ω p log pdx and S(q) = ∫ Ω q log qdy be the (negative) entropy of p and q, respectively. Then, the minimax problem (1) can be written in short as
min p∈P(Ω) max q∈P(Ω)
E(p, q) + β−1S(p)− β−1S(q). (2)
Remark 1. Strictly speaking, in (1) we should distinguish probability distributions and their density function (if exist), and the entropy should also be defined using the Radon-Nikodym derivative with canonical measure. In this paper, since p and q indeed have density functions because of the entropy regularization, we shall abuse the notation by using p and q to represent both probability distributions and their density functions.
The entropy regularizations in (1) and (2) make the problem strongly convex in p and strongly concave in q. Hence, there exists a unique Nash equilibrium for the problem. Such results are shown for example by the following theorem from (Domingo-Enrich et al., 2020). Theorem 1. (Theorem 4 of (Domingo-Enrich et al., 2020)) Assume Ω is a compact Polish metric space equipped with canonical Borel measure, and that K is a continuous function on Ω × Ω. Then, problem (2) has a unique Nash equilibrium given by the solution of the following fixed-point problem:
p(x) = 1 Zp exp(−βU(x, q)), q(x) = 1 Zq exp(βV (y, p)), (3)
where Zp and Zq are normalization constants to make sure p and q are probability distributions, and U and V are defined as
U(x, q) = δE(p, q)
δp (x) = ∫ Ω K(x, y)q(y)dy, V (y, p) = δE(p, q) δq (y) = ∫ Ω K(x, y)p(x)dx.
Considering the efficiency in high-dimensional cases, a natural dynamics of interest to find the Nash equilibrium for (2) is the gradient descent-ascent flow under the Wasserstein metric,
∂tpt = ∇ · ( pt∇(U(x, qt) + β−1 log pt) ) ,
∂tqt = ∇ · ( qt∇(−V (y, pt) + β−1 log qt) ) , (4)
because it can be easily discretized into a Langevin gradient descent-ascent method by treating the PDEs as Fokker-Planck equations of SDEs. When β−1 is sufficiently large, (4) can be proven to converge linearly to the unique MNE of (2) (Eberle et al., 2019). However, when β−1 is small, whether (4) converges remains open. This hinders the application of (4) because in practice the entropy terms are usually used as regularization and are kept small. (We realize that it is proven in Domingo-Enrich & Bruna (2022) when our work is under review.)
In (4), the dynamics of p and q have the same speed. In this work, instead, we study a quasistatic Wasserstein gradient descent dynamics, which can be understood as a limiting dynamics when the speed of q becomes faster and faster compared with that of p. In this case, at any time t, we assume qt reaches at the equilibrium of the maximizing problem instantaneously by fixing p = pt in (2). That is to say, at any time t, qt is determined by
qt = q[pt] := arg max q∈P(Ω)
E(pt, q)− β−1S(q). (5)
On the other hand, pt follows the Wasserstein gradient descent flow with qt = q[pt] at the equilibrium: ∂tpt = ∇ · ( pt∇ ( δ(E(pt, q[pt])− β−1S(q[pt]))
δpt + β−1 log pt
)) . (6)
The following theorem shows qt = q[pt] can be explicitly written as a Gibbs distribution depending on pt, and thus the free energy in (6) can be simplified to depend on a partition function related with pt. Theorem 2. Assume K is continuous on the compact set Ω and β > 0. Then, for fixed pt the maximization problem (5) has a unique solution
q[pt](y) := 1
Zq(pt) exp(βV (y, pt)), (7)
where Zq(p) is a normalization factor, Zq(p) := ∫
exp(βV (y, p))dy. Moreover, the dynamics (6) for pt can be written as
∂tpt = ∇ · ( pt∇ ( δβ−1 logZq(pt)
δpt + β−1 log pt
)) . (8)
Let Fp,β(p) := β−1 logZq(p) +β−1S(p). By Theorem 2, the dynamics (8) of pt is the Wasserstein gradient descent flow for minimizing Fp,β(p). By the Proposition 3 below, Fp,β is strongly convex with respect to p. Therefore, it is possible to prove global convergence for the dynamics (8), and thus the convergence for the quasistatic Wasserstein gradient flow for the minimax problem (2). Proposition 3. For any probability distributions p1, p2 in P(Ω), and any λ ∈ [0, 1], we have
Fp,β(λp1 + (1− λ)p2) < λFp,β(p1) + (1− λ)Fp,β(p2).
In practice the partition function logZq(pt) in (8) seems hard to approximate, especially when Ω is in high dimensional spaces. However, we show in the following proposition that the variation of the partition function with respect to pt can be written as a simple form involving qt. This property will be used to derive a particle method in Section 5 Proposition 4. For any p ∈ P(Ω), we have
δβ−1 logZq(p)
δp = U(·, q[p]), (9)
where q[p] is defined in (7). Therefore, the dynamics (8) is equivalent with ∂tpt = ∇ · ( pt∇ ( U(x, q[pt]) + β −1 log pt )) . (10)

heading: 4 CONVERGENCE ANALYSIS
text: In this section, we analyze the convergence of the quasistatic dynamics (7), (8). First, we make the following assumptions on K. Assumption 1. Assume K ∈ C∞(Ω×Ω), which means K has continuous derivatives of any order (with respect to both x and y).
Since Ω is compact, assumption 1 implies boundedness and Lipschitz continuity of any derivatives of K.
Now, we state our main theorem, which shows the convergence of QSWGF to the Nash equilibrium. Theorem 5. (main theorem) Assume Assumption 1 holds for K. Then, starting from any initial p0, q0 ∈ P(Ω), the dynamics (7), (8) has a unique solution (pt, qt)t≥0, and the solution converges weakly to the unique Nash equilibrium of (2), (p∗, q∗), which satisfies the fixed point problem (3).
Theorem 5 guarantees convergence of the quasistatic Wasserstein gradient flow for any β, giving theoretical endorsement to the discretized algorithm that we will introduce in the next section. Note that the initialization q0 in the theorem is not important, because we assume q achieves equilibrium immediately after the initialization. Remark 2. The assumption on K’s smoothness can be made weaker. For example, during the proof, up to 4-th order derivatives of K is enough to give sufficient regularity to the solution of the dynamics. We make the strong assumption partly to prevent tedious technical analysis so as to focus on the idea and insights.
Proof sketch We provide some main steps and ideas of the proof of the main theorem in this section. The detailed proof is put in the appendix.
By the last section, since qt is always at equilibrium, we only need to considering a Wasserstein gradient descent flow for Fp,β(p). Therefore, we can build our analysis based on the theories in (Mei et al., 2018) and (Jordan et al., 1998). However, compared with the analysis therein, our theory deals with a new energy term—β−1 logZq(p), which has not been studied by previous works. From now on, let Ep,β(p) = β−1 logZq(p), and Ψ(·, p) = δEp,β(p)δp . By simple calculation we have
Ψ(x, p) = U(x, q[p]) = 1
Zq(p) ∫ Ω K(x, y) exp (∫ Ω βK(x, y)p(x)dx ) dy. (11)
First, we study the free energy Fp,β(p), and show that it has a unique minimizer which satisfies a fixed point condition. This is the result of the convexity of Fp,β . We have the following lemma. Lemma 1. Assume Assumption 1 holds for K. Then, Fp,β has a unique minimizer p∗ that satisfies
Fp,β(p ∗) = inf
p∈P(Ω) Fp,β(p
∗).
Moreover, p∗ is the unique solution of the following fixed point problem,
p∗ = 1
Z exp (−βΨ(x, p∗)), (12)
where Z is the normalization factor.
Next, we want to show that any trajectory given by dynamics (10) will converge to the unique minimizer of Fp,β . To achieve this, we first study the existence, uniqueness, and regularity of the solution to (8), i.e. the trajectory indeed exists and is well behaved. Related results are given by the following lemma. Lemma 2. Assume Assumption 1 holds for K. Then, starting from any initial p0 ∈ P(Ω), the weak solution (pt)t≥0 to (8) exists and is unique. Moreover, (pt) is smooth on (0,∞)× Ω.
The proof of Lemma 2 is based on Proposition 5.1 of (Jordan et al., 1998). Especially, the existence part is proven using the JKO scheme proposed in (Jordan et al., 1998). We consider a sequence of probability distributions given by the following discrete iteration schemes with time step h,
ph0 = p0, p h k = arg min
p∈P(Ω)
{ 1
2 W 22 (p, p h k−1) + hFp,β(p)
} , k > 0,
where W2(p, q) means the 2-Wasserstein distance between probability distributions p and q. Let (pht )t≥0 be the piecewise constant interpolations of (p h k)k≥0 on time. We show (p h t ) converges weakly (after taking a subsequence) to a weak solution of (8) as h tends to 0. Details are given in the appendix.
Finally, noting that Fp,β is a Lyapunov function of the dynamics (8), we have the following lemma showing the convergence of (pt)t≥0 to the solution of the Boltzmann fixed point problem (12). This finishes the proof of the main theorem. Lemma 3. Let (pt)t≥0 be the solution of (8) from any initial p0 ∈ P(Ω). Let p∗ be the unique minimizer of Fp,β given by (12). Then, pt converges to p∗ weakly as t→∞.
As a byproduct, since our convergence results does not impose requirement on β, if one is interested in the minimax problem without entropy regularization,
min p∈P(Ω) max q∈P(Ω) E(p, q), (13)
then, Theorem 5 in (Domingo-Enrich et al., 2020) ensures that the quasistatic dynamics converges to approximate Nash equilibrium of (13) as long as β−1 is small enough. Specifically, a pair of probability distributions (p, q) is called -Nash equilibrium of (13) if
sup q′∈P(Ω) E(p, q′)− inf p′∈P(Ω) E(p′, q) ≤ .
Then, we have the following theorem as a direct results of Theorem 5 in (Domingo-Enrich et al., 2020):
Theorem 6. Let CK be the bound of K that satisfies |K(x, y)| ≤ CK for any x, y ∈ Ω, and let Lip(K) be the Lipschitz constant of K. For any > 0, let δ = /(2Lip(K)), and let Vδ be the volume of a ball with radius δ in Ω. Then, as long as
β > 4 log
( 2(1− Vδ)
Vδ
( 4CK − 1 )) ,
there exists T > 0 which depends on , such that for any t > T , the solution pt, qt of the dynamics (7) and (8) at t satisfies
sup q′∈P(Ω)
E(pt, q ′)− inf p′∈P(Ω) E(p′, qt) ≤ .

heading: 5 THE QUASISTATIC LANGEVIN GRADIENT DESCENT-ASCENT METHOD
text: It is well known that PDEs with the form ∂tp(t, x) = ∇ · (p(t, x)µ(t, x)) + λ∆p(t, x) are Fokker-Planck equations for SDEs dXt = −µ(t,Xt)dt + √
2λdWt, and the solution for the PDE characterizes the law of Xt—the solution of the SDE—at any time. This result connects the Wasserstein gradient flow with SDE, and gives a natural particle discretization to approximate the continuous Wasserstein gradient flow. For example, the Wasserstein gradient descent-ascent flow dynamics (4) is the Fokker-Planck equation of the SDEs
dXt = −∇xU(Xt, qt)dt+ √ 2β−1dWt
dYt = ∇yV (Yt, pt)dt+ √ 2β−1dW ′t ,
where pt and qt are the laws of Xt and Yt, respectively, and Wt and W ′t are two Brownian motions. Note that we have
∇xU(x, q) = ∫
Ω
∇xK(x, y)q(y)dy, ∇yV (y, p) = ∫
Ω
∇yK(x, y)p(x)dx.
Therefore, i.i.d. picking X(i)0 ∼ p0 and Y (i) 0 ∼ q0 for i = 1, 2, ..., n, the particle update scheme, named Langevin Gradient Descent-Ascent (LGDA),
X (i) k+1 = X (i) k −
h
n n∑ j=1 ∇xK(X(i)k , Y (j) k ) + √ 2hβ−1ξ (i) k ,
Y (i) k+1 = Y (i) k +
h
n n∑ j=1 ∇yK(X(j)k , Y (i) k ) + √ 2hβ−1ζ (i) k , (14)
approximately solves the SDEs, and thus the empirical distributions of X(i)k and Y (i) k approximate the solutions of (4) when n is large. Here, ξ(i)k and ζ (i) k are i.i.d. samples from the standard Gaussian.
Quasistatic Langevin gradient descent method Similarly, the dynamics (8) for p is the FokkerPlanck equation for the SDE
dXt = −∇Ψ(x, pt)dt+ √
2β−1dWt, (15) where pt is the law of Xt. By proposition 4 we have Ψ(x, pt) = U(x, q[pt]). Hence, (15) can be written as
dXt = −∇xU(Xt, q[pt])dt+ √
2β−1dWt, (16) with q[pt] at the equilibrium of the maximization problem (5), which can be attained by solving the SDE
dYt = ∇yV (Yt, pt)dt+ √
2β−1dW ′t (17) for sufficiently long time. This motivates us to design a quasistatic particle method as a discretization for the quasistatic Wasserstein gradient flow. Specifically, the method consists of an inner loop and an outer loop. The method starts from some particles X(i)0 and Y (i) 0 , i = 1, 2, ..., n, sampled i.i.d. from p0 and q0, respectively. Then, at the k-th step, the inner loop conducts enough iterations on the Y particles to solve (17) with pt fixed (i.e. with the X particles fixed), which drives the empirical distribution of {Y (i)k }ni=1 near equilibrium before each update of the outer loop. Next, the outer loop updates X(i)k according the SDE (16). The algorithm is summarized in Algorithm 1.
Algorithm 1: Quasistatic Langevin gradient descent method (QSLGD) input : nx, ny, k0, k1, k2, T ∈ N+, hx, hy > 0, p0, q0 ∈ P(Ω) output: Final particles (X(i)T , Y (i) T ) n i=1
1 Sample (X(i)0 ) nx i=1 i.i.d. from p0, and (Y (i) 0 ) ny i=1 i.i.d. from q0; 2 Y (i) 0,0 ← Y (i) 0 for i = 1, 2, ..., ny; 3 for s← 1 to k0 do 4 Y
(i) 0,s ← Y (i) 0,s−1 + hy nx nx∑ j=1 ∇yK(X(j)0 , Y (i) 0,s−1) + √ 2hyβ−1ξ ; /* remark 3.1 */
5 end 6 Y
(i) 0 ← Y (i) 0,k0
for i = 1, 2, ..., ny; 7 for t← 1 to T do 8 Y
(i) t−1,0 ← Y (i) t−1 for i = 1, 2, ..., ny;
9 for s← 1 to k1 do 10 Y
(i) t−1,s ← Y (i) t−1,s−1 + hy nx nx∑ j=1 ∇yK(X(j)t−1, Y (i) t−1,s−1) + √ 2hyβ−1ξ;
11 end 12 for s← 1 to k2 do 13 Y
(i) t−1,s+k1 ← Y (i) t−1,s+k1−1 + hy nx nx∑ j=1 ∇yK(X(j)t−1, Y (i) t−1,s+k1−1) + √ 2hyβ−1ξ;
14 Ŷ ((s−1)ny+i) t−1 ← Y (i) t−1,s+k1 , for i = 1, 2, ..., ny; 15 end
16 X (i) t ← X (i) t−1 + hx k2ny k2ny∑ j=1 ∇yK(X(i)t−1, Ŷ (j) t−1) + √ 2hxβ−1ξ ; /* remark 3.2 */ 17 Y (i) t ← Y (i) t−1,k1+k2 , for i = 1, 2, ..., ny; 18 end
Remark 3. Generally speaking, Algorithm 1 consists of two nested loops. The inner loop solves Y particles to equilibrium in each step of the outer loop, while the outer loop makes one iteration every time, using the equilibrium Y particles. In the following are some additional explanation for the Algorithm:
1 line 4: at the beginning of the algorithm, we conduct k0 additional inner iterations for Y , where k0 may be a large number. This is because at the beginning the Y particles are far from equilibrium. In later outer iterations, since each time the X particles only move for a small distance, the Y particles are close to the equilibrium. Therefore, k1 and k2 need not to be large.
2 line 17: In each inner loop, we conduct k1 + k2 inner iterations for the Y particles, and collect those from the last k2 iterations. We use these k2n particles in the update of X particles to approximate the distribution q[p]. We assume during the last k2 inner iterations the Y particles are at equilibrium. One can take k2 to be 1 if ny is large enough, while taking large k2 allows smaller number of Y particles.

heading: 5.1 EXAMPLES
text: In this section, we apply the quasistatic Langevin gradient descent method to several problems.
1-dimensional game on torus We first consider a problem with x and y on the 1-dimensional torus. Specifically, we consider
K(x, y) = sin(2πx) sin(2πy),
where x, y ∈ R/Z. It is easy to show that, with this K and a positive β, at the Nash equilibrium of the problem (1) p and q are both uniform distributions. We take initial distributions p0 and q0
to be the uniform distribution on [0, 1/4]. Figure 1 shows the comparison of the quasistatic particle method with LGDA for different β, step length, and number of particles. In the experiments, all quasistatic methods take k0 = 1000 and k2 = 1, with different k1 shown in the legends. For each experiment, we conduct 300000, 150000, 60000, 30000 outer iterations for LGDA, QS2, QS5, and QS10, respectively. We take different different numbers of iterations for different methods in the consideration of different number of inner iterations. The error is then computed after the last iteration, measured by the KL divergence of the empirical distribution given by particles and the uniform distribution (both in the forms of histograms with 10 equi-length bins). Each point in the figures is an average of 5 experiments.
Seen from the left figure, the QSLGD has comparable performance than LGDA when β is small, in which case diffusion dominates the dynamics, while it performs much better than LGDA when β is large. We can also see better tolerance to large β when more inner iterations are conducted. This shows the advantage of the QSLGD over LGDA when the regularization stength is weak. The middle figures shows slightly better performance of the QSLGD when the step length η (both ηx and ηy) is small. However, when η is big, LGDA tends to give smaller error. The results may be caused by the instability of the inner loop when η is big. It also guides us to pick small step length when applying the proposed method. Finally, the right figure compares the influence of the number of particles when β = 100 and η = 0.01, in which case the two methods perform similarly. We can see that the errors for both methods scale in a 1/n rate as the number of particles n changes.
Polynomial games on spheres In the second example, we consider a polynomial games on sphere similar to that studied in (Domingo-Enrich et al., 2020),
K(x, y) = xTA0x+ x TA1y + y TA2y + y TA3(x 2), (18)
where x, y ∈ Sd−1 and (x2) is the element-wise square of x. In this problem, we consider the Nash equilibrium of minp maxq E(p, q). Hence, we take big β (small β−1) and compare the Nikaido and Isoda (NI) error of the solutions found by different methods (Nikaidô & Isoda, 1955). The NI error is defined by
NI(p, q) := sup q′∈P(Ω)
E(pt, q ′)− inf p′∈P(Ω) E(p′, qt),
which is also used in Theorem 6. The left panel of Figure 2 shows the NI errors of the solutions found by different methods with different dimensions, we see comparable performance of the QSLGD with LGDA.
GANs Finally, we test our methods on the training of GANs. We train GANs to learn Gaussian mixtures. The results after training are shown in the middle and right panels of Figure 2, where Gaussian mixtures with 4 and 8 modes are learned, respectively. We train GANs with 5 generators and 5 discriminators, and take k0 = 100, k1 = 5, k2 = 1. The results show that the mixture of GANs trained by QSLGD can learn Gaussian mixtures successfully.
In the right panel of Figure 2, we show the results of learning high dimensional Gaussian mixtures. In the d-dimensional experiment, the Gaussian mixture has d modes centered at e1, e2, ..., ed with
standard deviation 0.1. Here, ei is the i-th unit vector in the standard basis of Rd. Model and algorithm with same hyper-parameters as above are used. In the figure, we measure the average squared distance of the generated data to the closest mode center along the training process. The figure shows that the average squared distance can be reduced to 0.3 − 0.5 after 10000 iterations. While the ideal value is 0.1, the current results still show that the learnt distribution concentrates at the mode centers. Better results may be obtained after longer training or careful hyper-parameter tuning.

heading: 6 DISCUSSION
text: In this paper, we study the quasistatic Wasserstein gradient flow for the mixed Nash equilibrium problem. We theoretically show the convergence of the continuous dynamics to the unique Nash equilibrium. Then, a quasistatic particle method is proposed by discretizing the continuous dynamics. The particle method consists of two nested loops, and conduct sufficient inner loop in each step of the outer loop. Numerical experiments show the effectiveness of the method. Comparison with LGDA shows the proposed method has advantage over LGDA when β is large (which is usually the case of interest), and performs as good as LGDA in most other cases.
Theoretical extensions are possible. For example, strong convergence results may be established by similar approaches taken in (Feng & Li, 2020). We leave this as future work.
In practice, the idea of nested loops is not new for minimax optimization problems. It is already discussed and utilized in the earliest works for GANs (Goodfellow et al., 2014a), and Wasserstein GANs (Arjovsky et al., 2017). In those works, the discriminator is updated for several steps each time the generator is updated. Our work is different from these works because we consider mixed Nash equilibrium and hence our method is particle based, while their method searches for pure Nash equilibrium.
Finally, though particle methods finding mixed Nash equilibria have stronger theoretical guarantees, applying these methods to the training of GANs faces the problem of computational cost. With both the generator and discriminator being large neural networks, training mixture of GANs with many generators and discriminators imposes formidable computational cost. Developing more efficient particle methods for GANs is an important future work.

heading: A PROOFS FOR SECTION 3

heading: A.1 PROOF OF THEOREM 2
text: Note that the free energy in (5) can be written as∫ Ω V (y, pt)q(y)dy − β−1S(q), (19)
in which the first term is linear with respect to q. Hence, a calculation with Lagrange multiplier shows (19) has a unique minimizer q[pt] with the form of a Gibbs distribution (e.g. see Chapter 4 of (Mezard & Montanari, 2009)):
q[pt](y) = 1
Zq(pt) exp(βV (y, pt)). (20)
Next, we consider the free energy for pt when q is at the equilibrium. By (19) we have
E(pt, q[pt]) = ∫ Ω V (y, pt)q[pt](y)dy
= 1
Zq(pt) ∫ Ω V (y, pt) exp(βV (y, pt))dy. (21)
On the other hand, we have β−1S(q[pt]) = β −1 ∫
Ω
1
Zq(pt) exp(βV (y, pt)) log
( 1
Zq(pt) exp(βV (y, pt))
) dy
= β−1 ∫
Ω
1
Zq(pt) exp(βV (y, pt)) (βV (y, pt)− logZq(pt)) dy
= 1
Zq(pt) ∫ Ω V (y, pt) exp(βV (y, pt))dy − β−1 logZq(pt). (22)
Combining (21) and (22), we obtain E(pt, q[pt]) + β
−1S(pt)− β−1S(q[pt]) = β−1 logZq(pt) + β−1S(pt). Therefore, the dynamics of pt is the Wasserstein gradient descent flow minimizing the free energy β−1 logZq(pt) + β −1S(pt), given by
∂tpt = ∇ · ( pt∇ ( δβ−1 logZq(pt)
δpt + β−1 log pt
)) .
This finishes the proof.

heading: A.2 PROOF OF PROPOSITION 3
text: Since S(p) is strongly convex, it suffices to show logZq(p) is convex. Recall that
logZq(p) = log (∫ Ω exp(βV (y, p))dy ) .
Note that V (·, p) = ∫ K(x, ·)p(x)dx is linear with respect to p, we have V (·, λp1 + (1− λ)p2) = λV (·, p1) + (1− λ)V (·, p2).
Hence, logZq(λp1 + (1− λ)p2) = log (∫
Ω
exp(βV (y, λp1 + (1− λ)p2))dy )
= log (∫ Ω exp(βλV (y, p1)) · exp(β(1− λ)V (y, p2)) )
≤ log ((∫ Ω exp(βV (y, p1)) )λ(∫ Ω exp(βV (y, p2)) )1−λ) = λ logZq(p1) + (1− λ) logZq(p2). (23)
The second last line is given by the Hölder inequality.

heading: A.3 PROOF OF PROPOSITION 4
text: The proposition follows from the following derivations.
δβ−1 logZq(p)
δp = β−1
1
Zq(p)
δZq(p)
δp
= β−1 1
Zq(p) ∫ Ω exp(βV (y, p))βK(x, y)dy
= ∫ Ω K(x, y) exp(βV (y, p)) Zq(p) dy
= ∫ Ω K(x, y)q[p](y)dy
= U(x, q[p]). (24)

heading: B PROOF OF THEOREM 5
text: In this section, we prove our main theorem. The proof will follow the sketch given in Section 4. Given the assumptions on K and the conclusions of Lemma 1 and Lemma 2, Lemma 3 is a direct result of Lemma 10.12 in (Mei et al., 2018), which we will ignore the proof. In the following, we show Lemma 1 and Lemma 2. Some techniques in the proof come from (Jordan et al., 1998) and (Mei et al., 2018).

heading: B.1 PROOF OF LEMMA 1
text: Our proof follows the proof of Proposition 4.1 in (Jordan et al., 1998). First, we show the existence of the minimizer for Fp,β . To see this, note that K is bounded on Ω. Assume CK is a constant such that |K(x, y)| ≤ CK for any x, y ∈ Ω. Then, we have
Ep,β(p) = ∫ Ω U(x, q[p])p(x)dx = ∫ Ω×Ω K(x, y)q[p](y)p(x)dxdy ≥ −CK ,
and
S(p) = ∫ Ω p(x) log p(x)dx ≥ ∫ Ω −1 e dx = −1 e .
This means Fp,β(p) is lower bounded, i.e. infp Fp,β(p) > −∞. Hence, we can find a sequence (pk) ∞ k=1 such that
lim k→∞ Fp,β(pk) = inf p Fp,β(p). Similar to (Jordan et al., 1998), we can show boundedness of { ∫
max{pk log pk, 0}dx} and { ∫ p2kdx}, which implies that (pk) is uniformly integrable, and thus there exists a weakly convergent subsequence of (pk).
Without loss of generality, assume pk ⇀ p∗ in L1(Ω). Then we need to show p∗ is a minimizer of Fp,β . By (Jordan et al., 1998), the entropy term satisfies
S(p∗) ≤ lim inf k→∞ S(pk).
Hence, the conclusion follows if Ep,β is continuous in the weak topology. To show this, first note that for any p ∈ P(Ω), we have ∫
Ω
exp(βV (y, p))dy ≥ e−βCK .
Because the function log(x) is 1/c-Lipschitz for x ∈ [c,∞], for any pk we have∣∣β−1 logZq(pk)− β−1 logZq(p∗)∣∣ ≤ β−1eβCK ∣∣∣∣∫ Ω ( eβV (y,pk) − eβV (y,p ∗) ) dy ∣∣∣∣ .
By similar boundedness argument, we have∣∣∣∣∫ Ω ( eβV (y,pk) − eβV (y,p ∗) ) dy ∣∣∣∣ ≤ ∫ Ω ∣∣∣eβV (y,pk) − eβV (y,p∗)∣∣∣ dy ≤ eβCK
∫ Ω β |V (y, pk)− V (y, p∗)| dy
≤ βeβCK ∫
Ω ∣∣∣∣∫ Ω K(x, y)(pk(x)− p∗(x))dx ∣∣∣∣ dy
Totally we have∣∣β−1 logZq(pk)− β−1 logZq(p∗)∣∣ ≤ e2βCK ∫ Ω ∣∣∣∣∫ Ω K(x, y)(pk(x)− p∗(x))dx ∣∣∣∣ dy.
Since K is bounded and Lipschitz, it is easy to show that
lim k→∞ ∫ Ω ∣∣∣∣∫ Ω K(x, y)(pk(x)− p∗(x))dx ∣∣∣∣ dy = 0.
Therefore, we have Fp,β(p
∗) ≤ lim inf k→∞ Fp,β(pk) = inf p Fp,β(p),
and thus Fp,β(p∗) = infp Fp,β(p).
Next, we show p∗ satisfies the fixed point condition
p∗ = 1
Z exp (−βΨ(x, p∗)). (25)
This follows the proof of Lemma 10.3 in (Mei et al., 2018), by first showing p∗ has full support on Ω, and then showing
Ψ(x, p∗) + β−1 log p∗(x)
is a constant.
Finally, we show p∗ is unique following Lemma 10.4 of (Mei et al., 2018). Specifically, we show the Boltzmann fixed point problem (25) only has one solution by the convexity of Fp,β . Assume (25) has two different solutions p1 and p2, i.e.
p1 = 1
Z(p1) exp (−βΨ(x, p1)), p2 =
1
Z(p2) exp (−βΨ(x, p2)).
Then, we have logZ(p1) = −βΨ(x, p1)− log p1(x), logZ(p2) = −βΨ(x, p2)− log p2(x).
Taking difference of the above two equations, and integrating with p1 − p2, we have
0 = ∫ Ω (logZ(p1)− logZ(p2))(p1(x)− p2(x))dx
= −β ∫
Ω
(Ψ(x, p1)−Ψ(x, p2))(p1(x)− p2(x))dx− ∫
Ω
(log p1(x)− log p2(x))(p1(x)− p2(x))dx.
(26) By the monotonicity of log x, the second term of (26) is non-negative, and takes zero only if p1 = p2. For the first term, recall that we have proven in Proposition 3 that Ep,β is convex, we then have
Ep,β(p1) ≥ Ep,β(p2) + ∫
Ω
Ψ(x, p2)(p1(x)− p2(x))dx,
and Ep,β(p2) ≥ Ep,β(p1) + ∫ Ω Ψ(x, p1)(p2(x)− p1(x))dx. Taking difference of the two equations gives∫ Ω (Ψ(x, p1)−Ψ(x, p2))(p1(x)− p2(x))dx ≥ 0. Therefore, (26) holds if and only if p1 = p2. This finishes the proof of uniqueness of p∗, and also completes the proof of Lemma 1.

heading: B.2 PROOF OF LEMMA 2
text: We show the existence of the weak solution using the JKO scheme used in (Jordan et al., 1998). Let dist(·, ·) be the distance metric on Ω. Then, the 2-Wasserstein distance is on P(Ω) is defined as
W2(p1, p2) :=
( inf
γ∈Γ(p1,p2) ∫ Ω×Ω dist(x1, x2)2dγ(x1, x2) )1/2 ,
where Γ(p1, p2) contains all couplings of p1 and p2, i.e. probability distributions on Ω×Ω with first and second marginals being p1 and p2, respectively. Then, for any h > 0, we consider the sequence of probability distributions (phk) ∞ k=0 obtained by the following iteration scheme:
ph0 = p0, p h k = arg min
p∈P(Ω)
{ 1
2 W 22 (p, p h k−1) + hFp,β(p)
} , k > 0, (27)
By similar argument of Proposition 4.1 in (Jordan et al., 1998), each minimization problem in (27) has a unique solution. Hence, (phk) ∞ k=0 is uniquely defined. Let p h t be the piecewise constant interpolation of (phk) on t, i.e. pht = p h k , for t ∈ [kh, (k + 1)h),
for k = 0, 1, 2, .... We now show that there exists a subsequence of hn → 0 and a pt such that phnt ⇀ pt on (0, T )×Ω for any T > 0 and pt is a weak solution of (8). This is proven in two steps:
1. The existence of weakly convergence subsequence, and
2. pht approximately satisfies the equation (8).
For the first point, we show uniform integrability by showing∫ Ω ‖x‖2phk(x)dx ≤ C (28) and ∫ Ω max{phk log phk , 0}dx ≤ C (29)
for any h and k ≥ 0, and an absolute constant C. Equation (28) follows directly from the compactness of Ω. For (29), note that for any h and k ≥ 0 we have
1 2 W 22 (p h k , p h k−1) + hFp,β(p h k) ≤ hFp,β(phk−1),
which implies Fp,β(p h k) ≤ Fp,β(phk−1). Therefore, ∫ Ω max{phk log phk , 0}dx ≤ S(phk) + ∫ Ω
∣∣min{phk log phk , 0}∣∣ dx ≤ S(phk) + ∫ Ω 1 e dx
≤ Fp,β(phk)− Ep,β(phk) + ∫
Ω
1 e dx
≤ Fp,β(ph0 )− Ep,β(phk) + ∫
Ω
1 e dx.
Since K is bounded, Ep,β is bounded, and thus the above expression is also bounded, which gives (29). With (28) and (29), there exists pt(x) and a sequence (hn) with hn → 0, such that phnt ⇀ pt in L1((0, T ) × Ω) for any T > 0. Moreover, pt ∈ P(Ω) for almost every T . By changing pt on a zero measure set of t, we can assume pt ∈ P(Ω) for any t ∈ (0,∞). With the same analysis as (Jordan et al., 1998), the weak convergence can happen for any t, i.e. phnt ⇀ pt in L
1(Ω) for any t ∈ (0,∞).
For the second point, similar to (Jordan et al., 1998), consider any vector field ξ ∈ C∞(Ω,Ω) and the corresponding flux Φτ given by
∂τΦτ = ξ(Φτ ), Φ0(x) = x,
and let qτ = Φτ ]phk , then we have
1
τ
(( 1
2 W 22 (p h k−1, qτ ) + hFp,β(qτ )
) − ( 1
2 W 22 (p h k−1, p h k) + hFp,β(p h k)
)) ≥ 0 (30)
for any τ > 0. We need to study the limit when τ → 0+. By the calculation in (Jordan et al., 1998) we have
lim sup τ→0+
1
τ
( 1
2 W 22 (p h k−1, qτ )−
1 2 W 22 (p h k−1, p 2 k)
) ≤ ∫
Ω×Ω (y − x)ξ(y)dγ(x, y), (31)
and d
dτ S(qτ ) ∣∣∣∣ τ=0 = − ∫ Ω phk∇ · ξdx, (32)
where the γ in (31) is the optimal transport between phk−1 and p h k . For the Ep,β term, we have
lim τ→0+
1
τ
( Ep,β(qτ )− Ep,β(phk) ) = lim τ→0+ 1 βτ log
[ ∫ Ω
exp(βV (y, qτ ))dy∫ Ω exp(βV (y, phk))dy ] = lim τ→0+ 1 βτ [ ∫ Ω exp(βV (y, qτ ))dy∫ Ω exp(βV (y, phk))dy − 1 ]
= 1
βZq(phk) lim τ→0+
1
τ (∫ Ω (exp(βV (y, qτ ))− exp(βV (y, phk)))dy )
= 1
βZq(phk) ∫ Ω exp(βV (y, phk))β (∫ Ω ∇xK(x, y) · ξ(x)phk(x)dx ) dy
= ∫ Ω ∇xΨ(x, phk) · ξ(x)phk(x)dx. (33)
Combining the above result with (31) and (32), taking both ξ and −ξ, we get from (30) that∫ Ω×Ω (y − x)ξ(y)dγ(x, y) + h ∫ Ω ∇xΨ(x, phk) · ξ(x)phk(x)dx− h β ∫ Ω phk∇ · ξdx = 0 (34)
for any ξ ∈ C∞(Ω,Ω). Then, following the derivation in (Jordan et al., 1998) (proof of Proposition 5.1), as well as the following control
N∑ k=1 W 22 (p h k−1, p h k) ≤ Ch
for any N that satisfies Nh ≤ T for some fixed T , we can integrate (34) over t by viewing phk as pht at appropriate t, and take the limit hn → 0 and show that pt is a weak solution. During the limit, we need to pay special attention to the second term, i.e. the following limit, which is not dealt with in the reference:
lim n→∞ ∫ T 0 ∫ Ω ∇xΨ(x, phnt ) · ξ(x)p hn t (x)dx = ∫ T 0 ∫ Ω ∇xΨ(x, pt) · ξ(x)pt(x)dx. (35)
We prove this by showing ∇xΨ(x, phnt ) converges to ∇xΨ(x, pt) uniformly. For any fixed t, recall that we have
phnt ⇀ pt.
Therefore, for any y ∈ Ω, we have∫ Ω K(x, y)phnt (x)dx→ ∫ Ω K(x, y)pt(x)dx.
Note that ∫
Ω K(x, y)p(x)dx is uniformly continuous with respect to y for any p ∈ P(Ω), we can conclude that ∫
Ω K(x, y)phnt (x)dx converges to ∫ Ω K(x, y)pt(x)dx uniformly over y. Hence, we
have Zq(phnt )→ Zq(pt), and q[p hn t ](y)→ q[pt](y) uniformly. This further implies that ∇xΨ(x, phnt ) = ∫
Ω
∇xK(x, y)q[phnt ](y)dy
converges to∇xΨ(x, pt) for all x ∈ Ω. This finishes the proof of (35), and also completes the proof that pt is a weak solution of of equation (8).
Now we have proven the existence of the weak solution. The regularity and uniqueness of the solution follows the same analysis of Proposition 5.1 in (Jordan et al., 1998).

Review summaries:
This paper focuses on the methodology of finding mixed Nash equilibrium. The authors first introduces the corresponding quasi static Wasserstein gradient flow dynamics for solving the problem, which is a limiting dynamics with $q$ seen as infinitely faster than $p$. Then the continuous dynamics is proved to be convergent under mild assumptions. Furthermore, by discretizing the above gradient flow, the authors proposes a practical algorithm for solving min-max optimization problems with GAN as a special example. Experimental results shows the superiority over the vanilla LGDA algorithm especially with large $\beta$’s. 

The paper under review studies the QuasiStatic Wasserstein Gradient Flow and its application to solve two-player zero-sum games. The corresponding Langevin SDE is derived for the Fokker-Planck equation, and the discretized algorithm (GD) is proposed in this new setting.  Example of the min-max problem on 1-d torus and spheres are provided. 

In this paper, the authors study the minimax problem for finding the mixed Nash equilibrium for mean-field two-player zero-sum games. This is an inf-sup problem in the Wasserstein space. In this paper, the authors consider the case where the sup problem is solved and apply the gradient descent flow for the result inf sup functional. Some theoretical analyses are conducted on these dynamics. They show its convergence to the mixed Nash equilibrium under some conditions. Some numerical schemes for the quasistatic Langevin gradient descent method with inner-outer iterations are presented. They show that the method works on training the mixture of GANs. 

This paper studies a minimax problem arising from finding the mixed Nash equilibrium for mean-field two player games. A quasi-static Wasserstein gradient flow is proposed to solve such optimization problems. By discretizing this flow, a simple while practical Langevin gradient descent algorithm is proposed.