1 Introduction↩︎

Zero-sum games provide a fundamental framework for modeling competitive interactions in which the objectives of two players are perfectly opposed. The mathematical formulation is the saddle-point problem, where one player minimizes a payoff function while the other maximizes it. The study of zero-sum games is grounded in the minimax theorem, first established by von Neumann (1928) [1], which guarantees the existence of an equilibrium under mild conditions.

The mathematical framework given by zero-sum games has important applications of modeling adversarial decision-making in economics, machine learning, and control systems. In economics, the pioneering work of von Neumann and Morgenstern [2] established the foundation of modern game theory through the minimax principle, with important applications in market competition, bargaining, auctions, and decision-making under uncertainty, while further developments by Kuhn and Tucker [3], Blackwell and Girshick [4], and Fudenberg and Tirole [5] strengthened its role in economic analysis and strategic behavior. In modern machine learning, Goodfellow et al. [6] introduced Generative Adversarial Networks (GANs), one of the most influential applications of zero-sum games, where a generator and a discriminator are trained through a minimax optimization process; this adversarial formulation has become central to deep generative modeling and has significantly advanced image generation, representation learning, and robust optimization [7]. In control theory and engineering, Isaacs [8] introduced zero-sum differential games as a foundation for robust control. These ideas have major applications in robotics, missile guidance, cyber-physical systems, and autonomous decision systems.

From a computational perspective, solving bilinear zero-sum games reduces to finding a saddle point of a structured objective; first-order methods are widely used due to their simplicity and scalability. Among the classical methods, the Arrow–Hurwicz and Uzawa methods proposed in [9] are the earliest and most fundamental ones, which simultaneously execute gradient descent on the primal minimizing variables and gradient ascent on the dual maximizing variables. We also have Optimistic Gradient Descent (OGD), which improves stability by incorporating predictive correction terms that reduce the cycling behavior common in standard gradient methods, and the extragradient method, which introduces an additional intermediate step to better handle saddle-point structure and is particularly effective for monotone variational inequalities. These two methods are applied to saddle-point problems [10], achieving the convergence rate O(1/T). Follow-the-Regularized-Leader (FTRL), developed from online learning theory [11], updates strategies by balancing cumulative past payoffs with regularization and is closely connected to no-regret learning and convergence in repeated games. More recently, Implicit Gradient Descent (IGD), proposed by Essid and Tabak [12], improves upon explicit methods by incorporating the opponent’s anticipated future response, yielding stronger local convergence and a natural transition toward Newton’s method near saddle points.

Beyond standard first-order methods, Bailey and Piliouras [13] introduced a different perspective by studying continuous-time FTRL dynamics through Hamiltonian systems rather than focusing only on equilibrium computation. They showed that in zero-sum games, continuous-time FTRL satisfies Hamilton’s equations and preserves an invariant energy function, similar to conservative physical systems such as planetary motion. This explains why trajectories often do not converge directly to the Nash equilibrium, but instead move along bounded periodic orbits around it. Motivated by this, we develop a new method by applying the implicit midpoint discretization scheme to the continuous-time FTRL dynamics, which preserves the Hamiltonian structure more accurately, leading to bounded orbits and more stable long-term behavior while maintaining the geometric properties of the original continuous system.

1.1 Our Contributions↩︎

We propose a new implicit update method for bilinear zero-sum games, named Implicit Midpoint Gradient Descent (2 ), which can be viewed as a principled formulation of implicit gradient descent. Our method is designed to exactly preserve the energy structure of the underlying continuous-time dynamics, providing a key explanation for its stability. We show that the resulting dynamics achieve convergence rates matching those of continuous-time Follow-the-Regularized-Leader (FTRL), establishing a tight correspondence between discrete-time updates and their continuous limit. Importantly, our method does not restrict the learning rate and remains stable for arbitrarily large step sizes. We demonstrate that the proposed method outperforms standard explicit gradient-based approaches, particularly in regimes where oscillations or divergence occur, highlighting the benefits of energy-preserving implicit discretizations in saddle-point optimization. Table 1 demonstrates the benefits of 2 compared to other baseline methods.

5pt 2pt

1.2 Outline of Paper↩︎

Section 2 introduces the two-player zero-sum game framework and reviews the associated continuous-time FTRL dynamics, highlighting their Hamiltonian structure. Section 3 presents our proposed discrete-time update rule and establishes its key properties, including bounded orbits and convergence to the Nash equilibrium. Section 4 provides a detailed two-dimensional analysis to illustrate the dynamics and build intuition for the general case. Section 5 reports numerical experiments comparing our method with standard approaches, demonstrating improved stability and performance. Section 6 concludes with a brief discussion on computational cost and implementation considerations.

2 Preliminaries↩︎

2.1 Two-Player Zero-Sum Games↩︎

We consider a two-player zero-sum game with strategy variables \(x_1 \in X_1 \subset \mathbb{R}^{k_1}\) and \(x_2 \in X_2 \subset \mathbb{R}^{k_2}\). The game is defined by the saddle-point problem \[\label{eq:minimax} \max_{x_1 \in X_1} \min_{x_2 \in X_2} u_1(x_1,x_2),\tag{1}\] where player 1’s utility function is \[u_1(x_1,x_2) = x_1^\top A x_2 + b_1^\top x_1 - b_2^\top x_2,\] with \(A \in \mathbb{R}^{k_1 \times k_2}\), \(b_1 \in \mathbb{R}^{k_1}\), and \(b_2 \in \mathbb{R}^{k_2}\) fixed. Since the game is zero-sum, player 2’s utility is \[u_2(x_1,x_2) = - u_1(x_1,x_2) = - x_1^\top A x_2 - b_1^\top x_1 + b_2^\top x_2.\] A pair of strategies \((x_1^\star, x_2^\star) \in X_1 \times X_2\) is called a Nash equilibrium if \[u_1(x_1, x_2^\star) \le u_1(x_1^\star, x_2^\star) \le u_1(x_1^\star, x_2)\] for all \(x_1 \in X_1\) and \(x_2 \in X_2\). In the unconstrained setting (\(X_i=\mathbb{R}^{k_i}\)), a Nash equilibrium satisfies the first-order optimality conditions \[\nabla_{x_1} u_1(x_1^\star,x_2^\star)=0, \qquad \nabla_{x_2} u_2(x_1^\star,x_2^\star)=0.\] As a result, the standard measure for the distance to the set of Nash equilibria is the \(\ell_2\) norm of the gradient. For most of the results in Sections 3–5, we assume that the game admits a Nash equilibrium, a natural condition. The following proposition provides two sufficient conditions under which this assumption is satisfied for the unconstrained bilinear zero-sum game considered in this paper:

In this paper, our main results focus on the unconstrained setting, i.e., \(X_i=\mathbb{R}^{k_i}\), but much of our analysis extends when \(X_i\) is convex and full-dimensional.

2.2 The Role of the Unconstrained Setting↩︎

In the unconstrained setting, Nash equilibria can be directly computed by solving the linear systems \(Ax_2^*+b_1=0\) and \(-A^\top x_1^* +b_2=0\). Thus, our purpose in studying this setting is not to argue that an iterative learning approach is preferable to direct linear algebra for solving a single game. Rather, the unconstrained setting serves as a canonical model for understanding the qualitative behavior of online optimization methods in normal-form (constrained) zero-sum games.

This role is common in the literature; methods such as optimistic gradient descent and alternating gradient descent were first analyzed in unconstrained bilinear games [10], [20], [21], where the dynamics of the methods can be isolated, before extending the results to constrained or normal-form settings [18], [22]. The reason the model remains informative is that constrained zero-sum dynamics often reduce locally to an unconstrained or affine-constrained system once the active set is identified. In particular, for normal-form games, if the support of a Nash equilibrium—the set of strategies with positive weight—is known, then the dynamics restricted to the corresponding face of the simplex evolve on the affine hull of that face. In this reduced subspace, the local behavior is determined by the corresponding unconstrained setting.

As a result, the unconstrained model should be viewed as the basic geometric testbed for local convergence of learning dynamics in the constrained setting.

2.3 FTRL Dynamics with Hamiltonian Structure↩︎

To find a Nash equilibrium for 1 , we study dynamics based on the Follow-the-Regularized-Leader (FTRL) algorithm. Let \(y_i(t)\) denote the cumulative payoff vector of player \(i\) at time \(t\), and let \(h_i : X_i \to \mathbb{R}\) be a strictly convex regularization function. In continuous time, the FTRL dynamics take the form \[y_i(t) = y_i(0) + \int_0^t \eta \, \nabla_{x_i} u_i(x(s)) \, ds, \qquad i = 1,2,\] \[x_i(t) = \arg\max_{x_i \in X_i} \left\{ \langle x_i , y_i(t) \rangle - h_i(x_i) \right\}, \qquad i=1,2.\] where \(\eta > 0\) is the learning rate and \(x(t) = (x_1(t),x_2(t))\).

We also define the time-averaged strategies as \(\bar{x}_i(T) := \frac{1}{T} \int_0^T x_i(t)\,dt\).

For bilinear zero-sum games, the continuous-time FTRL dynamics generate bounded orbits; the trajectories remain confined to a finite region of the state space for all time rather than diverging to infinity. Bailey and Piliouras (2019) [13] find that the dynamics on a two-agent zero-sum game are Hamiltonian with an invariant energy function. Therefore, if the initial strategy is not a Nash equilibrium, the trajectory does not converge to equilibrium but instead stays a fixed positive Bregman distance away and moves around it on a constant-distance orbit.

The notion of a Hamiltonian system originates from classical mechanics and was introduced by William Hamilton to describe the time evolution of physical systems through energy-preserving dynamics.

Knowing that the continuous-time FTRL dynamics have such nice properties, there exists a simple proof for time-averaged strategies convergent to Nash.

We remark it is straightforward to extend the result when \(b_i\neq \vec{0}\) via the variable substitution introduced in Theorem 1. From Lemma 1, we know that continuous-time FTRL is highly robust to learning-rate selection, with the learning rate only affecting the speed of convergence but not stability. Also, the Hamiltonian interpretation of zero-sum game dynamics provides a powerful perspective for algorithm design and analysis. Since continuous-time FTRL can be viewed as a Hamiltonian system, we can leverage techniques developed for Hamiltonian dynamics to simulate or discretize these flows efficiently. For instance, symplectic integrators, which are widely used in Hamiltonian mechanics to preserve geometric structure and energy over long time horizons, can be adapted to design discrete-time FTRL algorithms that inherit the stability and convergence properties of the continuous system. By combining the geometric insights from Hamiltonian dynamics with the structure of FTRL, we can develop novel algorithms that are both robust to step size choices and faithful to the underlying game geometry, potentially outperforming classical gradient-based methods in zero-sum games.

2.4 Gradient Based Methods↩︎

In this section, we study several classical gradient-based methods for online optimization in games. In particular, we view these algorithms as discrete approximations of the continuous FTRL learning dynamics introduced in the previous section. Under this perspective, many characteristic properties of classical methods emerge naturally from the underlying continuous-time learning dynamics they approximate.

Note that in the unconstrained setting \(X_i = \mathbb{R}^{k_i}\), choosing the quadratic regularizer \(h_i(x_i)=\frac{1}{2}\|x_i\|^2\) makes the maximization problem \[x_i(t)=\arg\max_{x_i\in \mathbb{R}^{k_i}}\left\{\langle x_i,y_i(t)\rangle-\frac{1}{2}\|x_i\|^2\right\}\] have the explicit solution \[x_i(t)=y_i(t)=x_i(0) + \int_{0}^t \eta \nabla_{x_i} u_i(x(s))ds, \qquad i = 1,2,\] yielding the continuous-time gradient descent dynamics. In this section, we explore common discrete-time gradient-based methods that can be expressed as an approximation of continuous-time gradient descent.

2.4.0.1 Gradient Descent (GD).

In a two-player game, gradient descent methods update each player’s strategy by moving in the direction that locally improves their own utility. For the minimax optimization problem 1 , the update rule is: \[x_1^{t+1} = x_1^t + \eta \nabla_{x_1} u_1(x_1^t,x_2^t),\] \[x_2^{t+1} = x_2^t + \eta \nabla_{x_2} u_2(x_1^t,x_2^t).\]

This update rule can be obtained by applying explicit Euler’s method to the continuous-time FTRL-induced gradient flow system [13]; it connects gradient-based learning in games to numerical discretizations of continuous-time learning dynamics.

Applying GDA with a constant step size is known to have last-iterate divergence [15], [17]. This divergence is quickly explained when examining GDA as an approximation of continuous-time gradient descent. Since the continuous-time gradient flow admits a Hamiltonian structure, the total system energy is strictly conserved over time. Explicit Euler is well-known to artificially inflate the system’s energy (distance to the set of Nash equilibria) in Hamiltonian systems [23].

This phenomenon also indicates that the rule fails to achieve last-iterate convergence. Nevertheless, \(O(1/\sqrt{T})\) time-average convergence can be recovered by employing decaying learning rates in the constrained setting [16]. More recently, it has been shown that arbitrary fixed learning rates achieve the same results in 2-agent, 2-strategy games [17].

2.4.0.2 Alternating Gradient Descent (AGD).

For alternating gradient descent, agents update their strategies sequentially: \[x_1^{t+1} = x_1^t + \eta \nabla_{x_1} u_1(x_1^t, x_2^t),\] \[x_2^{t+1} = x_2^t + \eta \nabla_{x_2} u_2(x_1^{t+1}, x_2^t).\] This coincides exactly with both the symplectic Euler and Verlet-Störmer discretizations of the continuous-time dynamics [20]. Therefore, AGD can be interpreted as a symplectic integrator for the underlying Hamiltonian dynamics.

Under unconstrained bilinear zero-sum settings,  [24] showed that AGD exhibits bounded orbits by exactly preserving a modified discrete invariant energy function at each iteration, where the discretization errors incurred in opposing quadrants analytically counterbalance each other. This is natural, since AGD is obtained via symplectic integration, a technique known to approximately preserve the energy from the continuous system [23].

While the conservation of energy implies AGD does not achieve last-iterate convergence, fast time-average convergence is successfully recovered. The empirical time-averaged strategies are guaranteed to converge to the exact Nash Equilibrium at an optimal rate of \(O(1/T)\), a property that generalizes to complex multi-agent network zero-sum games [20]. [22] extends this analysis to normal-form games (strategies are probability vectors) when there is a fully-mixed Nash equilibrium. This conservation of energy can be further exploited to characterize the set of Nash equilibria after a finite number of iterations [25].

Furthermore, AGD offers an advantage regarding its learning rate tolerance [20]. It remains structurally stable under a much larger step size than standard anticipatory updates, satisfying the upper bound \(\eta < \frac{2}{\|A\|}\); in contrast, the other standard gradient-based method, optimistic gradient descent, requires much smaller learning rates \(\eta < \frac{1}{2\|A\|}\).

2.4.0.3 Implicit Gradient Descent (IGD).

Implicit gradient methods evaluate the utility gradients at the next iterate rather than the current one. The update rule is given by \[x_1^{t+1} = x_1^t + \eta \nabla_{x_1} u_1(x_1^{t+1},x_2^{t+1}),\] \[x_2^{t+1} = x_2^t + \eta \nabla_{x_2} u_2(x_1^{t+1},x_2^{t+1}).\]

This method can be interpreted as an implicit Euler discretization of the continuous-time gradient flow system. [12] referred to this approach as implicit twisted gradient descent and studied it in the context of minimax optimization problems. The implicit update can be written compactly as \[z^{t+1}=z^t-\eta JF(z^{t+1}),\] with \[z=\begin{pmatrix}x_1\\x_2\end{pmatrix}, \qquad F(z)=\begin{pmatrix}-\nabla_{x_1} u_1\\ -\nabla_{x_2}u_1\end{pmatrix}, \qquad J= \begin{pmatrix} I & 0\\ 0 & -I \end{pmatrix}.\] They performed a first-order Taylor expansion of the future gradient: \[F(z^{t+1}) \approx F(z^t)+H(z^t)(z^{t+1}-z^t),\] where \(H\) denotes the Hessian matrix of the payoff function. This leads to the practical update rule \[z^{t+1} = z^t-\eta (J+\eta H(z^t))^{-1}F(z^t).\] The resulting method interpolates between explicit gradient descent and Newton’s method: for small learning rates it behaves similarly to standard gradient descent/ascent, while for large learning rates it approaches a Newton-type iteration. Furthermore, Essid et al. [12] proved that the algorithm converges locally in a neighborhood of a strict local minimax point. To prevent convergence to critical points of the wrong type, they also proposed an adaptive learning-rate strategy together with acceptance conditions tailored to the minimax structure.

3 Implicit Midpoint Gradient Descent↩︎

We propose a new rule based on the geometric insights from the continuous-time Hamiltonian dynamics. We call the rule Implicit Midpoint Gradient Descent (IMGD) since it’s obtained through the implicit midpoint integration method, which is a one-step, second-order accurate numerical integrator for ordinary differential equations, and the method is a special case of a symplectic Runge-Kutta method, meaning that it approximately preserves the geometric structure of Hamiltonian systems. We show this new rule obtains bounded orbits and time-average convergence to the set of Nash equilibria with arbitrary learning rates.

Our proposed rule, IMGD, can be viewed as an averaging of GDA (explicit Euler) and IGD (implicit Euler). We express this update rule in three ways: (2 ) is written consistently with the notation used to approximate dynamics as given in Section 2.4. (3 ) is expressed in stacked notation, which will significantly improve the readability of our analysis. Finally, the expression in Proposition 2, shows that (2 ) can be expressed as two half-steps followed by a normalization factor to maintain bounded orbits. This formulation is useful in practical settings as it expresses (2 ) as a sequence of separable updates based only on local information available to each agent, i.e., agents do not have to share information when updating.

\[x_i^{t+1}= x_i^t + \eta \frac{\nabla_{x_i}u_i(x^{t+1}) + \nabla_{x_i}u_i(x^{t})}{2}, \qquad i=1,2 \label{eq:IMGD}\tag{2}\]

To simplify our analysis, we also define our rule using stacked-variable notation and gradient vector field: \[x= \begin{pmatrix} x_1\\ x_2 \end{pmatrix}, \qquad \begin{pmatrix} \nabla_{x_1}u_1(x)\\ \nabla_{x_2}u_2(x) \end{pmatrix}=Gx+c,\qquad G= \begin{pmatrix} 0 & A\\ -A^\top & 0 \end{pmatrix}, \qquad c= \begin{pmatrix} b_1\\ b_2 \end{pmatrix}.\] yielding \[\label{eqn:StackedIMGD} \begin{align} x^{t+1} &= x^t + \frac{\eta}{2} \Bigl( (Gx^{t+1}+c) + (Gx^t+c) \Bigr)\\ &= \left(I-\frac{\eta}{2}G\right)^{-1} \left( \left(I+\frac{\eta}{2}G\right)x^t +\eta c \right)\\ &=\Phi_\eta x^t + D \end{align}\tag{3}\] where \(\Phi_\eta = \big(I - \tfrac{\eta}{2} G\big)^{-1} \big(I + \tfrac{\eta}{2} G\big), \quad D =\big(I - \tfrac{\eta}{2} G\big)^{-1} \eta c.\) In the proof of Proposition 2, we establish that \(\big(I - \tfrac{\eta}{2} G\big)^{-1}\) is well-defined.

Finally, we present (2 ) as two half steps followed by a normalization factor. This presentation shows that (2 ) can be expressed as a separable update rule.

3.1 Fixed Points for 2↩︎

First, we define the fixed points for (2 ), and we show they are equivalent to the set of Nash equilibria, which is an essential property for any learning algorithm.

3.2 Obtaining 2 via Symplectic Integration↩︎

Next, we show this update rule is naturally derived by applying the implicit midpoint integration technique, a known symplectic integrator, to continuous-time gradient descent. With symplectic integrators approximately preserving geometric structure and energy in Hamiltonian mechanics over long time horizons when learning rates are sufficiently small [26], we expect 2 to behave similarly to continuous-GD. Remarkably, in Section 3.3, we will uncover that 2 perfectly preserves energy with arbitrary learning rates.

We remark that the resulting discretization is algebraically identical to the trapezoidal rule for approximating dynamical systems in this (bilinear) setting. We prefer the interpretation as the implicit midpoint rule as it is a symplectic integrator and therefore preserves key geometric features of the continuous-time learning dynamics, allowing the discrete-time method to inherit important qualitative properties of continuous-time systems. The trapezoidal rule is not always a symplectic integrator and therefore is less likely to preserve properties in nonlinear systems.

3.3 Preserving the Orbits from Continuous-GD↩︎

In bilinear zero-sum games, the underlying continuous-time dynamics are governed by a skew-symmetric operator, leading to energy-preserving rotational behavior. Consequently, trajectories remain on level sets of a conserved quantity and are inherently bounded. This suggests that a well-designed discretization should inherit this qualitative property.

2 achieves this exactly: it preserves a discrete energy function, ensuring that the iterates evolve on invariant level sets and remain bounded for all step sizes. In contrast, alternating gradient descent, which is obtained from the symplectic Verlet-Störmer integration technique, only approximately preserves this energy [20], [22]; the AGD update introduces perturbations, which can lead to energy growth and divergence of the iterates when learning rates are too large.

3.4 Convergence↩︎

We now establish the ergodic convergence of the proposed method.

In particular, the learning rate \(\eta\) improves the convergence rate in the sense that the leading term scales as \(O(1/(\eta T))\) and there is no restriction on the learning rate \(\eta > 0\) for stability or convergence of the method. However, as \(\eta \to \infty\), the rate saturates at \(O\!\left(\frac{1}{T}\right).\) Thus, increasing \(\eta\) beyond a certain scale does not further improve the convergence rate in this analysis.

Moreover, when there exist multiple Nash equilibria, the solution set forms an affine manifold rather than an isolated point. In this setting, convergence can no longer be expected to a unique point. Instead, a natural notion of convergence is to the closest equilibrium within the solution set.

The proof is identical to the same result for Continuous-time GD in [27]. The key to the proof is that, since the strategies remain equidistant to every Nash equilibrium (Theorem 1), both the strategies and the time-average of the strategies are contained in an affine space containing only the equilibrium \(P(x^0)\), as depicted in Figure 1. Since the time-average of the strategies converges to the set of Nash (Theorem 2), and since \(P(x^0)\) is the equilibrium contained in an affine space containing the time-average of the strategies, the time-average must converge to \(P(x^0)\). For completeness, we include a near-verbatim of the proof from [27] in the appendix.

4 Two-Step Averaging and Decomposable Geometries↩︎

With the powerful properties of 2 established, we now turn our attention to its learning-rate selection and geometric behavior. Surprisingly, we find that with \(\eta\to\infty\), 2 can achieve time-average convergence after one iteration.

We prove this by decomposing the update rule into separable 2D dynamics, and base our analysis on these 2D projections. In Section 4.1, we introduce the decomposition of 2 into separable 2D dynamics; in Section 4.2, we find that 2 reduces to a planar rotation in 2D; in Section 4.3 we prove the fast convergence behavior in 2D; and in Section 4.4, we extend the result to higher dimensions and give the proof of Theorem 4.

4.1 Decomposition of 2 into Separable 2D Dynamics↩︎

First, to decompose a system of higher dimensions into separable 2-dimensional systems, we introduce the canonical form of a real skew-symmetric matrix:

4.1.0.1 Real Schur decomposition [28]

Let \(G \in \mathbb{R}^{k\times k}\), \(G^\top = -G\), \(k=k_1+k_2\). There exists an orthogonal matrix \(Q \in \mathbb{R}^{k\times k}\), i.e., \(Q^\top Q = 1\), such that \[G = Q J Q^\top ,\] where \(J\) is block diagonal of the form \[J = \begin{pmatrix} J(\omega_1) & & & \\ & \ddots & & \\ & & J(\omega_l) & \\ & & & 0 \end{pmatrix}, \quad J(\omega_j)= \begin{pmatrix} 0 & \omega_j \\ -\omega_j & 0 \end{pmatrix}, \qquad \omega_j \in \mathbb{R}.\]

and \(\lambda = \pm i\omega_j (j=1,\dots,l)\) are eigenvalues, zero eigenvalues corresponding to the zero block in \(J\).

The structure of J provides insight into analyzing the dynamics on separable 2 × 2 subspaces. And we find that 2 dynamics can decompose into separable 2D dynamics:

4.2 Planar Rotations in 2D↩︎

Now that we have established the system can be decomposed into separable 2D systems, we study the dynamics in 2D. We begin by establishing that dynamics rotate around the Nash equilibrium with constant velocity that is determined as a function of the learning rate.

Consider the update rule in Theorem 5, in every 2D systems, with \(G= \begin{pmatrix} 0 & \omega\\ -\omega & 0 \end{pmatrix}\), \(c= \begin{pmatrix} b_1\\ b_2 \end{pmatrix}\) and \(D =\big(I - \tfrac{\eta}{2} G\big)^{-1} \eta c\), the 2 update is

\[\label{eq:updatein2}x^{t+1} = \Phi_\eta x^t+D, \quad \Phi_{\eta} = \begin{pmatrix} 1 & -\frac{\eta \omega}{2}\\ \frac{\eta \omega}{2} & 1 \end{pmatrix}^{-1} \begin{pmatrix} 1 & \frac{\eta \omega}{2}\\ -\frac{\eta \omega}{2} & 1 \end{pmatrix} =\frac{1}{1+\frac{\eta^2 \omega^2}{4}} \begin{pmatrix} 1-\frac{\eta^2 \omega^2}{4} & \eta \omega\\ -\eta \omega & 1-\frac{\eta^2 \omega^2}{4} \end{pmatrix}.\tag{4}\]

Figure [fig:rotation] illustrates this geometric interpretation. Further, as \(\eta \to \infty\), the rotation approaches \(\pi\) as depicted in Figure [fig:largreta]. Formally, \(|\theta| = |2\arctan\!\left(\frac{\eta \omega}{2}\right)| \to \pi\).

4.3 Convergence After One Iteration in 2D↩︎

We can directly connect the angular movement to ergodic convergence rates after a single iteration. Specifically, since the rotation approaches \(\pi\) as \(\eta\to \infty\), the dynamics approximately start from the starting location in the second iteration. We show that larger learning rates drive the two-step average iterate closer to the Nash equilibrium, with approximation error proportional to \(\frac{1}{\eta}\):

4.4 Convergence After One Iteration in Arbitrary Dimension↩︎

Combining Theorems 4.1 and 4.3, we show 2 obtains time-average convergence after 2 iterations in higher dimensions.

5 Numerical Experiments↩︎

In this section, we evaluate the performance of 2 by comparing it with two other standard rules, Alternating Gradient Descent (AGD) and Optimistic Gradient Descent (OGD). In Section 5.1, we first discuss the choice of learning rate and run several experiments for 2 as it has no limits for \(\eta\). Then in Section 5.2, we run these three methods together for both fixed iterations and fixed time, and compare the outcomes. We don’t run dedicated experiments for analyzing \(\bar x^2\) since there’s no baseline, as 2 is the only known method to yield arbitrarily good approximations in a single iteration; however, the two-step time-average convergence still naturally appears in our analysis of IMGD for large learning rates.

For our benchmark methods, AGD is a classical gradient-based method as introduced in Section 2.4, and OGD is also a widely studied first-order method for zero-sum games, which incorporates information from the previous gradient to anticipate the rotational behavior of game dynamics: \[\begin{align} x_1^{t+1} &= x_1^t + 2\eta \nabla_{x_1}u_1(x_1^t,x_2^t) - \eta \nabla_{x_1}u_1(x_1^{t-1},x_2^{t-1}), \\ x_2^{t+1} &= x_2^t + 2\eta \nabla_{x_2}u_2(x_1^t,x_2^t) - \eta \nabla_{x_2}u_2(x_1^{t-1},x_2^{t-1}). \end{align}\]

For settings in the experiments, we solve a two-player zero-sum game: \[\begin{align} \max_{x_1 \in X_1} \; \min_{x_2 \in X_2} \; x_1^\top A x_2 + b_1 \cdot x_1 - b_2 \cdot x_2, \end{align}\] where \(x_1 \in X_1 \subset \mathbb{R}^{k_1}\) , \(x_2 \in X_2 \subset \mathbb{R}^{k_2}\), \(A \in \mathbb{R}^{k_1 \times k_2}\) and \(b_1 \in \mathbb{R}^{k_1}\), \(b_2 \in \mathbb{R}^{k_2}\) are fixed vectors. The payoff matrix \(A \in \mathbb{R}^{k_1 \times k_2}\) has entries independently sampled from the uniform distribution \([-1,1]\), while the vectors \(b_1 \in \mathbb{R}^{k_1}\) and \(b_2 \in \mathbb{R}^{k_2}\) are independently drawn from the standard normal distribution. The problem dimensions are set as \(k_1 = k_2 = 20\), giving a total dimension \(k = 40\).

All numerical experiments were conducted on a Lenovo ThinkPad L14 Gen 5 laptop equipped with an AMD Ryzen 5 PRO 7535U processor (2.90 GHz) and 32 GB of RAM, running a 64-bit Windows operating system.

5.1 Choice of Learning Rate↩︎

The theoretical learning rate bounds for OGD and AGD respectively are \(\eta_{\text{OGD}} \le 1/2 \|A\|\) [10] and \(\eta_{\text{AGD}} < 2/\|A\|\) [20], where \(\|A\|\) denotes the spectral norm of \(A\). We have matrix \(A \in \mathbb{R}^{k_1 \times k_2}\) sampled from the uniform distribution \([-1,1]\), so \(||A||\leq k_1\) (\(k_1 = k_2\)) [20]. Following the theoretical step-size scaling, we parameterize the learning rates as \(\eta_{\mathrm{AGD}} = c\,\frac{2}{k_1}, \quad \eta_{\mathrm{OGD}} = c\,\frac{1}{2k_1},\) where \(c \in \{2^{0},2^{-1},2^{-2},\ldots\}\). We perform a parameter sweep over these values, with the complete tuning results reported in Appendix 9. Among the tested values, we find that \(c=2^0\) consistently yields the best empirical performance. Therefore, we set the learning rate to be \[\eta_{\mathrm{AGD}}=\frac{2}{k_1}, \qquad \eta_{\mathrm{OGD}}=\frac{1}{2k_1}.\]

For 2 , there is no explicit upper bound on the learning rate. To select a suitable step size, we similarly parameterize it as \(\eta_\mathrm{IMGD} = \frac{c}{k_1}\) where \(c > 0\) is a tunable scalar. We consider a wide range of constants \(c \in \{1, 3, 5, 10, 20, 50, 80, 10^5, 10^9, 10^{16}\}\), and repeat each experiment 30 times. At each iteration, we compute the time-averaged iterate and record the distance to the set of Nash equilibria, \(\|G\bar{x}^T +c\|\), (in log scale).

The following figures present the convergence curves for 1000 and 10 iterations. Figure 3 (a) is highly oscillatory over the \(1000\)-iteration horizon, and Figure 3 (b) clearly reveals the underlying convergence behavior in \(10\) iterations.

From Figure 3, it is noticeable that 2 converges after 2 iterations for an extremely large learning rate, as indicated by Theorem 4. However, these good estimates also make it difficult to inspect the worst-case convergence. Thus, we examine long-term behavior by replacing each point with the maximum over future iterations to suppress oscillations and highlight convergence trends. Running the algorithm for 1000 iterations each time, the result is shown in Figure 4.

As \(c\) increases, we observe that the convergence becomes faster, and very large step sizes (e.g., \(c \ge 10^5\)) still exhibit stable behavior and achieve substantially lower error. This highlights a key advantage of 2 : it can tolerate, and even benefit from, very large step sizes, unlike existing methods.

As indicated by Theorem 2, the improvement saturates beyond a certain point, and extremely large values of \(c\) provide diminishing returns while potentially introducing numerical sensitivity. Therefore, for a fair and representative comparison across methods, we select \(c = 10\), as shown in Figure 4, which achieves a good balance between convergence speed, stability, and robustness. So we set the learning rate of 2 to be \[\eta_{\mathrm{IMGD}} = \frac{10}{k_1}.\]

5.2 Experiments↩︎

In the experiments, we evaluate the convergence behavior when all algorithms are run for a fixed number of iterations or time.

5.2.1 Fixed Iteration↩︎

We set the initial point \(x^0 \in \mathbb{R}^n\) to be sampled uniformly from \([-1,1]^n\). All algorithms are initialized from the same point within each trial to ensure a fair comparison. Each algorithm is run for a fixed number of iterations \(T = 1000,2000,3000\). We evaluate performance using the time-averaged iterate and measure its distance to a Nash equilibrium \(x^\star\). We repeat the experiment over 30 independent runs with different random initializations. Consistent with our theoretical analysis, we measure the distance to set of Nash equilibria with \(||G\bar{x}^T+c||\) and report a 95% confidence interval for this measured deviation from Nash.

From Figure 5 and Table 2, across all iteration budgets, 2 consistently achieves the smallest distance to equilibrium and exhibits steady improvement as the number of iterations increases. In contrast, both OGD and AGD converge more slowly and show significantly larger residual errors.

As the number of iterations increases from 1000 to 3000, the performance gap becomes more pronounced according to Table 3. The relative distance ratios indicate that AGD is approximately \(3\sim 5\times\) worse than 2 , while OGD can be more than an order of magnitude worse, reaching over \(10\times\) at 3000 iterations. These results highlight the superior stability and convergence efficiency of 2 in this setting.

5.2.2 Fixed Time↩︎

We also evaluate all algorithms under a fixed wall-clock time budget. Specifically, each method is allowed to run for \(t=0.05s , 0.2s, 0.5s\), and the number of iterations is determined dynamically within the time constraint. For 2 , we include the preprocessing time required to construct the matrices used in each iteration within the total time budget. This ensures a fair comparison by accounting for both setup cost and iterative computation.

Other settings are the same as fixed iteration. Performance is reported as absolute value of the time-averaged iterate gradient, together with a 95% confidence interval. Results are shown in Figure 6, Table  4 and Table  5.

Across all time budgets, 2 consistently achieves the smallest distance to equilibrium and shows steady improvement as the time budget increases. AGD and OGD converge more slowly with significantly larger residual errors. Relative performance shows that AGD is approximately \(4\sim5\times\) worse than 2 , while OGD can be \(20\times\) worse, highlighting the superior stability and efficiency of 2 .

These experiments indicate that, on the tested bilinear instances, IMGD achieves substantially smaller residuals than AGD and OGD under both fixed-iteration and fixed-time budgets.

6 Discussion↩︎

6.1 Computational Costs↩︎

The main additional computational cost of IMGD, relative to explicit methods such as AGD and OGD, is a one-time inverse computation. Specifically, to effectively implement IMGD, we precompute \(\phi_\eta\) and \(D\) in the update \(x^{t+1}=\Phi_\eta x^t+D\).

For dense games with total dimension \(k=k_1+k_2\), this initialization has cost \(O(k^3)\) using standard dense linear algebra. After initialization, each IMGD iteration requires a dense matrix-vector multiplication and vector addition, giving per-iteration cost \(\Theta(k^2)\). This is the same per-iteration order as standard dense implementations of AGD and OGD, whose updates also require matrix-vector products involving \(A\), \(A^\top\), or \(G\).

Consequently, the inverse cost is naturally amortized over the run. If an algorithm is run for \(T\) iterations, then the iterative cost is \(\Theta(Tk^2)\), while the one-time inverse cost is \(O(k^3)\). Therefore, in the common regime where \(T\) is significantly larger than \(k\), the initialization cost is asymptotically dominated by the cumulative iteration cost. This is the relevant regime for many online optimization experiments, where the number of gradient updates is typically much larger than the problem dimension.

This comparison is also consistent with our numerical experiments. In the fixed-time experiments, the preprocessing time required by IMGD is included in the total wall-clock budget, so the reported performance already accounts for the inverse computation rather than only the cheaper post-initialization iterations. Despite this accounting, IMGD remains substantially more accurate than AGD and OGD under the tested time budgets.

Finally, the absence of an upper stability restriction on \(\eta\) reduces another source of computational overhead. Methods such as AGD and OGD require learning rates satisfying bounds depending on \(\|G\|\), and practical implementations often require additional tuning or norm estimation. IMGD is well-defined for every \(\eta>0\), so its computational cost should be viewed as a tradeoff: a one-time inverse computation in exchange for learning-rate robustness, exact orbit preservation, and strong empirical performance.

6.2 Extensions to Normal-Form Games↩︎

As discussed in Section 2.2, the unconstrained bilinear setting serves as a canonical model for the local behavior of learning dynamics in normal-form zero-sum games. Once the equilibrium support is identified, the dynamics are confined to the affine hull of the corresponding face of the simplex. After eliminating the affine equality constraints, the resulting dynamics are governed by an unconstrained bilinear game, so the analysis developed in this paper applies directly on the reduced subspace.

For explicit methods such as AGD and OGD, this reduction is largely conceptual: the update rule itself remains unchanged, while the analysis is performed on the active face. IMGD introduces an additional challenge. The implicit midpoint update requires computing an inverse involving the game matrix. In the constrained setting, this inverse is formed only over the coordinates corresponding to the current support of each player’s strategy. Consequently, whenever the support changes, the reduced system changes as well, requiring the inverse to be recomputed.

Recomputing this inverse after every support change would limit the computational advantage obtained by preprocessing the inverse once. This suggests that IMGD is likely to be less effective during the initial exploratory phase of learning, when supports change frequently.

Instead, the structure of IMGD naturally suggests a hybrid approach. A first-stage learning algorithm, such as AGD, OGD, or another support-identification procedure, could be used to rapidly identify the equilibrium support. Once a candidate support has stabilized, IMGD could be initialized on the induced subgame to exploit its favorable geometric properties and rapidly converge while preserving bounded orbits and allowing arbitrarily large learning rates.

This perspective raises several interesting research questions. Can equilibrium support be identified reliably during learning? When is it safe to transition from an exploratory learning dynamic to IMGD? Can support changes be detected cheaply enough that inverse updates remain computationally practical? Addressing these questions may lead to learning algorithms that combine the fast exploratory behavior of explicit methods with the strong stability guarantees of structure-preserving implicit dynamics.

7 Conclusion↩︎

We introduced implicit midpoint gradient descent (IMGD), an implicit update rule for unconstrained bilinear two-player zero-sum games. The method is obtained by applying the implicit midpoint integrator to the continuous-time gradient dynamics, yielding a discrete-time learning rule that preserves the underlying energy structure exactly. We proved that IMGD preserves the Euclidean distance to every Nash equilibrium, has bounded orbits for every learning rate \(\eta>0\), and achieves \(O(1/T)\) ergodic convergence without an upper stability restriction on the step size.

We also showed that the dynamics decompose into independent two-dimensional rotations under the real Schur decomposition of the skew-symmetric game matrix. This geometric viewpoint explains the behavior of the method in the large-learning-rate regime and shows that a two-step average approaches the Nash residual at rate \(O(1/\eta)\). Numerical experiments further demonstrate that IMGD can substantially outperform AGD and OGD in both fixed-iteration and fixed-time comparisons on the tested bilinear instances.

Overall, the results suggest that structure-preserving discretizations provide a useful algorithmic principle for zero-sum game dynamics. The method trades a one-time inverse computation for learning-rate robustness and exact orbit preservation, motivating future work on sparse implementations, adaptive learning-rate selection, and support-identification methods for constrained normal-form games.

8 Proof of Theorem 3: Time-Average Convergence↩︎

9 Learning Rate Selection for AGD and OGD↩︎

Consider \(\eta_{\mathrm{AGD}} = c\,\frac{2}{k_1}, \quad \eta_{\mathrm{OGD}} = c\,\frac{1}{2k_1},\) where \(c \in \{2^{0},2^{-1},2^{-2},2^{-3},2^{-4}\}\). For each algorithm, we evaluate all candidate values of \(c\) on the experimental instances used in the paper. The corresponding convergence curves are shown in Figures 7 and 8. In both cases, \(c=1\) consistently provides the best overall empirical performance.

References↩︎