\nabla_{y} f) = (0,0)$) and assume furthermore that $\lambda_{\min} {\coloneqq}\min \left(\lambda_{\min} \left(\eta D_{xx}^2 f + \bar{D}\right), \lambda_{\min} \left(-\eta D_{yy}^2 f + \bar{D} \right) \right) > 0$ and $f \in C^2(\mathbb{R}^{m + n})$ with Lipschitz continuous mixed Hessian. Then there exists a neighbourhood $\mathcal{U}$ of $(x^*,y^*)$, such that CGD started in $(x_{1},y_{1}) \in \mathcal{U}$ converges to a point in $\mathcal{U}$ at an exponential rate that depends only on $\lambda_{\min}$.

The results on local stability for existing modifications of GDA, including those of [@mescheder2017numerics; @daskalakis2017training; @mertikopoulos2019optimistic] (see also [@liang2018interaction]) all require the stepsize to be chosen inversely proportional to an upper bound on $\sigma_{\max} (D_{xy}^2f)$ and indeed we will see in our experiments that the existing methods are prone to divergence under strong interactions between the two players (large $\sigma_{\max}(D_{xy}^2f)$). In contrast to these results, our convergence results *only improve* as the interaction between the players becomes stronger.

**Why not use $D_{xx}^2f$ and $D_{yy}^2g$?:** The use of a bilinear approximation that contains some, but not all second order terms is unusual and begs the question why we do not include the diagonal blocks of the Hessian in Equation  resulting in the damped and regularized Newton’s method $$\label{eqn:newton}
    \begin{pmatrix}
         \Delta x\\
         \Delta y
    