# Local convergence of min-max algorithms to differentiable equilibrium on
Riemannian manifold

Sixin Zhang
Université de Toulouse, INP, IRIT
2, rue Camichel, BP 7122, 31071 Toulouse Cedex 7, France
<EMAIL_ADDRESS>

###### Abstract

We study min-max algorithms to solve zero-sum differentiable games on
Riemannian manifold. The notions of differentiable Stackelberg equilibrium and
differentiable Nash equilibrium in Euclidean space are generalized to
Riemannian manifold, through an intrinsic definition which does not depend on
the choice of local coordinate chart of manifold. We then provide sufficient
conditions for the local convergence of the deterministic simultaneous
algorithms $\tau$-GDA and $\tau$-SGA near such equilibrium, using a general
methodology based on spectral analysis. These algorithms are extended with
stochastic gradients and applied to the training of Wasserstein GAN. The
discriminator of GAN is constructed from Lipschitz-continuous functions based
on Stiefel manifold. We show numerically how the insights obtained from the
local convergence analysis may lead to an improvement of GAN models.

##  1 Introduction

Riemannian min-max problem has attracted a lot of research attention in recent
years, with various machine learning applications including robust PCA [22],
robust neural network training [20] and generative adversarial networks (GANs)
[17]. This problem is formalized as a two-player zero-sum game, where the
variables of each player are constrained on the Riemannian manifold
$\mathcal{M}_{1}$ and $\mathcal{M}_{2}$,

$\min_{x\in\mathcal{M}_{1}}\max_{y\in\mathcal{M}_{2}}f(x,y).$

When $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$ are Euclidean, gradient-based
methods such as gradient-descent-ascent (GDA) [10, 25], extra-gradient [15,
27], optimistic mirror descent [29], Hamiltonian-gradient descent [26] and
symplectic gradient-adjustment (SGA) [4] are often considered to solve this
problem. When $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$ are Riemannian, how to
extend existing algorithms from Euclidean space to Riemannian manifold and how
to analyze their convergence become an interesting topic in recent years. In
Table 1, we summarize several related works on the Riemannian min-max problem.
In particular, we observe that under suitable assumptions of the problem, one
can obtain the global convergence of a suitable algorithm towards certain
limit set such as Nash equilibrium. One major issue is that these assumptions
typically do not hold in applications such as GANs [33]. However, when a min-
max problem is non-convex and non-concave, even in the Euclidean case, the
limit set of GDA algorithm can be very complicated [18]. As a consequence, we
aim to develop suitable concepts of limit set and algorithms with local
convergence guarantees for an interesting class of problems [21].

In this article, we study the convergence of simultaneous min-max algorithms
in order to find differentiable Stackelberg equilibrium (DSE) or
differentiable Nash equilibrium (DNE). When $\mathcal{M}_{1}$ and
$\mathcal{M}_{2}$ are Euclidean, these two limit sets are generic among local
Stackelberg or Nash equilibrium of smooth $f$ [13]. The main contribution of
this article is to generalize these concepts to the Riemannian min-max
problem, and to analyze the local convergence of two representative
deterministic min-max algorithms to such equilibrium. We further develop these
algorithms using stochastic-gradients and analyze numerically their behavior
in the training of Wasserstein GANs.

Section 2 provides the definition of DSE and DNE when $\mathcal{M}_{1}$ and
$\mathcal{M}_{2}$ are finite-dimensional Riemannian manifolds. By using a
local coordinate chart near such equilibrium, we generalize these concepts in
Euclidean space [13, 38], and show that their definitions only reply on
intrinsic quantities of the manifolds. Using the same local coordinate chart,
we then provide a methodology to analyze the local convergence of
deterministic simultaneous algorithms to a stationary point in Section 3.1. It
is based on the classical Ostrowski Theorem to analyze the stability of a
fixed point in one-step stationary methods [31]. This methodology is general
in the sense that the problem is reduced to the spectral analysis of a
Jacobian matrix in Euclidean space, which is at the heart of analyzing various
min-max algorithms [10, 3, 12, 40, 24, 11]. We focus on two algorithms
$\tau$-GDA and $\tau$-SGA, where the variables $(x,y)$ are updated
simultaneously at each iteration with a learning rate ratio $\tau$. In Section
3.2, sufficient conditions on $\tau$ are given to ensure the local convergence
of $\tau$-GDA to DSE and DNE. One issue of $\tau$-GDA is its slow convergence
rate when there is strong rotational dynamics [24]. This is a known phenomenon
near a Nash equilibrium in bilinear games due to the competition between the
two players [4]. To improve the convergence rate, first-order methods such as
extra-gradient and optimistic mirror descent are often used to correct the
gradient direction of $\tau$-GDA. Recently, these methods are extended to
Riemannian manifold [41, 19, 37]. However, their computations rely on the
exponential map and parallel transport which can be costly [1]. In Section
3.3, we consider the $\tau$-SGA algorithm which addresses the same rotational
problem based on auto-differentiation. It naturally extends the SGA algorithm
[14, 4, 23] to Riemannian manifold. We further provide an asymptotic
convergence analysis of $\tau$-SGA to DSE which is valid when $\tau$ is large.
In Section 4, we apply stochastic $\tau$-GDA and $\tau$-SGA to train
Wasserstein GAN. The underlying min-max problem is Riemannian since we impose
Stiefel manifold constraints on $y$ to construct Lipschitz-continuous
discriminators. We consider two examples, where one of the discriminators is
build upon the wavelet scattering transform to capture discriminative
information in MNIST dataset [8]. We illustrate a local convergence dilemma
faced by $\tau$-GDA, and show how $\tau$-SGA may lead to an improved local
convergence and generator quality of GAN compared to $\tau$-GDA. Section 5
concludes.

Table 1: Related works of Riemannian min-max problem Ref. | Class of problem | Limit set | Algorithm
---|---|---|---
[41] | Geodesically convex $\mathcal{M}_{1}$,$\mathcal{M}_{2}$ | Nash equilibrium | Extra-gradient
|  $f$ is geodesically convex/concave (quasi), semi-continuous (lower/upper)
[22] | +Bounded $\mathcal{M}_{1}$,$\mathcal{M}_{2}$ | Nash equilibrium | Extra-gradient
|  $f$ is geodesically convex/concave, smooth
[17] | Complete $\mathcal{M}_{1}$,$\mathcal{M}_{2}$ | Stationary point of $f$ | Hamiltonian
|  $(x,y)\mapsto\|\mbox{grad}_{x}f\|_{x}^{2}+\|\mbox{grad}_{y}f\|_{y}^{2}$ is Polyak-Łojasiewicz | (consensus)
[20] | Euclidean $\mathcal{M}_{2}$ | Stationary point of $\max_{y}$ f | GDA
|  $f$ is strongly concave in $y$ | (accelerated)
This work |  $f$ is twice continuously differentiable | Differentiable equilibrium | GDA, SGA

##  2 Differentiable equilibrium on Riemannian manifold

In this article, we assume that
$f:\mathcal{M}_{1}\times\mathcal{M}_{2}\to\mathbb{R}$ is twice continuously
differentiable on the product manifold $\mathcal{M}_{1}\times\mathcal{M}_{2}$.
The dimension of $\mathcal{M}_{1}$ is $d_{1}$, and the dimension of
$\mathcal{M}_{2}$ is $d_{2}$. This section introduces the notions of
differentiable Stackelberg differentiable and differentiable Nash equilibrium,
and gives some simple examples to illustrate their difference.

###  2.1 Differentiable Stackelberg equilibrium (DSE)

To define the notion of DSE on Riemannian manifold, we first introduce the
definition of DSE in a local coordinate chart of
$\mathcal{M}_{1}\times\mathcal{M}_{2}$. We then show that it is equivalent to
an intrinsic definition which does not depend on the choice of the coordinate
chart.

##### Local definition of DSE

Let $(O_{1}\times O_{2},\varphi_{1}\times\varphi_{2})$ be a local coordinate
chart around $(x^{\ast},y^{\ast})\in\mathcal{M}_{1}\times\mathcal{M}_{2}$. We
rewrite the function $f(x,y)$ on $(x,y)\in O_{1}\times O_{2}$ using this chart
by

$\bar{f}(u_{1},u_{2})=f(\varphi^{-1}_{1}(u_{1}),\varphi^{-1}_{2}(u_{2})).$

The next definition of DSE based on $\bar{f}$ is consistent with the standard
definition of DSE in Euclidean space [13, 38]. Let $\partial$ and
$\partial^{2}$ denote the first-order and second-order partial derivatives of
a function on Euclidean space.

######  Definition 2.1.

We say that $(x^{\ast},y^{\ast})$ is a DSE of $f$ if
$(u_{1}^{\ast},u_{2}^{\ast})=(\varphi_{1}(x^{\ast}),\varphi_{2}(y^{\ast}))$ is
a DSE of $\bar{f}$, i.e.

$\displaystyle\partial_{u_{1}}\bar{f}(u_{1}^{\ast},u_{2}^{\ast})=0,\quad\partial_{u_{2}}\bar{f}(u_{1}^{\ast},u_{2}^{\ast})=0,$
(1)
$\displaystyle-\partial_{u_{2}u_{2}}^{2}\bar{f}(u_{1}^{\ast},u_{2}^{\ast})\quad\mbox{p.d
(abbrev. of positive definite)},$ (2)
$\displaystyle\left[\partial_{u_{1}u_{1}}^{2}\bar{f}-\partial_{u_{2}u_{1}}^{2}\bar{f}\cdot(\partial_{u_{2}u_{2}}^{2}\bar{f})^{-1}\cdot\partial_{u_{1}u_{2}}^{2}\bar{f}\right](u_{1}^{\ast},u_{2}^{\ast})\quad\mbox{p.d}.$
(3)

The condition (1) alone means that $(u_{1}^{\ast},u_{2}^{\ast})$ is a
stationary point of $\bar{f}$. The conditions (1),(2),(3) imply that
$(u_{1}^{\ast},u_{2}^{\ast})$ is a local minimax point of $\bar{f}$ [21, 38].
See further discussions in Appendix A.

##### Intrinsic definition of DSE

We denote the tangent space of $\mathcal{M}_{1}$ at $x\in\mathcal{M}_{1}$ by
$T_{x}\mathcal{M}_{1}$ and the tangent space of $\mathcal{M}_{2}$ at
$y\in\mathcal{M}_{2}$ by $T_{y}\mathcal{M}_{2}$. The following concepts are
standard in manifold optimization, and they are used in Riemannian min-max
problem [17].

$\bullet$ Riemannian gradient: $\mbox{grad}_{x}f(x,y)\in
T_{x}\mathcal{M}_{1},\mbox{grad}_{x}f(x,y)\in T_{y}\mathcal{M}_{2}.$
$\bullet$ Riemannian Hessian: $\mbox{Hess}_{x}f(x,y):T_{x}\mathcal{M}_{1}\to
T_{x}\mathcal{M}_{1},\mbox{Hess}_{y}f(x,y):T_{y}\mathcal{M}_{2}\to
T_{y}\mathcal{M}_{2}.$
$\bullet$ Riemannian cross-gradient:
$\mbox{grad}^{2}_{xy}f(x,y):T_{x}\mathcal{M}_{1}\to
T_{y}\mathcal{M}_{2},\mbox{grad}^{2}_{yx}f(x,y):T_{y}\mathcal{M}_{2}\to
T_{x}\mathcal{M}_{1}.$

######  Definition 2.2.

We say that $(x^{\ast},y^{\ast})$ is a DSE of $f$ if

$\displaystyle\mbox{grad}_{x}f(x^{\ast},y^{\ast})=0,\quad\mbox{grad}_{y}f(x^{\ast},y^{\ast})=0,$
(4) $\displaystyle-\mbox{Hess}_{y}f(x^{\ast},y^{\ast})\quad\mbox{p.d},$ (5)
$\displaystyle[\mbox{Hess}_{x}f-\mbox{grad}^{2}_{yx}f\cdot(\mbox{Hess}_{y}f)^{-1}\cdot\mbox{grad}^{2}_{xy}f](x^{\ast},y^{\ast})\quad\mbox{p.d}.$
(6)

The following result shows that the above two definitions are equivalent, no
matter how one choose the local coordinate chart nearby $(x^{\ast},y^{\ast})$,
or the Riemannian metric on $\mathcal{M}_{1}\times\mathcal{M}_{2}$.

######  Theorem 2.1.

The local definition of DSE in Definition 2.1 is equivalent to the intrinsic
definition of DSE in Definition 2.2.

The proof is given in Appendix B. The main idea of the proof is to use the
local coordinate chart $(O_{1}\times O_{2},\varphi_{1}\times\varphi_{2})$ to
represent the smooth vector fields $(x,y)\mapsto\mbox{grad}_{x}(x,y)$ and
$(x,y)\mapsto\mbox{grad}_{y}(x,y)$ around the point $(x^{\ast},y^{\ast})$, so
as to show the equivalence between (1)-(3) and (4)-(6).

###  2.2 Differentiable Nash equilibrium (DNE) and examples

In game theory, one is often interested in finding a Nash equilibrium since it
maintains a symmetry between the role of the two players. Based on the
intrinsic definition of DSE, we immediately obtain an intrinsic definition of
DNE. We then provide some simple examples of DSE and DNE.

######  Definition 2.3.

We say that $(x^{\ast},y^{\ast})$ is a DNE of $f$ if

$\displaystyle\mbox{grad}_{x}f(x^{\ast},y^{\ast})=0,\quad\mbox{grad}_{y}f(x^{\ast},y^{\ast})=0,$
(7)
$\displaystyle-\mbox{Hess}_{y}f(x^{\ast},y^{\ast})\quad\mbox{p.d},\quad\mbox{Hess}_{x}f(x^{\ast},y^{\ast})\quad\mbox{p.d}.$
(8)

The definition of DNE generalizes an equivalent notion “strongly local min-max
point” in Euclidean space [10, Definition 1.6]. It is clear that if
$(x^{\ast},y^{\ast})$ is a DNE, then it is also a DSE.

####  2.2.1 Example 1: DSE

Consider $f(x,y)=\langle y,Ax-b\rangle$, with
$x\in\mathcal{M}_{1}=\mathbb{R}^{d_{1}}$ and $y\in\mathcal{M}_{2}=S^{d_{2}}$.
The manifold $S^{d_{2}}$ is the unit sphere embedded in
$\mathbb{R}^{d_{2}+1}$, endowed with the Euclidean metric
$\langle\cdot,\cdot\rangle$ on $\mathbb{R}^{d_{2}+1}$. Let $A^{+}$ denote the
pseudo-inverse of $A$. We next provide a sufficient condition on the existence
of DSE.

######  Proposition 2.1.

Assume $b\not\in\mbox{Range}(A)$, $\mbox{Ker}(A)=\\{0\\}$. Let
$x^{\ast}=A^{+}b$, $y^{\ast}=\frac{Ax^{\ast}-b}{\|Ax^{\ast}-b\|}$, then
$(x^{\ast},y^{\ast})$ is a DSE of the $f$ in Example 1.

The proof is given in Appendix C. Since $\mathcal{M}_{1}$ is Euclidean, we
note that $\partial^{2}_{xx}f(x^{\ast},y^{\ast})=0$. This implies that the
$(x^{\ast},y^{\ast})$ is not DNE. But the $k$-th eigenvalue
$\lambda_{k}(\mathbf{A})$ of $\mathbf{A}=-\mbox{Hess}_{y}f(x^{\ast},y^{\ast})$
equals to $\|Ax^{\ast}-b\|>0$, since from the proof
$\langle\mathbf{A}[\eta^{\ast}],\eta^{\ast}\rangle=\|Ax^{\ast}-b\|\|\eta^{\ast}\|^{2}$
for any $\eta^{\ast}\in T_{y^{\ast}}\mathcal{M}_{2}$.

####  2.2.2 Example 2: DSE

Consider $f(x,y)=\langle y,Ax-b\rangle-\frac{\kappa}{2}\|Ax\|^{2}$, with
$x\in\mathcal{M}_{1}=\mathbb{R}^{d_{1}}$ and $y\in\mathcal{M}_{2}=S^{d_{2}}$.
Compared to the $f$ in Example 1, we add a quadratic function of $Ax$ with a
curvature parameter $\kappa>0$. We next show that if $\kappa$ is close to
zero, we can still find a DSE near the DSE in Proposition 2.1.

######  Proposition 2.2.

Assume $b\not\in\mbox{Range}(A)$, $\mbox{Ker}(A)=\\{0\\}$. There exists
$\kappa_{0}>0$ such that for any $0<\kappa<\kappa_{0}$, there is a number $c$
close to $1$, and $x^{\ast}=cA^{+}b$,
$y^{\ast}=\frac{Ax^{\ast}-b}{\|Ax^{\ast}-b\|}$, so that $(x^{\ast},y^{\ast})$
is a DSE of the $f$ in Example 2.

The proof is given in Appendix D. From the proof, we have
$\mathbf{C}=\partial^{2}_{xx}f(x,y)=-\kappa A^{\intercal}A$. The spectral
radius of $\mathbf{C}$ equals to its operator norm
$\|\mathbf{C}\|=\kappa\|A^{\intercal}A\|>0$. As Example 1, all the eigenvalues
of $\mathbf{A}=-\mbox{Hess}_{y}f(x^{\ast},y^{\ast})$ equal to
$\|Ax^{\ast}-b\|=\|cAA^{+}b-b\|$. These quantities will be used later to
analyze the local convergence of $\tau$-GDA and $\tau$-SGA.

####  2.2.3 Example 3: DNE

Consider $f(x,y)=\frac{1}{2}\|Ax+y-b\|^{2}$ with
$x\in\mathcal{M}_{1}=\mathbb{R}^{d_{1}}$ and $y\in\mathcal{M}_{2}=S^{d_{2}}$.
$\mathcal{M}_{2}$ is the same embedded sub-manifold of $\mathbb{R}^{d_{2}+1}$
as above. We next provide a sufficient condition on the existence of DNE. The
proof is given in Appendix E.

######  Proposition 2.3.

Assume $b\not\in\mbox{Range}(A)$, $\mbox{Ker}(A)=\\{0\\}$. Let
$x^{\ast}=A^{+}b$, $y^{\ast}=\frac{Ax^{\ast}-b}{\|Ax^{\ast}-b\|}$, then
$(x^{\ast},y^{\ast})$ is a DNE of the $f$ in Example 3.

##  3 Min-max algorithms for differentiable equilibrium

Gradient-based min-max algorithms such as GDA and SGA are often used to find
local Nash equilibrium when $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$ are
Euclidean [10, 23]. We study how to extend these algorithms to Riemannian
manifold and then analyze their local convergence to DSE and DNE.

###  3.1 Local convergence of deterministic simultaneous algorithms

We analyze the local convergence of one-step deterministic algorithms to a
stationary point. As the update rule of each algorithm in this family does not
change over iteration, we represent it in a local coordinate chart and then
apply the proof technique of the classical Ostrowski Theorem.

A deterministic simultaneous algorithm is defined by two vector fields
$(x,y)\mapsto\xi_{1}(x,y)$, $(x,y)\mapsto\xi_{2}(x,y)$ on $\mathcal{M}_{1}$
and $\mathcal{M}_{2}$, and a choice of manifold retractions. Initialized at a
point on $\mathcal{M}_{1}\times\mathcal{M}_{2}$, the algorithm modifies the
variables $(x,y)$ of the two players, through an update rule
$\mathbf{T}:\mathcal{M}_{1}\times\mathcal{M}_{2}\to\mathcal{M}_{1}\times\mathcal{M}_{2}$
of the following form,

$\mathbf{T}(x,y)=(\mathcal{R}_{1,x}(\xi_{1}(x,y)),\mathcal{R}_{2,y}(\xi_{2}(x,y))),$

where $\mathcal{R}_{1,x}:T_{x}\mathcal{M}_{1}\to\mathcal{M}_{1}$ (resp.
$\mathcal{R}_{2,y}:T_{y}\mathcal{M}_{2}\to\mathcal{M}_{2}$) denotes the
restriction of a retraction $\mathcal{R}_{1}$ at $x\in\mathcal{M}_{1}$ (resp.
retraction of $\mathcal{R}_{2}$ at $y\in\mathcal{M}_{2}$). For example, on the
$\mathcal{M}_{1}=\mathbb{R}^{d_{1}}$ and $\mathcal{M}_{2}=S^{d_{2}}$ in
Section 2.2, we will take $\mathcal{R}_{1,x}(\delta)=x+\delta$ for
$\delta\in\mathbb{R}^{d_{1}}$, and $\mathcal{R}_{2,y}(\eta)$ to be the
projection of the vector $y+\eta$ in $\mathbb{R}^{d_{2}+1}$ to the sphere
$\mathcal{M}_{2}$.

We next define what it means for an update rule $\mathbf{T}$ to be locally
convergent to a stationary point $(x^{\ast},y^{\ast})$, which implies that
$(x^{\ast},y^{\ast})$ is a point of attraction [31, Definition 10.1.1].

######  Definition 3.1 (locally convergent).

Let $(x^{\ast},y^{\ast})$ be a stationary point satisfying (4). There exists a
local stable region $S_{\delta}\subset\mathcal{M}_{1}\times\mathcal{M}_{2}$,
which is an open set containing $(x^{\ast},y^{\ast})$ such that started from
$(x(0),y(0))\in S_{\delta}$, the update rule generates a sequence of points
$(x(t+1),y(t+1))=\mathbf{T}(x(t),y(t))$ belonging to $S_{\delta},\forall t\geq
0$, and moreover $\lim_{t\to\infty}(x(t),y(t))=(x^{\ast},y^{\ast})$.

We next give a sufficient condition for $\mathbf{T}$ to be locally convergent
to $(x^{\ast},y^{\ast})$. It is based on the spectral radius of the following
linear transformation on $T_{x^{\ast}}\mathcal{M}_{1}\times
T_{y^{\ast}}\mathcal{M}_{2}$,

$\mathbf{T}^{\prime}(x^{\ast},y^{\ast})=\left(\begin{array}[]{ll}I+\nabla_{x}\xi_{1}(x^{\ast},y^{\ast})&D_{y}\xi_{1}(x^{\ast},y^{\ast})\\\
D_{x}\xi_{2}(x^{\ast},y^{\ast})&I+\nabla_{y}\xi_{2}(x^{\ast},y^{\ast})\end{array}\right),$
(9)

where $D_{x}$ and $\nabla_{x}$ (resp. $D_{y}$ and $\nabla_{y}$) denote the
differential operator (in a tangent map) and the Riemannian connection on
$\mathcal{M}_{1}$ (resp. on $\mathcal{M}_{2}$).

######  Theorem 3.1.

Let $(x^{\ast},y^{\ast})\in\mathcal{M}_{1}\times\mathcal{M}_{2}$ be a
stationary point of $f$. If $\xi_{1}$ and $\xi_{2}$ are differentiable on
$\mathcal{M}_{1}\times\mathcal{M}_{2}$, and the spectral radius of
$\mathbf{T}^{\prime}(x^{\ast},y^{\ast})$ is strictly smaller than one, then
$\mathbf{T}$ is locally convergent to $(x^{\ast},y^{\ast})$.

###### Proof.

The idea of the proof is illustrated in Figure 1. Inspired by the proof
technique in [7, Theorem 4.19], we first identify a set $\bar{S}$ so that the
induced dynamics $\bar{\mathbf{T}}$ is well-defined on $\bar{S}$. Let us
consider the set $S=\mathbf{T}^{-1}(O_{1}\times O_{2})\cap(O_{1}\times O_{2})$
and its local coordinate domain
$\bar{S}=(\varphi_{1}\times\varphi_{2})(S)\subset\bar{O}_{1}\times\bar{O}_{2}$.
We can therefore write both $(x,y)\in S\subset O_{1}\times O_{2}$ and
$\mathbf{T}(x,y)\in O_{1}\times O_{2}$ in local coordinate, i.e.
$\forall(u_{1},u_{2})\in\bar{S}$,

$\bar{\mathbf{T}}(u_{1},u_{2})=\varphi_{1}\times\varphi_{2}\circ(\mathbf{T}(\varphi^{-1}_{1}(u_{1}),\varphi^{-1}_{2}(u_{2}))).$

Note that $S$ and $\bar{S}$ are non-empty open sets since $\mathbf{T}$ is
continuous (by the continuity of the vector fields $\xi_{1}$, $\xi_{2}$ and
the retractions $\mathcal{R}_{1}$ and $\mathcal{R}_{2}$) and
$(x^{\ast},y^{\ast})$ is a fixed point of $\mathbf{T}$. We next identify a
local stable region $\bar{S}_{\delta}\subset\bar{S}$ around
$(u_{1}^{\ast},u_{2}^{\ast})$ (the local coordinate of $(x^{\ast},y^{\ast})$),
so that if $(u_{1},u_{2})\in\bar{S}_{\delta}$, then
$\bar{\mathbf{T}}(u_{1},u_{2})\in\bar{S}_{\delta}\subset\bar{O}_{1}\times\bar{O}_{2}$.
Let $S_{\delta}=(\varphi_{1}\times\varphi_{2})^{-1}\circ(\bar{S}_{\delta})$,
and assume $(x(0),y(0))\in S_{\delta}$, then the sequence $(x(t),y(t))\in
S_{\delta}\subset S,\forall t\geq 1$ by recursion. To construct such a region,
the key is to verify that the spectral radius $\rho$ of the Jacobian matrix
$\bar{\mathbf{T}}^{\prime}(u_{1}^{\ast},u_{2}^{\ast})$ is strictly smaller
than one (in Appendix F, we show that the eigenvalues of
$\mathbf{T}^{\prime}(x^{\ast},y^{\ast})$ are the same as
$\bar{\mathbf{T}}^{\prime}(u_{1}^{\ast},u_{2}^{\ast})$), so that according to
the proof of Ostrowski Theorem in [31, Section 10.1.3], for an arbitrary
$\epsilon>0$, there exists a norm $\|\cdot\|_{\epsilon}$ on
$\mathbb{R}^{d_{1}}\times\mathbb{R}^{d_{2}}$ and $\delta>0$, satisfying

$\|\bar{\mathbf{T}}(u_{1},u_{2})-(u_{1}^{\ast},u_{2}^{\ast})\|_{\epsilon}\leq(\rho+2\epsilon)\|(u_{1},u_{2})-(u_{1}^{\ast},u_{2}^{\ast})\|_{\epsilon},\quad\forall(u_{1},u_{2})\in\bar{S}_{\delta},$
(10)

where
$\bar{S}_{\delta}=\\{(u_{1},u_{2})|\|(u_{1},u_{2})-(u_{1}^{\ast},u_{2}^{\ast})\|_{\epsilon}<\delta\\}\subset\bar{S}$.
Let’s choose $\epsilon$ so that $\rho+2\epsilon<1$, then from (10), if
$(u_{1}(0),u_{2}(0))\in\bar{S}_{\delta}$, the sequence $(u_{1}(t),u_{2}(t))$
stays in $\bar{S}_{\delta}$ and converges to $(u_{1}^{\ast},u_{2}^{\ast})$ as
$t\to\infty$. As $\varphi_{1}\times\varphi_{2}$ is a continuous bijection
(homeomorphism) from $O_{1}\times O_{2}$ to $\bar{O}_{1}\times\bar{O}_{2}$, we
have equivalently that if $(x(0),y(0))\in
S_{\delta}=(\varphi_{1}\times\varphi_{2})^{-1}\circ(\bar{S}_{\delta})$, the
sequence $(x(t),y(t))$ will stay in $S_{\delta}$ and converges to
$(x^{\ast},y^{\ast})$. ∎

(a)

$\bar{S}$$\bar{O}_{1}\times\bar{O}_{2}$$S$$O_{1}\times
O_{2}$$\bar{\mathbf{T}}$$\varphi_{1}\times\varphi_{2}$$\mathbf{T}$$\varphi_{1}\times\varphi_{2}$

(b)

.$(u_{1}^{\ast},u_{2}^{\ast})$$\bar{S}_{\delta}$$\bar{S}$$\bar{O}_{1}\times\bar{O}_{2}$

Figure 1: Local convergence of a deterministic simultaneous algorithm. (a):
The update rule $\mathbf{T}$ on the manifold
$\mathcal{M}_{1}\times\mathcal{M}_{2}$ induces a local update rule
$\bar{\mathbf{T}}$ in the local coordinate system, defined by a chart
$(O_{1}\times O_{2},\varphi_{1}\times\varphi_{2})$ around
$(x^{\ast},y^{\ast})$. (b): The induced update rule $\bar{\mathbf{T}}$ is
defined in the local coordinate on the set $\bar{S}$ around
$(u_{1}^{\ast},u_{2}^{\ast})=\varphi_{1}\times\varphi_{2}\circ(x^{\ast},y^{\ast})$.
It contains a stable region $\bar{S}_{\delta}$.

###  3.2 Simultaneous gradient-descent-ascent algorithm ($\tau$-GDA)

The $\tau$-GDA algorithm uses the Riemannian gradients $\mbox{grad}_{x}f(x,y)$
and $\mbox{grad}_{y}f(x,y)$ to update $x$ and $y$ simultaneously. The local
convergence of deterministic $\tau$-GDA to DSE and DNE is known in Euclidean
space [10, 12, 21, 24]. This section extends these results to Riemannian
manifold. We obtain a sharp lower-bound of $\tau$ for $\tau$-GDA to be locally
convergent to DSE, based on a refinement of the spectral analysis in Euclidean
space [24].

#####  $\tau$-GDA algorithm

In deterministic setting, the update rule $\mathbf{T}$ of $\tau$-GDA is
determined by

$\displaystyle\xi_{1}(x,y)=-\gamma\mbox{grad}_{x}f(x,y)),\quad\xi_{2}(x,y)=\tau\gamma\mbox{grad}_{y}f(x,y)).$
(11)

Note that $\gamma>0$ and we use a ratio $\tau>0$ to adjust the learning rate
(step size) of the Riemannian gradients. The deterministic $\tau$-GDA can be
readily extended to stochastic $\tau$-GDA by using an unbiased estimation of
the Riemannian gradients [22, 20].

##### Local convergence of deterministic $\tau$-GDA

Based on Theorem 3.1, we obtain the following local convergence of $\tau$-GDA
to DSE on Riemannian manifold. We rewrite the Riemannian Hessian and cross-
gradient in the following linear transform on
$T_{x^{\ast}}\mathcal{M}_{1}\times T_{y^{\ast}}\mathcal{M}_{2}$ using
$(\mathbf{A},\mathbf{B},\mathbf{C})$, i.e.

$\left(\begin{array}[]{ll}-\mathbf{C}&-\mathbf{B}\\\
\tau\mathbf{B}^{\intercal}&-\tau\mathbf{A}\end{array}\right)=\left(\begin{array}[]{ll}-\mbox{Hess}_{x}(x^{\ast},y^{\ast})&-\mbox{grad}^{2}_{yx}f(x^{\ast},y^{\ast})\\\
\tau\mbox{grad}^{2}_{xy}f(x^{\ast},y^{\ast})&\tau\mbox{Hess}_{y}(x^{\ast},y^{\ast})\end{array}\right).$
(12)

######  Theorem 3.2.

Assume $(x^{\ast},y^{\ast})$ is a DSE of $f$. If
$\tau>\frac{\max_{k}|\lambda_{k}(\mathbf{C})|}{\lambda_{min}(\mathbf{A})}$,
then there exists $\gamma^{\bullet}>0$ such that for any
$0<\gamma<\gamma^{\bullet}$, the $\tau$-GDA (defined by (11)) is locally
convergent to $(x^{\ast},y^{\ast})$.

The proof is given Appendix G. When $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$
are Euclidean, the range
$\\{\tau\in\mathbb{R}_{+}|\tau>\|\partial_{xx}^{2}f(x^{\ast},y^{\ast})\|/\lambda_{min}(-\partial_{yy}^{2}f(x^{\ast},y^{\ast}))\\}$
in Theorem 3.2 is broader compared to the range given in [24, Theorem 4.2].
Moreover, this range is sharp as one can construct a counter-example [24,
Theorem 4.1] to show that if $\tau$ is outside this range, the spectral radius
of the Jacobian matrix (9) is strictly larger than one for any $\gamma>0$.

We next analyze the local convergence of $\tau$-GDA to DNE. Theorem 3.2 can be
readily applied. However, we can obtain a broader range of $\tau$. In
Euclidean space, the local convergence to DNE in analyzed in [10, Theorem
1.10] when $\tau=1$, and generalized to $\tau>0$ in [21, Proposition 26].

######  Theorem 3.3.

Assume $(x^{\ast},y^{\ast})$ is a DNE of $f$. If $\tau>0$, then there exists
$\gamma^{\bullet}>0$ such that for any $0<\gamma<\gamma^{\bullet}$, the
$\tau$-GDA (defined by (11)) is locally convergent to $(x^{\ast},y^{\ast})$.

The proof is given in Appendix H. Different to the proof of Theorem 3.2, we
analyze the spectral radius of the Jacobian matrix (9) based on the Ky Fan
inequality [39, Lemma 4.20] as in [10].

###  3.3 Symplectic gradient-adjustment method ($\tau$-SGA)

In Euclidean space, SGA modifies a vector field $\xi$ in an update rule
$\mathbf{T}(x,y)=(x,y)+\gamma\xi(x,y)$ so that it allows one to cross the curl
of $\xi$ which would otherwise results in a strong rotational dynamics near a
stationary point [14]. We apply this idea to $\tau$-GDA on Riemannian manifold
and analyze its local convergence to DSE when $\tau$ is large.

#####  $\tau$-SGA algorithm

SGA adjusts the vector field $\xi$ using an orthogonal vector field which is
constructed from the anti-symmetric part of the Jacobian matrix of $\xi$. In
$\tau$-GDA, $\xi(x,y)=(-\delta(x,y),\tau\eta(x,y))$ where
$\delta(x,y)=\mbox{grad}_{x}f(x,y)$ and $\eta(x,y)=\mbox{grad}_{y}f(x,y)$. The
adjustment of $\tau$-SGA is based on the Riemannian cross-gradients
$\tilde{\mathbf{B}}(x,y)=\mbox{grad}^{2}_{yx}f(x,y)$ and
$\tilde{\mathbf{B}}^{\intercal}(x,y)=\mbox{grad}^{2}_{xy}f(x,y)$, which
results in the following vector fields in the update rule,

$\displaystyle\xi_{1}(x,y)$
$\displaystyle=-\gamma\left(\delta(x,y)+\mu\frac{(\tau+1)\tau}{2}\tilde{\mathbf{B}}(x,y)[\eta(x,y)]\right),$
(13) $\displaystyle\xi_{2}(x,y)$
$\displaystyle=\gamma\left(\tau\eta(x,y)-\mu\frac{\tau+1}{2}\tilde{\mathbf{B}}^{\intercal}(x,y)[\delta(x,y)]\right),$
(14)

where $\mu\in\mathbb{R}$ is a hyper-parameter. In Appendix I, we detail the
derivation of the update rule, and verify that the adjusted direction
$(-\tau\tilde{\mathbf{B}}[\eta],-\tilde{\mathbf{B}}^{\intercal}[\delta])$ is
orthogonal to $\xi=(-\delta,\tau\eta)$ at each point $(x,y)$. In Appendix J,
we discuss how to perform (stochastic) gradient adjustment using auto-
differentiation when $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$ are Euclidean
embedded sub-manifolds [1, Chapter 3.3].

##### Local convergence of deterministic $\tau$-SGA: asymptotic analysis

We study the local convergence of the deterministic $\tau$-SGA to DSE in an
asymptotic regime where $\tau\to\infty$. This simplifies the update rule of
$\tau$-SGA in order to adapt the proof of Theorem 3.2. In this regime, to make
the term $\tilde{\mathbf{B}}[\eta]$ comparable to the term $\delta$ in (13),
we set $\mu\sim\frac{1}{\tau^{2}}$. As a consequence,
$\mu\frac{\tau+1}{2}\sim\frac{1}{\tau}$ and the update rule in (14) can be
simplified to the $\xi_{2}$ in (11). The next theorem makes sense for
$\tau$-SGA when $\tau$ is chosen sufficiently large.

######  Theorem 3.4.

Assume $(x^{\ast},y^{\ast})$ is a DSE of $f$. If
$\tau>\min\left(\frac{\max_{k}|\lambda_{k}(\mathbf{C})|}{\lambda_{min}(\mathbf{A})},\frac{\max_{k}|\lambda_{k}(\mathbf{C}+\theta\mathbf{B}\mathbf{B}^{\intercal})|}{\lambda_{min}(\mathbf{A})}\right)$,
and $\mu=\theta\frac{2}{\tau(\tau+1)}$ with
$0\leq\theta\leq\frac{1}{\lambda_{max}(\mathbf{A})}$, then there exists
$\gamma^{\bullet}>0$ such that for any $0<\gamma<\gamma^{\bullet}$, the
asymptotic $\tau$-SGA ($\xi_{1}$ defined by (13), $\xi_{2}$ defined by (11))
is locally convergent to $(x^{\ast},y^{\ast})$.

The proof is given in Appendix K. Contrary to Theorem 3.2, the valid range of
$\tau$ depends also on the choice of $\theta$. If the spectral radius of
$\mathbf{C}$ is larger than that of
$\mathbf{C}+\theta\mathbf{B}\mathbf{B}^{\intercal}$ with a suitable choice of
$\theta$, then a broader range of $\tau$ could be used in $\tau$-SGA compared
to $\tau$-GDA. In Figure 2, we compare the local convergence rate of
$\tau$-GDA and $\tau$-SGA in the Example 2 of Section 2.2, where
$A=[1;1;1]\in\mathbb{R}^{3\times 1}$, $b=[1;1;0.99]\in\mathbb{R}^{3}$ and
$\kappa=0.1$. In this case, we find numerically that $x^{\ast}\approx 0.9975$.
The initial point of each algorithm is $x=A^{+}b$,$y=(Ax-b)/\|Ax-b\|$. In
Figure 2(a), we verify the range of $\tau$ for $\tau$-GDA convergence.
According to Theorem 3.2 and the discussion after Proposition 2.2, a valid
range of $\tau$ should be larger than
$\frac{\max_{k}|\lambda_{k}(\mathbf{C})|}{\lambda_{min}(\mathbf{A})}=\kappa\|A^{\intercal}A\|/\|Ax^{\ast}-b\|\approx
36.18$. Figure 2(a) shows that when $\tau=10$, $\tau$-GDA can diverge. When
$\tau=50$, $\tau$-GDA converges slowly to the value
$\max_{y}f(x^{\ast},y)=-0.141$. This is a convergence dilemma of $\tau$-GDA:
to achieve local convergence, $\tau$ needs to be large, but the rate can be
very slow. In Figure 2(b) and (c), we study $\tau$-SGA using the same learning
rate $\tau\gamma$ for $y$. Figure 2(b) validates that the gap between the
asymptotic $\tau$-SGA and $\tau$-SGA is not so big even if $\tau=10$.
Therefore Theorem 3.4 is valid for $\tau$-SGA. When $\theta=0.15$, we find
that
$\frac{\max_{k}|\lambda_{k}(\mathbf{C}+\theta\mathbf{B}\mathbf{B}^{\intercal})|)|}{\lambda_{min}(\mathbf{A})}\approx
16.7$, which suggests that $\tau$-SGA can be locally convergent with a smaller
$\tau$ compared to $\tau$-GDA. Figure 2(c) confirms this analysis and it also
shows that a faster rate can be achieved by $\tau$-SGA. Surprisingly,
$\tau$-SGA quickly converges towards $-0.141$ even if $\tau=10$.

(a)

(b)

(c)

Figure 2:  Evolution of $f(x(t),y(t))$ as a function of iteration $t$ in
deterministic $\tau$-GDA and $\tau$-SGA on Example 2
($\gamma=0.001/\tau,\theta=0.15$). a). $\tau$-GDA with $\tau=10$ vs.
$\tau=50$. b). $\tau$-SGA vs. asymptotic $\tau$-SGA with $\tau=10$ and
$\tau=50$. c). $\tau$-GDA at $\tau=50$ vs. $\tau$-SGA at
$\tau\in\\{10,30,50\\}$.

The local convergence of $\tau$-SGA to DNE in Euclidean space is analyzed in
[23, Theorem 10] when $\tau=1$. However, we could not easily extend this
result to $\tau>0$, since when $\tau\neq 1$, a DNE is not necessarily a stable
fixed point of the vector field $\xi=(-\delta,\tau\eta)$, c.f [23, Definition
4]. It is thus unclear if there is a valid range of $\mu$ such that for any
$\tau>0$, $\tau$-SGA is locally convergent to DNE.

##  4 Application to Wasserstein GAN

Based on the local convergence analysis of $\tau$-GDA and $\tau$-SGA, we apply
these algorithms in the stochastic setting to train Wasserstein GAN. In the
min-max problem of GAN, we assume $\mathcal{M}_{1}$ is Euclidean and
$\mathcal{M}_{2}$ is Riemannian. The manifold constraint on $y$ makes the GAN
discriminator Lipschitz-continuous, so that $\max_{y\in\mathcal{M}_{2}}f(x,y)$
gives a lower-bound estimation of the Wasserstein-1 distance (up to
multiplicative constant) between the data and the GAN generator distribution
[5, 30].

###  4.1 Basic setup of Wasserstein GAN

In Wasserstein GAN111 We assume that the data is generated by the random
variable $\phi_{data}\in\mathbb{R}^{d}$. The idea of GAN is to approximate the
probability distribution of $\phi_{data}$ using a generator $G_{x}$
parameterized by $x\in\mathcal{M}_{1}$. It defines a generative model
$\phi_{x}=G_{x}(Z)$, which transform a random variable $Z\in\mathbb{R}^{p}$ to
$\phi_{x}\in\mathbb{R}^{d}$. The generator is optimized so that the
distribution of $\phi_{x}$ approximates that of $\phi_{data}$. This is
achieved through the feedback of a discriminator function $D_{y}$,
parameterized by $y\in\mathcal{M}_{2}$., we are interested in the following
min-max problem

$f(x,y)=\mathbb{E}(D_{y}(\phi_{data}))-\mathbb{E}(D_{y}(\phi_{x})),$

where the discriminator $D_{y}:\mathbb{R}^{d}\to\mathbb{R}$ is Lipschitz-
continuous. The expectations in $f$ are taken with respect to the random
variable $\phi_{data}$ and $\phi_{x}$.222In the stochastic setting, the
expectation of $D_{y}(\phi_{data})$ (resp. $D_{y}(\phi_{x})$) is estimated at
each iteration using a batch of samples of data in the training set (resp. a
batch of samples of $Z$). The number of samples in a batch is called batch
size. We consider modeling two examples of $\phi_{data}$ using low-dimensional
generator (dimension $p<d$). More details are given in Appendix L.

##### Gaussian distribution by linear generator

We want to model $\phi_{data}\sim\mathcal{N}(0,\Sigma)$ on $d=5$. We focus on
$p=4$ to build a PCA-like model. The covariance matrix $\Sigma$ is diagonal
with a small eigenvalue, i.e. $(1,2^{2},3^{2},4^{2},0.01)$. We choose the
generator $G_{x}(Z)=A_{x}Z$, where $A_{x}\in\mathbb{R}^{d\times p}$ and $Z$ is
isotropic normal. For the discriminator, we denote the Stiefel manifold of $k$
orthogonal vectors on $\mathbb{R}^{n}$ by $St(k,n)$, and choose (using a
smooth non-linearity $\rho$, c.f. Appendix L.1)

$D_{y}(\phi)=\langle v_{y},\rho(W_{y}\phi)\rangle,$

where $W_{y}\in St(k,d)$, $v_{y}\in St(1,k)$. Since $A_{x}$ is the only
generator parameter, $\mathcal{M}_{1}=\mathbb{R}^{d\times p}$. The
discriminator parameter $y:=(W_{y},v_{y})$, i.e.
$\mathcal{M}_{2}=St(k,d)\times St(1,k)$.

##### MNIST (digit 0) by DCGAN generator

We aim to model the digit 0 in the MNIST dataset ($d=28\times 28$) using the
DCGAN architecture ($p=128$) in WGAN-GP [16]. The space $\mathcal{M}_{1}$ is
Euclidean with dimension $d_{1}=1556673$. Unlike the CNN discriminator in
WGAN-GP, we consider a hybrid Scattering CNN discriminator which has a smaller
number of trainable parameters [32],

$D_{y}(\phi)=\langle v_{y},\rho(w_{y}\star P(\phi)+b_{y})\rangle.$

The scattering transform $P:\mathbb{R}^{d}\to\mathbb{R}^{I\times n\times n}$
is Lipschitz-continuous [8] and it has no trainable parameter. The scattering
features $P(\phi)\in\mathbb{R}^{I\times n\times n}$ are computed from an image
$\phi$ of size $\sqrt{d}\times\sqrt{d}$. We then apply an orthogonal
convolutional layer to $P(\phi)$ [9]. It has the kernel orthogonality [2], by
reshaping $w_{y}$ (size $J\times I\times k\times k$) into a matrix
$W_{y}\in\mathbb{R}^{J\times(Ik^{2})}$ and imposing that $W_{y}\in
St(J,Ik^{2})$. The bias parameter $b_{y}\in\mathbb{R}^{J}$ and the output of
the convolutional layer is reshaped into a vector in $\mathbb{R}^{JN^{2}}$. As
in the Gaussian case, we further assume $v_{y}\in St(1,JN^{2})$. In summary,
$y:=(W_{y},b_{y},v_{y})$, i.e.
$\mathcal{M}_{2}=St(J,Ik^{2})\times\mathbb{R}^{J}\times St(1,JN^{2})$.

###  4.2 Results on Gaussian distribution by linear generator

To study the local convergence, we first obtain a reasonable solution $(x,y)$
in terms of model quality (to be specified) by solving the min-max problem
using $\tau$-GDA. We then study the local convergence of $\tau$-GDA
initialized from this solution. In Figure 3(a), we show how the function value
$f(x(t),y(t))$ changes over iteration. It is estimated using a batch of $1000$
training samples. For each choice of $\tau\in\\{1,10,100\\}$, we set
$\gamma=0.02/\tau$. This ensures a fixed learning rate $0.02$ for the
discriminator. We observe that when $\tau=1$, $f(x(t),y(t))$ has a huge
oscillation, but when $\tau=10$ and $100$, it has a much smaller oscillation
amplitude. To investigate the underlying reason, we compute in Figure 3(b),
the evolution of an “angle” quantity every one thousand iterations. The angle
at $(x,y)$ is defined as $\frac{\langle
v_{y},\delta(x,y)\rangle}{\|\delta(x,y)\|}$, where
$\delta(x,y)=\mathbb{E}(\rho(W_{y}\phi_{data}))-\mathbb{E}(\rho(W_{y}\phi_{x}))$
is estimated from the same batch of samples during the training. When the
angle is close to $1$, it indicates that $v_{y}$ is solved to be optimal. We
observe that when $\tau=1$, the angle oscillates between positive and negative
values, suggesting that $v_{y(t)}$ is detached from $\delta(x(t),y(t))$.
Therefore the minimization of $f$ over $x$ does not get a good feedback to
improve the model. This is confirmed in Figure 3(c), where we compute the EMD
distance [34] every one thousand iterations, between the empirical measure of
$\phi_{data}$ and $\phi_{x(t)}$ (using $2000$ validation samples). We observe
that the EMD distance increases over $t$ when $\tau=1$. On the contrary, it
stays close to one when $\tau=10$ and $100$. When $\tau=100$, we find that the
covariance error $\|\Sigma-A_{x(t)}A_{x(t)}^{\intercal}\|$ stays around $0.3$
over all iterations. This error is between the smallest eigenvalue of $\Sigma$
(which is $0.01$, the PCA error) and the second smallest eigenvalue of
$\Sigma$. Therefore to use a large enough $\tau$ is crucial to avoid large
model error, as it is necessary to ensure a positive angle. This phenomenon is
consistent with the local convergence of $\tau$-GDA to DSE in Figure 2(a).

(a) $f(x(t),y(t))$

(b) angle$(x(t),y(t))$

(c) EMD$(x(t))$

Figure 3:  Evolution of $f(x(t),y(t))$, angle$(x(t),y(t))$, and EMD$(x(t))$
distance over iteration $t$ in stochastic $\tau$-GDA on Wasserstein GAN of
Gaussian distribution. Top: $\tau=100,\gamma=0.0002$. Middle:
$\tau=10,\gamma=0.002$, Bottom: $\tau=1,\gamma=0.02$.  Table 2:  Last
iteration measures in stochastic $\tau$-GDA and $\tau$-SGA on Wasserstein GAN
of MNIST (digit 0). We report the $f$, angle, and FID scores computed at
$(x(T),y(T))$. $\tau$-GDA is trained for $T=2\times 10^{5}$ iterations with
$\gamma=0.05/\tau$. $\tau$-SGA is trained longer with $\gamma=0.01$, $\tau=5$,
$\theta=0.075$. Only significant digits are reported with respect to the
standard deviation (shown in $(\pm)$) of the statistical estimation of these
quantities (see footnote 3 and 4).

|  $\tau$-GDA |
---|---|---
$\tau$ | $f$ | angle | FID (train) | FID (val)
$5$ | -0.09 | -0.26 | 21 |  29 ($\pm 0.9$)
$10$ | 0.09 | 0.11 | 34 |  36 ($\pm 0.6$)
$30$ | 0.036 | 0.4 | 13 |  20 ($\pm 0.95$)

|  $\tau$-SGA |
---|---|---
$T(10^{6})$ | $f$ | angle | FID (train) | FID (val)
$1$ | 0.02 | 0.16 | 12 |  21 ($\pm 0.9$)
$2$ | 0.018 | 0.2 |  9.3 ($\pm 0.2$) |  15 ($\pm 0.5$)
$3$ | 0.016 | 0.2 |  8.2 ($\pm 0.2$) |  14 ($\pm 0.3$)

###  4.3 Results on MNIST (digit 0) by DCGAN generator

We study the convergence of $\tau$-GDA and $\tau$-SGA to improve the training
of a pre-trained model obtained by alternating $\tau$-GDA method (c.f.
Appendix L.5). To measure the quality of the pre-trained model, we compute the
FID scores [36]. It is estimated on both training and validation samples of
MNIST.333 To compute FID (train) (resp. FID (val)), we generate the same
number of fake samples as the training samples (resp. as the validation
samples), and report an average value by re-sampling the fake samples $10$
times.  We obtain a pre-trained model with FID (train) $12$ and FID (val)
$18$. The results of $\tau$-GDA and $\tau$-SGA using batch size $128$,
initialized from this pre-trained model, are reported in Table 2. For
$\tau$-GDA, we choose $\gamma=0.05/\tau$ among $\tau\in\\{5,10,30\\}$. We
observe that when $\tau=30$, $\tau$-GDA has a more stable dynamics than
$\tau=5$ and $10$ because the angle stays around a positive constant.444To
compute the angle, we estimate $\delta(x,y)=\mathbb{E}(\rho(w_{y}\star
P(\phi_{data})+b_{y}))-\mathbb{E}(\rho(w_{y}\star P(\phi_{x})+b_{y}))$. We
compute the $f$ and angle values by re-sampling $\phi_{x}$ ten times (similar
to the estimation of FID). When $\tau=5$ and $\tau=10$, sometimes we observe a
positive angle at $t=T$ in $\tau$-GDA because it oscillates around zero over
$t$.  On the other hand, a larger $\tau$ tends to slowdown the reduction of
$f$ as shown in Figure 2(a). Due to the stochastic gradients, we observe a
stagnation of $f$ over iteration when $\tau=30$. Indeed, the FID scores at
$T=3\times 10^{6}$ are similar to those at $T=2\times 10^{5}$. Facing such a
convergence dilemma, we evaluate the performance of $\tau$-SGA using $\tau=5$
such that $\tau\gamma$ remains the same. This is the case where $\tau$-GDA is
not convergent due to oscillating angles. We find that both the angle and the
FID scores are significantly improved with a suitable choice of $\theta$ and
$T$ in $\tau$-SGA.

##  5 Conclusion

In this article, we extend the notions of DSE and DNE to Riemannian manifold
from Euclidean space, and then study the local convergence of $\tau$-GDA and
$\tau$-SGA in theory. Our theoretical results are based on the spectral
analysis of an intrinsic Jacobian matrix analyzed in local coordinate. This
methodology can be served as a basis to analyze other algorithms. However,
this local analysis approach may not be suitable to analyze stochastic
algorithms since they typically have a more global behavior. For this reason,
we study stochastic $\tau$-GDA and $\tau$-SGA numerically, where we observe
certain phenomenon similar to our theoretical results, such as the importance
of the choice of a suitable $\tau$ to ensure the local convergence to DSE.
From this insight, we are able to improve the convergence of stochastic
$\tau$-GDA in Wasserstein GANs, leading to a better generative model.

## References

  * Absil et al. [2008] P.-A. Absil, R Mahony, and Rodolphe Sepulchre.  _Optimization Algorithms on Matrix Manifolds_.  Princeton University Press, Princeton, 2008.
  * Achour et al. [2022] El Mehdi Achour, François Malgouyres, and Franck Mamalet.  Existence, stability and scalability of orthogonal convolutional neural networks.  _Journal of Machine Learning Research_ , 23:1–56, jan 2022\.
  * Azizian et al. [2020] Waïss Azizian, Ioannis Mitliagkas, Simon Lacoste-Julien, and Gauthier Gidel.  A Tight and Unified Analysis of Gradient-Based Methods for a Whole Spectrum of Differentiable Games.  In _International Conference on Artificial Intelligence and Statistics_ , volume 108 of _Proceedings of Machine Learning Research_ , pages 2863–2873, Virtual Event, 2020.
  * Balduzzi et al. [2018] David Balduzzi, Sebastien Racaniere, James Martens, Jakob Foerster, Karl Tuyls, and Thore Graepel.  The Mechanics of n-Player Differentiable Games.  In _International Conference on Machine Learning_ , volume 80, pages 354–363, Stockholmsmässan, Stockholm Sweden, 2018.
  * Biau et al. [2021] Gérard Biau, Maxime Sangnier, and Ugo Tanielian.  Some theoretical insights into Wasserstein GANs.  _Journal of Machine Learning Research_ , 22:1–45, 2021\.
  * Biswas et al. [2022] Koushik Biswas, Sandeep Kumar, Shilpak Banerjee, and Ashish Kumar Pandey.  Smooth Maximum Unit: Smooth Activation Function for Deep Networks using Smoothing Maximum Technique.  In _Conference on Computer Vision and Pattern Recognition_ , pages 784–793, New Orleans, LA, USA, 2022.
  * Boumal [2023] Nicolas Boumal.  _An Introduction to Optimization on Smooth Manifolds_.  Cambridge University Press, 2023.
  * Bruna and Mallat [2013] J Bruna and S Mallat.  Invariant Scattering Convolution Networks.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 35(8):1872–1886, 2013.
  * Cisse et al. [2017] Moustapha Cisse, Piotr Bojanowski, Edouard Grave, Yann Dauphin, and Nicolas Usunier.  Parseval networks: improving robustness to adversarial examples.  In _International Conference on Machine Learning_ , pages 854–863, Sydney, Australia, 2017.
  * Daskalakis and Panageas [2018] Constantinos Daskalakis and Ioannis Panageas.  The Limit Points of (Optimistic) Gradient Descent in Min-Max Optimization.  In _Advances in neural information processing systems_ , pages 9256–9266, Red Hook, NY, USA, 2018.
  * de Montbrun and Renault [2022] Étienne de Montbrun and Jérôme Renault.  Optimistic Gradient Descent Ascent in Zero-Sum and General-Sum Bilinear Games.  _arXiv preprint arXiv:2208.03085_ , 2022.
  * Fiez and Ratliff [2021] Tanner Fiez and Lillian J Ratliff.  Local convergence analysis of gradient descent ascent with finite timescale separation.  In _International Conference on Learning Representation_ , 2021.
  * Fiez et al. [2020] Tanner Fiez, Benjamin Chasnov, and Lillian Ratliff.  Implicit learning dynamics in stackelberg games: Equilibria characterization, convergence analysis, and empirical study.  In _International Conference on Machine Learning_ , pages 3133–3144, Virtual Event, 2020. PMLR.
  * Gemp and Mahadevan [2018] Ian Gemp and Sridhar Mahadevan.  Global Convergence to the Equilibrium of GANs using Variational Inequalities.  _arXiv preprint arXiv:1808.01531_ , 2018.
  * Gidel et al. [2019] Gauthier Gidel, Hugo Berard, Gaëtan Vignoud, Pascal Vincent, and Simon Lacoste-Julien.  A Variational Inequality Perspective on Generative Adversarial Networks.  In _International Conference on Learning Representations_ , New Orleans, LA, USA, 2019.
  * Gulrajani et al. [2017] Ishaan Gulrajani, Faruk Ahmed, Martin Arjovsky, Vincent Dumoulin, and Aaron C Courville.  Improved training of wasserstein gans.  In _Advances in neural information processing systems_ , volume 30, pages 5767–5777, Long Beach, CA, USA, 2017.
  * Han et al. [2023] Andi Han, Bamdev Mishra, Pratik Jawanpuria, Pawan Kumar, and Junbin Gao.  Riemannian Hamiltonian methods for min-max optimization on manifolds.  _SIAM Journal on Optimization_ , 33(3):1797–1827, 2023.
  * Hsieh et al. [2021] Ya-Ping Hsieh, Panayotis Mertikopoulos, and Volkan Cevher.  The Limits of Min-Max Optimization Algorithms: Convergence to Spurious Non-Critical Sets.  In _International Conference on Machine Learning_ , volume 139, pages 4337–4348, Virtual Event, 2021.
  * Hu et al. [2023] Zihao Hu, Guanghui Wang, Xi Wang, Andre Wibisono, Jacob Abernethy, and Molei Tao.  Extragradient Type Methods for Riemannian Variational Inequality Problems.  _arXiv preprint arXiv:2309.14155_ , 2023.
  * Huang and Gao [2023] Feihu Huang and Shangqian Gao.  Gradient Descent Ascent for Minimax Problems on Riemannian Manifolds.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 1(7):8466 – 8476, 2023.
  * Jin et al. [2020] Chi Jin, Praneeth Netrapalli, and Michael Jordan.  What is local optimality in nonconvex-nonconcave minimax optimization?  In _International Conference on Machine Learning_ , pages 4880–4889, Virtual Event, 2020.
  * Jordan et al. [2022] Michael Jordan, Tianyi Lin, and Emmanouil-Vasileios Vlatakis-Gkaragkounis.  First-Order Algorithms for Min-Max Optimization in Geodesic Metric Spaces.  In _Advances in Neural Information Processing Systems_ , volume 35, pages 6557–6574, New Orleans, Louisiana, USA, 2022.
  * Letcher et al. [2019] Alistair Letcher, David Balduzzi, Sébastien Racanière, James Martens, Jakob Foerster, Karl Tuyls, and Thore Graepel.  Differentiable Game Mechanics.  _Journal of Machine Learning Research_ , 20(84):1–40, 2019.
  * Li et al. [2022] Haochuan Li, Farzan Farnia, Subhro Das, and Ali Jadbabaie.  On convergence of gradient descent ascent: A tight local analysis.  In _International Conference on Machine Learning_ , pages 12717–12740, Baltimore, Maryland USA, 2022. PMLR.
  * Lin et al. [2020] Tianyi Lin, Chi Jin, and Michael I Jordan.  On Gradient Descent Ascent for Nonconvex-Concave Minimax Problems.  In _International Conference on Machine Learning_ , pages 6083–6093, Virtual Event, 2020.
  * Loizou et al. [2020] Nicolas Loizou, Hugo Berard, Alexia Jolicoeur-Martineau, Pascal Vincent, Simon Lacoste-Julien, and Ioannis Mitliagkas.  Stochastic Hamiltonian Gradient Methods for Smooth Games.  In _International Conference on Machine Learning_ , volume 119, pages 6370–6381, Virtual Event, 2020.
  * Mahdavinia et al. [2022] Pouria Mahdavinia, Yuyang Deng, Haochuan Li, and Mehrdad Mahdavi.  Tight analysis of extra-gradient and optimistic gradient methods for nonconvex minimax problems.  _Advances in Neural Information Processing Systems_ , 35:31213–31225, 2022.
  * Manton [2002] Jonathan H Manton.  Optimization algorithms exploiting unitary constraints.  _IEEE Transactions on Signal Processing_ , 50(3):635–650, 2002.
  * Mertikopoulos et al. [2019] Panayotis Mertikopoulos, Bruno Lecouat, Houssam Zenati, Chuan-Sheng Foo, Vijay Chandrasekhar, and Georgios Piliouras.  Optimistic mirror descent in saddle-point problems: Going the extra (gradient) mile.  In _International Conference on Learning Representations_ , New Orleans, LA, USA, 2019.
  * Müller et al. [2019] Jan Müller, Reinhard Klein, and Michael Weinmann.  Orthogonal wasserstein gans.  _arXiv preprint arXiv:1911.13060_ , 2019.
  * Ortega and Rheinboldt [1970] J M Ortega and W C Rheinboldt.  _Iterative Solution of Nonlinear Equations in Several Variables_.  New York : Academic Press, 1970.
  * Oyallon et al. [2017] Edouard Oyallon, Eugene Belilovsky, and Sergey Zagoruyko.  Scaling the Scattering Transform: Deep Hybrid Networks.  In _International Conference on Computer Vision_ , Venice, Italy, 2017\.
  * Razaviyayn et al. [2020] Meisam Razaviyayn, Tianjian Huang, Songtao Lu, Maher Nouiehed, Maziar Sanjabi, and Mingyi Hong.  Nonconvex min-max optimization: Applications, challenges, and recent theoretical advances.  _IEEE Signal Processing Magazine_ , 37(5):55–66, 2020.
  * Rubner et al. [2000] Yossi Rubner, Carlo Tomasi, and Leonidas J Guibas.  The earth mover’s distance as a metric for image retrieval.  _International journal of computer vision_ , 40:99–121, 2000.
  * Rudin [1976] W Rudin.  _Principles of Mathematical Analysis_.  McGraw-Hill Publishing Company, 3rd edition, 1976.
  * Seitzer [2020] Maximilian Seitzer.  pytorch-fid: FID Score for PyTorch.  https://github.com/mseitzer/pytorch-fid, 2020.
  * Wang et al. [2023] Xi Wang, Deming Yuan, Yiguang Hong, Zihao Hu, Lei Wang, and Guodong Shi.  Riemannian Optimistic Algorithms.  _arXiv preprint arXiv:2308.16004_ , 2023.
  * Wang et al. [2020] Yuanhao Wang, Guodong Zhang, and Jimmy Ba.  On Solving Minimax Optimization Locally: A Follow-the-Ridge Approach.  In _International Conference on Learning Representations_ , Addis Ababa, Ethiopia, 2020.
  * Zhan [2002] Xingzhi Zhan.  _Matrix Inequalities_.  Springer Berlin, Heidelberg, 2002.
  * Zhang et al. [2022] Guodong Zhang, Yuanhao Wang, Laurent Lessard, and Roger B Grosse.  Near-optimal local convergence of alternating gradient descent-ascent for minimax optimization.  In _International Conference on Artificial Intelligence and Statistics_ , pages 7659–7679, Virtual Event, 2022.
  * Zhang et al. [2023] Peiyuan Zhang, Jingzhao Zhang, and Suvrit Sra.  Sion’s Minimax Theorem in Geodesic Metric Spaces and a Riemannian Extragradient Algorithm.  _SIAM Journal on Optimization_ , 33(4), 2023.

##  Appendix A Implicit function theorem on Riemannian manifold

We state the implicit function theorem in the manifold case which allows one
to understand why $(x^{\ast},y^{\ast})$ is a local minimax point. Similar to
the local definition of DSE, this theorem implies the existence of a solution
$(x,h(x))\approx(x^{\ast},y^{\ast})$ s.t. $\mbox{grad}_{y}f(x,h(x))=0$. As a
consequence, $h(x)$ is a strict local maximum of $y\mapsto f(x,y)$, thanks to
the continuity of $\mbox{Hess}_{y}f$ on $\mathcal{M}_{1}\times\mathcal{M}_{2}$
and the continuity of $h$ near $x^{\ast}$ (i.e. provided that $x$ is close
enough to $x^{\ast}$, we have the positive definiteness of
$-\mbox{Hess}_{y}f(x,h(x))$).

Denote the set of continuously differentiable vector fields on
$\mathcal{M}_{2}$ by $\mathcal{X}(\mathcal{M}_{2})$. Consider a parameterized
vector field $\xi$ defined at each $x\in\mathcal{M}_{1}$ with
$\xi(x,\cdot)\in\mathcal{X}(\mathcal{M}_{2})$. Denote the Riemannian
connection on $\mathcal{M}_{2}$ by $\nabla_{y}$. By definition, for each
$(x,y)\in\mathcal{M}_{1}\times\mathcal{M}_{2}$, $\xi(x,y)\in
T_{y}\mathcal{M}_{2}$, and $\nabla_{y}\xi(x,y)$ is a linear map from
$T_{y}\mathcal{M}_{2}$ to $T_{y}\mathcal{M}_{2}$.

######  Theorem A.1.

Assume $\xi$ is continuously differentiable on an open set
$E\subset\mathcal{M}_{1}\times\mathcal{M}_{2}$. Let $(x^{\ast},y^{\ast})\in E$
be a solution of

$\xi(x,y)=0,\quad(x,y)\in E.$

If $\nabla_{y}\xi(x^{\ast},y^{\ast})$ is invertible on
$T_{y^{\ast}}\mathcal{M}_{2}$, then there exists an open set $U_{1}\times
U_{2}\subset E$ and an open set $W_{1}\subset U_{1}$ such that

$\forall x\in W_{1},\exists!\,y\in U_{2}\quad\mbox{s.t.}\quad\xi(x,y)=0.$

Let the unique $y=h(x)$, i.e. $h:W_{1}\to U_{2}$, such that

$y^{\ast}=h(x^{\ast}),\quad\xi(x,h(x))=0,\quad\forall x\in W_{1},$

then $h$ is continuously differentiable on $W_{1}$. Let $\delta\in
T_{x^{\ast}}\mathcal{M}_{1}$, and denote the tangent map of $h$ at $x^{\ast}$
by $D_{x}h(x^{\ast})$, we have

$D_{x}h(x^{\ast})[\delta]=-\nabla_{y}\xi(x^{\ast},y^{\ast})^{-1}\cdot
D_{x}\xi(x^{\ast},y^{\ast})[\delta].$

The proof is omitted due to its length. The main idea of the proof is similar
to that of Theorem 2.1, by using a local coordinate chart around the point
$(x^{\ast},y^{\ast})$ to represent a smooth vector field. We then apply the
proof technique of the implicit function theorem in Euclidean space [35,
Theorem 9.28].

##  Appendix B Proof of Theorem 2.1

Let $\mathcal{X}(\mathcal{M}_{1})$ (resp. $\mathcal{X}(\mathcal{M}_{2})$) be
the set of continuously differentiable vector fields on $\mathcal{M}_{1}$
(resp. $\mathcal{M}_{2}$). Since $f$ is twice continuously differentiable, we
have for any $y\in\mathcal{M}_{2}$ (resp. $x\in\mathcal{M}_{1}$),
$\mbox{grad}_{x}f(\cdot,y)\in\mathcal{X}(\mathcal{M}_{1})$ (resp.
$\mbox{grad}_{y}f(x,\cdot)\in\mathcal{X}(\mathcal{M}_{2})$).

To show that the local definition of DSE is equivalent to its intrinsic
definition, we shall write the Riemannian gradients $\mbox{grad}_{x}f(x,y)$
and $\mbox{grad}_{y}f(x,y)$ in the local coordinate chart $(O_{1}\times
O_{2},\varphi_{1}\times\varphi_{2})$. We write $g_{1}$ and $g_{2}$ to denote
the Riemannian metric on $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$.

For each $x\in O_{1}$ (resp. $y\in O_{2}$), let $\\{E_{1,i}(x)\\}_{i\leq
d_{1}}$ (resp. $\\{E_{2,j}(y)\\}_{j\leq d_{2}}$) be the canonical basis of the
tangent space $T_{x}\mathcal{M}_{1}$ (resp. $T_{y}\mathcal{M}_{2}$), defined
by the tangent map $D\varphi^{-1}_{1}(\varphi_{1}(x))[e_{1,i}]$ of the
canonical basis $\\{e_{1,i}\\}_{i\leq d_{1}}$ on $\mathbb{R}^{d_{1}}$ (resp.
$D\varphi^{-1}_{2}(\varphi_{2}(y))[e_{2,j}]$ of $\\{e_{2,j}\\}_{j\leq d_{2}}$
on $\mathbb{R}^{d_{2}}$).

Let the local coordinate of $(x,y)\in O_{1}\times O_{2}$ be
$(u_{1},u_{2})=(\varphi_{1}(x),\varphi_{2}(y))\in\bar{O}_{1}\times\bar{O}_{2}=\varphi_{1}(O_{1})\times\varphi_{2}(O_{2})\subset\mathbb{R}^{d_{1}}\times\mathbb{R}^{d_{2}}$.
We can represent

$\displaystyle\xi_{1}(x,y)=\mbox{grad}_{x}f(x,y)$ $\displaystyle=\sum_{i\leq
d_{1}}\xi_{1,i}(x,y)E_{1,i}(x),$ (15)
$\displaystyle\xi_{2}(x,y)=\mbox{grad}_{y}f(x,y)$ $\displaystyle=\sum_{j\leq
d_{2}}\xi_{2,j}(x,y)E_{2,j}(y).$ (16)

From (15) and (16), the local coordinates of the Riemannian gradients on
$\bar{O}_{1}\times\bar{O}_{2}$ are

$\displaystyle\bar{\xi}_{1}(u_{1},u_{2})$ $\displaystyle=\sum_{i\leq
d_{1}}\xi_{1,i}(\varphi^{-1}_{1}(u_{1}),\varphi^{-1}_{2}(u_{2}))e_{1,i},$ (17)
$\displaystyle\bar{\xi}_{2}(u_{1},u_{2})$ $\displaystyle=\sum_{j\leq
d_{2}}\xi_{2,j}(\varphi^{-1}_{1}(u_{1}),\varphi^{-1}_{2}(u_{2}))e_{2,j}.$ (18)

##### Equivalence between (1) and (4)

Let the local Riemannian metric matrix at $x=\varphi^{-1}_{1}(u_{1})\in O_{1}$
and $y=\varphi^{-1}_{2}(u_{2})\in O_{2}$ be

$\bar{g}_{1}(u_{1})=(g_{1}(E_{1,i}(x),E_{1,i^{\prime}}(x)))_{i,i^{\prime}\leq
d_{1}},\quad\bar{g}_{2}(u_{2})=(g_{2}(E_{2,j}(y),E_{2,j^{\prime}}(y)))_{j,j^{\prime}\leq
d_{2}}.$ (19)

From [1, chap. 3.6], we have for
$(u_{1},u_{2})\in\bar{O}_{1}\times\bar{O}_{2}$,

$\displaystyle\bar{\xi}_{1}(u_{1},u_{2})=\bar{g}_{1}(u_{1})^{-1}\cdot\partial_{u_{1}}\bar{f}(u_{1},u_{2}),$
(20)
$\displaystyle\bar{\xi}_{2}(u_{1},u_{2})=\bar{g}_{2}(u_{2})^{-1}\cdot\partial_{u_{2}}\bar{f}(u_{1},u_{2}).$
(21)

The equivalence between the DSE condition (1) and (4) follows from (20) and
(21), as $\bar{\xi}_{1}(u_{1}^{\ast},u_{2}^{\ast})=0$ i.f.f. (if and only if)
$\mbox{grad}_{x}f(x^{\ast},y^{\ast})=0$ (similarly for the relationship
between $\bar{\xi}_{2}$ and $\mbox{grad}_{y}f$).

##### Equivalence between (2) and (5)

Let $\eta^{\ast}\in T_{y^{\ast}}\mathcal{M}_{2}$, with its local coordinate
$\bar{\eta}^{\ast}=D\varphi_{2}(y^{\ast})[\eta^{\ast}]$. To show that the
positive definiteness of $-\mbox{Hess}_{y}f(x^{\ast},y^{\ast})$ is equivalent
to the positive definiteness of
$-\partial^{2}_{u_{2}u_{2}}\bar{f}(u_{1}^{\ast},u_{2}^{\ast})$, it is
sufficient to verify that

$\langle\bar{\eta}^{\ast},\partial^{2}_{u_{2}u_{2}}\bar{f}(u_{1}^{\ast},u_{2}^{\ast})\bar{\eta}^{\ast}\rangle=\langle\eta^{\ast},\mbox{Hess}_{y}f(x^{\ast},y^{\ast})[\eta^{\ast}]\rangle_{y^{\ast}},$
(22)

where $\langle\cdot,\cdot\rangle$ denotes the inner product on Euclidean
space, and $\langle\cdot,\cdot\rangle_{y^{\ast}}$ denotes the Riemannian inner
product on $T_{y^{\ast}}M$.

Recall that the Riemannian Hessian on $\mathcal{M}_{2}$ is defined by the
Riemannian connection
$\nabla_{y}:\mathcal{X}(\mathcal{M}_{2})\times\mathcal{X}(\mathcal{M}_{2})\to\mathcal{X}(\mathcal{M}_{2})$
with $\mbox{Hess}_{y}f(x,y)=\nabla_{y}\xi_{2}(x,y)$. By the R-$linearity$ and
Leibniz law of the connection, we have

$\displaystyle\nabla_{y}\xi_{2}(x^{\ast},y^{\ast})[\eta^{\ast}]$
$\displaystyle=\nabla_{y}(\sum_{j}\xi_{2,j}(x^{\ast},y^{\ast})E_{2,j}(y^{\ast}))[\eta^{\ast}]$
$\displaystyle=\sum_{j}\nabla_{y}(\xi_{2,j}(x^{\ast},y^{\ast})E_{2,j}(y^{\ast}))[\eta^{\ast}]$
$\displaystyle=\sum_{j}(D_{y}\xi_{2,j}(x^{\ast},y^{\ast})[\eta^{\ast}])E_{2,j}(y^{\ast})+\xi_{2,j}(x^{\ast},y^{\ast})\nabla_{y}(E_{2,j}(y^{\ast}))[\eta^{\ast}]$
(23) $\displaystyle=\sum_{j}\left(\frac{\partial{\bar{\xi}_{2,j}}}{\partial
u_{2}}(u_{1}^{\ast},u_{2}^{\ast})\bar{\eta}^{\ast}\right)E_{2,j}(y^{\ast}).$
(24)

We have obtained (24) from (23) because any DSE $(x^{\ast},y^{\ast})$ is a
stationary point of $f$, therefore $\xi_{2,j}(x^{\ast},y^{\ast})=0$ for each
$j\leq d_{2}$. Moreover, the tangent map
$D_{y}\xi_{2,j}(x^{\ast},y^{\ast})[\eta^{\ast}]=\frac{\partial{\bar{\xi}_{2,j}}}{\partial
u_{2}}(u_{1}^{\ast},u_{2}^{\ast})\bar{\eta}^{\ast}$ by definition [1, chap.
3.5.4].

Therefore the local coordinate of the tangent vector
$\nabla_{y}\xi_{2}(x^{\ast},y^{\ast})[\eta^{\ast}]\in
T_{y^{\ast}}\mathcal{M}_{2}$ is $\frac{\partial{\bar{\xi}_{2}}}{\partial
u_{2}}(u_{1}^{\ast},u_{2}^{\ast})\bar{\eta}^{\ast}\in\mathbb{R}^{d_{2}}$. It
follows from (21),(24), and the Riemannian inner product [1, chap. 3.6] that

$\displaystyle\langle\eta^{\ast},\mbox{Hess}_{y}f(x^{\ast},y^{\ast})[\eta^{\ast}]\rangle_{y^{\ast}}$
$\displaystyle=\left\langle\bar{\eta}^{\ast},\bar{g}_{2}(u_{2}^{\ast})\cdot\frac{\partial{\bar{\xi}_{2}}}{\partial
u_{2}}(u_{1}^{\ast},u_{2}^{\ast})\bar{\eta}^{\ast}\right\rangle$ (25)
$\displaystyle=\left\langle\bar{\eta}^{\ast},\partial^{2}_{u_{2}u_{2}}\bar{f}(u_{1}^{\ast},u_{2}^{\ast})\bar{\eta}^{\ast}\right\rangle.$

Therefore (22) holds. This verifies the equivalence between (2) and (5).

##### Equivalence between (3) and (6)

Assume $(x^{\ast},y^{\ast})$ is a DSE of $f$, from the equivalence between
(1),(2) and (4),(5), we deduce that $(x^{\ast},y^{\ast})$ (resp.
$(u_{1}^{\ast},u_{2}^{\ast})$) is a stationary point of $f$ (resp. $\bar{f}$).
Moreover, $\mbox{Hess}_{y}f(x^{\ast},y^{\ast})$ and
$\partial^{2}_{u_{2}u_{2}}\bar{f}(u_{1}^{\ast},u_{2}^{\ast})$ are negative
definite.

To simplify the notation, let

$\displaystyle
A=-\partial^{2}_{u_{2}u_{2}}\bar{f}(u_{1}^{\ast},u_{2}^{\ast}),\quad
B=\partial^{2}_{u_{2}u_{1}}\bar{f}(u_{1}^{\ast},u_{2}^{\ast}),\quad
C=\partial^{2}_{u_{1}u_{1}}\bar{f}(u_{1}^{\ast},u_{2}^{\ast}).$ (26)

The condition (3) is equivalent to the positive definiteness of the matrix
$C+BA^{-1}B^{\intercal}$.

Let $\delta^{\ast}\in T_{x^{\ast}}\mathcal{M}_{1}$, with its local coordinate
$\bar{\delta}^{\ast}=D\varphi_{1}(x^{\ast})[\delta^{\ast}]$. As in (22), it
suffices to verify that

$\langle\bar{\delta}^{\ast},(C+BA^{-1}B^{\intercal})\bar{\delta}^{\ast}\rangle=\langle\delta^{\ast},[\mbox{Hess}_{x}f-\mbox{grad}^{2}_{yx}f\cdot(\mbox{Hess}_{y}f)^{-1}\cdot\mbox{grad}^{2}_{xy}f](x^{\ast},y^{\ast})[\delta^{\ast}]\rangle_{x^{\ast}}.$

Based on the proof of (22), we also have
$\langle\bar{\delta}^{\ast},C\bar{\delta}^{\ast}\rangle=\langle\delta^{\ast},\mbox{Hess}_{x}f(x^{\ast},y^{\ast})[\delta^{\ast}]\rangle_{x^{\ast}}$.
It remains to verify that

$\langle\bar{\delta}^{\ast},BA^{-1}B^{\intercal}\bar{\delta}^{\ast}\rangle=\langle\delta^{\ast},[-\mbox{grad}^{2}_{yx}f\cdot(\mbox{Hess}_{y}f)^{-1}\cdot\mbox{grad}^{2}_{xy}f](x^{\ast},y^{\ast})[\delta^{\ast}]\rangle_{x^{\ast}}.$
(27)

Let
$\eta^{\ast}=[\mbox{Hess}_{y}f^{-1}\cdot\mbox{grad}^{2}_{xy}f](x^{\ast},y^{\ast})[\delta^{\ast}]$,
we compute the local coordinate $\bar{\eta}^{\ast}$ of $\eta^{\ast}$ by
converting the following equation

$\mbox{Hess}_{y}f(x^{\ast},y^{\ast})[\eta^{\ast}]=\mbox{grad}^{2}_{xy}f(x^{\ast},y^{\ast})[\delta^{\ast}],$
(28)

into local coordinate. From (24), we have
$\mbox{Hess}_{y}f(x^{\ast},y^{\ast})[\eta^{\ast}]=\sum_{j}\left(\frac{\partial{\bar{\xi}_{2,j}}}{\partial
u_{2}}(u_{1}^{\ast},u_{2}^{\ast})\bar{\eta}^{\ast}\right)E_{2,j}(y^{\ast})$.

By (16) and (18), we obtain

$\displaystyle\mbox{grad}^{2}_{xy}f(x^{\ast},y^{\ast})[\delta^{\ast}]$
$\displaystyle=D_{x}\mbox{grad}_{y}f(x^{\ast},y^{\ast})[\delta^{\ast}]=\sum_{j}D_{x}\xi_{2,j}(x^{\ast},y^{\ast})[\delta^{\ast}]E_{2,j}(y^{\ast})$
$\displaystyle=\sum_{j}\left(\frac{\partial\bar{\xi}_{2,j}}{\partial
u_{1}}(u_{1}^{\ast},u_{2}^{\ast})\bar{\delta}^{\ast}\right)E_{2,j}(y^{\ast}).$
(29)

It follows that (28) is equivalent to

$\displaystyle\frac{\partial{\bar{\xi}_{2}}}{\partial
u_{2}}(u_{1}^{\ast},u_{2}^{\ast})\bar{\eta}^{\ast}=\frac{\partial\bar{\xi}_{2}}{\partial
u_{1}}(u_{1}^{\ast},u_{2}^{\ast})\bar{\delta}^{\ast}.$ (30)

We can compute the right hand side of (27) in the local coordinate as in (29)
and (25),

$\displaystyle\langle\delta^{\ast},-\mbox{grad}^{2}_{yx}f(x^{\ast},y^{\ast})[\eta^{\ast}]\rangle_{x^{\ast}}$
$\displaystyle=-\langle\delta^{\ast},D_{y}\mbox{grad}_{x}f(x^{\ast},y^{\ast})[\eta^{\ast}]\rangle_{x^{\ast}}$
$\displaystyle=-\left\langle\delta^{\ast},\sum_{i}\left(\frac{\partial\bar{\xi}_{1,i}}{\partial
u_{2}}(u_{1}^{\ast},u_{2}^{\ast})\bar{\eta}^{\ast}\right)E_{1,i}(x^{\ast})\right\rangle_{x^{\ast}}$
$\displaystyle=-\left\langle\bar{\delta}^{\ast},\bar{g}_{1}(u_{1}^{\ast})\cdot\frac{\partial\bar{\xi}_{1}}{\partial
u_{2}}(u_{1}^{\ast},u_{2}^{\ast})\bar{\eta}^{\ast}\right\rangle.$

The left hand side of (27) can computed based on (20) and (21), which results
in

$\displaystyle A$
$\displaystyle=-\bar{g}_{2}(u_{2}^{\ast})\cdot\frac{\partial{\bar{\xi}_{2}}}{\partial
u_{2}}(u_{1}^{\ast},u_{2}^{\ast}),$ (31) $\displaystyle B$
$\displaystyle=\bar{g}_{1}(u_{1}^{\ast})\cdot\frac{\partial{\bar{\xi}_{1}}}{\partial
u_{2}}(u_{1}^{\ast},u_{2}^{\ast}),$ (32) $\displaystyle C$
$\displaystyle=\bar{g}_{1}(u_{1}^{\ast})\cdot\frac{\partial{\bar{\xi}_{1}}}{\partial
u_{1}}(u_{1}^{\ast},u_{1}^{\ast}).$

Furthermore, the symmetry of the Hessian matrix of $\bar{f}$ implies that

$B^{\intercal}=\partial^{2}_{u_{1}u_{2}}\bar{f}(u_{1}^{\ast},u_{2}^{\ast})=\bar{g}_{2}(u_{2}^{\ast})\cdot\frac{\partial\bar{\xi}_{2}}{\partial
u_{1}}(u_{1}^{\ast},u_{2}^{\ast}).$ (33)

As a consequence, (30) implies that
$\bar{\eta}^{\ast}=-A^{-1}B^{\intercal}\bar{\delta}^{\ast}$, and

$\langle\bar{\delta}^{\ast},BA^{-1}B^{\intercal}\bar{\delta}^{\ast}\rangle=-\left\langle\bar{\delta}^{\ast},\bar{g}_{1}(u_{1}^{\ast})\cdot\frac{\partial{\bar{\xi}_{1}}}{\partial
u_{2}}(u_{1}^{\ast},u_{2}^{\ast})\bar{\eta}^{\ast}\right\rangle.$

Therefore (27) holds.

##  Appendix C Proof of Proposition 2.1

We verify the intrinsic definition of DSE for $(x^{\ast},y^{\ast})$. For the
first-order condition in (4), we compute based on [7, Proposition 3.61],

$\partial_{x}f(x,y)=A^{\intercal}y,\quad\mbox{grad}_{y}f(x,y)=(I-yy^{\intercal})(Ax-b).$
(34)

From the identity involving the pseudo-inverse
$A^{\intercal}AA^{+}=A^{\intercal}$, we deduce that $A^{\intercal}y^{\ast}=0$.
As $b\not\in\mbox{Range}(A)$, we have $Ax^{\ast}-b\neq 0$ and therefore
$y^{\ast}$ is parallel to the (non-zero) vector $Ax^{\ast}-b$. From above, the
first-order condition holds, i.e.
$\partial_{x}f(x^{\ast},y^{\ast})=0,\mbox{grad}_{y}f(x^{\ast},y^{\ast})=0$.

From (34), we first verify the second-order condition (5). From [7, Corollary
5.16], we compute for $\eta\in T_{y}\mathcal{M}_{2}$,

$\displaystyle\mbox{Hess}_{y}f(x,y)[\eta]$
$\displaystyle=(I-yy^{\intercal})(\langle\partial_{y}\mbox{grad}_{y}f(x,y),\eta\rangle)$
$\displaystyle=(I-yy^{\intercal})(-(y\eta^{\intercal}+\eta
y^{\intercal})(Ax-b))$ $\displaystyle=\langle y,Ax-b\rangle(-\eta).$

We verify that $-\mbox{Hess}_{y}f(x^{\ast},y^{\ast})$ is d.p, because for non-
zero $\eta^{\ast}\in T_{y^{\ast}}\mathcal{M}_{2}$, we have
$\|\eta^{\ast}\|^{2}>0$. Moreover, $Ax^{\ast}-b\neq 0$, therefore

$\langle-\mbox{Hess}_{y}f(x^{\ast},y^{\ast})[\eta^{\ast}],\eta^{\ast}\rangle=\langle
y^{\ast},Ax^{\ast}-b\rangle\|\eta^{\ast}\|^{2}=\|Ax^{\ast}-b\|\|\eta^{\ast}\|^{2}>0.$

We now check the second-order condition (6). It is clear that
$\partial^{2}_{xx}f(x,y)=0$. We next show
$\delta^{\ast}\mapsto\mbox{grad}^{2}_{xy}f(x^{\ast},y^{\ast})[\delta^{\ast}]$
is an injection from $T_{x^{\ast}}\mathcal{M}_{1}$ to
$T_{y^{\ast}}\mathcal{M}_{2}$. We compute from (34),

$\mbox{grad}^{2}_{xy}f(x^{\ast},y^{\ast})[\delta^{\ast}]=\sum_{i\leq
d_{1}}\partial_{x_{i}}\mbox{grad}_{y}f(x^{\ast},y^{\ast})\delta^{\ast}_{i}=(I-y^{\ast}y^{\ast\intercal})A\delta^{\ast}.$

If $\delta^{\ast}\neq 0$, then $A\delta^{\ast}\neq 0$ because
$\mbox{Ker}(A)=\\{0\\}$. But $A\delta^{\ast}$ is not parallel to $y^{\ast}$,
since $y^{\ast}$ is along the direction $Ax^{\ast}-b$ which is not in the
range of $A$. This proves that
$\mbox{grad}^{2}_{xy}f(x^{\ast},y^{\ast})[\delta^{\ast}]\neq 0$ if
$\delta^{\ast}\neq 0$.

The above injection property implies that for $\delta^{\ast}\neq 0$,
$\eta^{\ast}=\mbox{grad}^{2}_{xy}f(x^{\ast},y^{\ast})[\delta^{\ast}]\neq 0$.
As $-\mbox{Hess}_{y}f(x^{\ast},y^{\ast})$ is d.p, we use the symmetry of the
Riemannian cross-gradients (see Proposition I.1) to obtain: If
$\delta^{\ast}\neq 0$, then

$\langle\delta^{\ast},[\mbox{Hess}_{x}f-\mbox{grad}^{2}_{yx}f\cdot(\mbox{Hess}_{y}f)^{-1}\cdot\mbox{grad}^{2}_{xy}f](x^{\ast},y^{\ast})[\delta^{\ast}]\rangle=\langle\eta^{\ast},-\mbox{Hess}_{y}f(x^{\ast},y^{\ast})^{-1}[\eta^{\ast}]\rangle>0.$

Therefore (6) holds.

##  Appendix D Proof of Proposition 2.2

As in the proof of Proposition 2.1 in Appendix C, we have

$\partial_{x}f(x,y)=A^{\intercal}y-\kappa
A^{\intercal}Ax,\quad\mbox{grad}_{y}f(x,y)=(I-yy^{\intercal})(Ax-b).$ (35)

To show the existence of DSE, we shall construct a solution of the form
$(x,y^{\ast}_{x})$, where $\partial_{x}f(x,y^{\ast}_{x})=0$. We consider
$y^{\ast}_{x}=(Ax-b)/\|Ax-b\|$ so that $\mbox{grad}_{y}f(x,y^{\ast}_{x})=0$.
We verify that $\mbox{Ker}(A)=\\{0\\}$ implies that $A^{\intercal}A$ is p.d.
It follows that $A^{+}=(A^{\intercal}A)^{-1}A^{\intercal}$, and that the
condition $\partial_{x}f(x,y^{\ast}_{x})=0$ is equivalent to

$\frac{A^{\intercal}(Ax-b)}{\|Ax-b\|}=\kappa
A^{\intercal}Ax\quad\Leftrightarrow\quad(1-\kappa\|Ax-b\|)x=A^{+}b.$

For $(x^{\ast},y^{\ast})=(x^{\ast},y^{\ast}_{x^{\ast}})$ to be a DSE, we
assume $x^{\ast}=cA^{+}b$ and we want to find a $c\in\mathbb{R}$ such that

$F(\kappa,c)=c(1-\kappa\|cAA^{+}b-b\|)-1=0.$

From Proposition 2.1, we have $F(0,1)=0$. We next apply the implicit function
theorem (c.f. Theorem A.1) to show the existence of $c$ close to $1$ if
$\kappa$ is close to $0$. It suffices to verify that $\frac{\partial
F}{\partial c}(0,1)=1$ based on

$\frac{\partial F}{\partial
c}(\kappa,c)=\left(1-\kappa\|cAA^{+}b-b\|\right)+c\left(-\kappa\frac{\langle
AA^{+}b,cAA^{+}b-b\rangle}{\|cAA^{+}b-b\|}\right).$

This suggests that there is an implicit function $h$ defined on an open
neighbor $W_{1}$ of $\kappa=0$ ($0\in W_{1}$) such that
$F(\kappa,h(\kappa))=0$. Moreover, $h(0)=1$. We next check that if
$c=h(\kappa)$ for some range $0<\kappa<\kappa_{0}$, then $(x^{\ast},y^{\ast})$
defined in the statement is a DSE of $f$. Since (35) hold at
$(x^{\ast},y^{\ast})$, it suffices to verify the second-order condition of
DSE.

We verify the second-order condition (5) by using the same proof of
Proposition 2.1. To check the second-order condition (6), we first compute
from (35)

$\eta^{\ast}=\mbox{grad}^{2}_{xy}f(x^{\ast},y^{\ast})[\delta^{\ast}]=\sum_{i\leq
d_{1}}\partial_{x_{i}}\mbox{grad}_{y}f(x^{\ast},y^{\ast})\delta^{\ast}_{i}=(I-y^{\ast}y^{\ast\intercal})A\delta^{\ast}.$

By following the proof of Proposition 2.1, we then compute

$\displaystyle S(\kappa)$
$\displaystyle=\langle\delta^{\ast},[\mbox{Hess}_{x}f-\mbox{grad}^{2}_{yx}f\cdot(\mbox{Hess}_{y}f)^{-1}\cdot\mbox{grad}^{2}_{xy}f](x^{\ast},y^{\ast})[\delta^{\ast}]\rangle$
$\displaystyle=-\kappa\|A\delta^{\ast}\|^{2}+\langle\eta^{\ast},-\mbox{Hess}_{y}f(x^{\ast},y^{\ast})^{-1}[\eta^{\ast}]\rangle$
$\displaystyle=-\kappa\|A\delta^{\ast}\|^{2}+\frac{1}{\|Ax^{\ast}-b\|}\|\eta^{\ast}\|^{2}.$

It is clear that $S(\kappa)$ is a continuous function of $\kappa$ because
$x^{\ast}=h(\kappa)A^{+}b$ and $y^{\ast}=(Ax^{\ast}-b)/\|Ax^{\ast}-b\|$ are
continuous with respect to $\kappa$. Moreover, $S(0)>0$ from the proof of
Proposition 2.1 (due to the injection property of
$\delta^{\ast}\mapsto\mbox{grad}^{2}_{xy}f(x^{\ast},y^{\ast})[\delta^{\ast}]$).
Therefore there exists $\kappa_{0}>0$ so that $S(\kappa)>0$ for
$0<\kappa<\kappa_{0}$ (i.e. (6) holds).

##  Appendix E Proof of Proposition 2.3

We verify the intrinsic definition of DNE for $(x^{\ast},y^{\ast})$. For the
first-order condition (7), we compute

$\partial_{x}f(x,y)=A^{\intercal}(Ax+y-b),\quad\mbox{grad}_{y}f(x,y)=(I-yy^{\intercal})(Ax+y-b).$
(36)

From the identity involving the pseudo-inverse
$A^{\intercal}AA^{+}=A^{\intercal}$, we deduce that
$A^{\intercal}(Ax^{\ast}+y^{\ast}-b)=0$. As $b\not\in\mbox{Range}(A)$, we have
$Ax^{\ast}-b\neq 0$ and therefore $y^{\ast}$ is parallel to the vector
$Ax^{\ast}-b$. From above, the first-order condition in (4) holds.

From (36), we next verify the second-order condition (8). It is clear that
$\partial^{2}_{xx}f(x,y)=A^{\intercal}A$ is p.d. From [7, Corollary 5.16], we
compute for $\eta\in T_{y}\mathcal{M}_{2}$,

$\displaystyle\mbox{Hess}_{y}f(x,y)[\eta]$
$\displaystyle=(I-yy^{\intercal})(\langle\partial_{y}\mbox{grad}_{y}f(x,y),\eta\rangle)$
$\displaystyle=(I-yy^{\intercal})(\eta-(y\eta^{\intercal}+\eta
y^{\intercal})(Ax+y-b))$ $\displaystyle=\eta(1-\langle y,Ax+y-b\rangle)$
$\displaystyle=\langle y,Ax-b\rangle(-\eta).$

We verify that $-\mbox{Hess}_{y}f(x^{\ast},y^{\ast})$ is d.p, because for non-
zero $\eta^{\ast}\in T_{y^{\ast}}\mathcal{M}_{2}$, we have
$\|\eta^{\ast}\|^{2}>0$. Moreover, $Ax^{\ast}-b\neq 0$, therefore

$\langle-\mbox{Hess}_{y}f(x^{\ast},y^{\ast})[\eta^{\ast}],\eta^{\ast}\rangle=\langle
y^{\ast},Ax^{\ast}-b\rangle\|\eta^{\ast}\|^{2}=\|Ax^{\ast}-b\|\|\eta^{\ast}\|^{2}>0.$

Therefore (8) holds.

##  Appendix F Proof of Theorem 3.1

It remains to prove that the spectral radius of the Jacobian matrix
$\bar{\mathbf{T}}^{\prime}(u_{1}^{\ast},u_{2}^{\ast})$ is strictly smaller
than one, we shall first compute this matrix
$\bar{\mathbf{T}}^{\prime}(u_{1}^{\ast},u_{2}^{\ast})$, and then relate it to
$\mathbf{T}^{\prime}(x^{\ast},y^{\ast})$.

First of all, we specify the induced dynamics $\bar{\mathbf{T}}$ of
$\mathbf{T}$ by using the local coordinate representation of the retractions
$\mathcal{R}_{1}$ and $\mathcal{R}_{2}$ near the stationary point
$(x^{\ast},y^{\ast})$. Let $(O_{1}\times O_{2},\varphi_{1}\times\varphi_{2})$
be the local chart. Since a retraction [1, Section 4.1] is a smooth function
from the tangent bundle of a manifold to the manifold itself, we can identify
an open subset $B_{1}$ of the tangent bundle of $\mathcal{M}_{1}$ (similarly
on $B_{2}$ for $\mathcal{M}_{2}$) such that $x\in O_{1}$ and
$\mathcal{R}_{1,x}(\xi_{1})\in O_{1}$ if $(x,\xi_{1})\in B_{1}$ This set
$B_{1}$ can then be mapped to an open set
$\bar{B}_{1}\subset\mathbb{R}^{d_{1}}\times\mathbb{R}^{d_{1}}$, using the
local chart of the tangent bundle ([1, Section 3.5.3]). On this subset
$\bar{B}_{1}$, we can define the local representation of $\mathcal{R}_{1}$ as
follows

$\bar{\mathcal{R}}_{1}:\bar{B}_{1}\to\bar{O}_{1},\quad\bar{\mathcal{R}}_{1}(u_{1},\bar{\xi}_{1})=\varphi_{1}\circ\mathcal{R}_{1,\varphi^{-1}_{1}(u_{1})}(D\varphi^{-1}_{1}(u_{1})[\bar{\xi}_{1}]),$

Similarly for $\mathcal{R}_{2}$, we can define
$\bar{\mathcal{R}}_{2}:\bar{B}_{2}\to\bar{O}_{2}$ on an open set
$\bar{B}_{2}\subset\mathbb{R}^{d_{2}}\times\mathbb{R}^{d_{2}}$. From the above
definition of $\bar{\mathcal{R}}_{1}$ and $\bar{\mathcal{R}}_{2}$, as well as
the local coordinate of $\xi_{1}(x,y)$ and $\xi_{2}(x,y)$ (as in (17),(18)),
we obtain the induced dynamics

$\bar{\mathbf{T}}(u_{1},u_{2})=\left(\begin{array}[]{l}\bar{\mathcal{R}}_{1}(u_{1},\bar{\xi}_{1}(u_{1},u_{2}))\\\
\bar{\mathcal{R}}_{2}(u_{2},\bar{\xi}_{2}(u_{1},u_{2}))\end{array}\right),$
(37)

which is well defined on the non-empty open set
$\\{(u_{1},u_{2})\in\bar{O}_{1}\times\bar{O}_{2}|(u_{1},\bar{\xi}_{1}(u_{1},u_{2}))\in\bar{B}_{1},(u_{2},\bar{\xi}_{2}(u_{1},u_{2}))\in\bar{B}_{2}\\}$.

From (37), we can compute the Jacobian matrix
$\bar{\mathbf{T}}^{\prime}(u_{1}^{\ast},u_{2}^{\ast})$ using the chain rule.
As $\bar{\xi}_{1}(u_{1}^{\ast},u_{2}^{\ast})=0$ and
$\bar{\xi}_{2}(u_{1}^{\ast},u_{2}^{\ast})=0$, we have

$\bar{\mathbf{T}}^{\prime}(u_{1}^{\ast},u_{2}^{\ast})=\left(\begin{array}[]{ll}I_{d_{1}}+\frac{\partial\bar{\xi}_{1}}{\partial
u_{1}}(u_{1}^{\ast},u_{2}^{\ast})&\frac{\partial\bar{\xi}_{1}}{\partial
u_{2}}(u_{1}^{\ast},u_{2}^{\ast})\\\ \frac{\partial\bar{\xi}_{2}}{\partial
u_{1}}(u_{1}^{\ast},u_{2}^{\ast})&I_{d_{2}}+\frac{\partial\bar{\xi}_{2}}{\partial
u_{2}}(u_{1}^{\ast},u_{2}^{\ast})\end{array}\right),$ (38)

since the retraction by definition satisfies
$\frac{\partial\bar{\mathcal{R}}_{1}}{\partial
u_{1}}(u_{1}^{\ast},0)=I_{d_{1}},\frac{\partial\bar{\mathcal{R}}_{1}}{\partial\bar{\xi}_{1}}(u_{1}^{\ast},0)=I_{d_{1}}$
(similarly for $\bar{\mathcal{R}}_{2}$).

We next verify that eigenvalues of the matrix in (38) are the same as
$\mathbf{T}^{\prime}(x^{\ast},y^{\ast})$. Let
$\delta^{\ast}=\sum_{i}\bar{\delta}_{i}^{\ast}E_{1,i}(x^{\ast})\in
T_{x^{\ast}}\mathcal{M}_{1}$,
$\eta^{\ast}=\sum_{j}\bar{\eta}_{j}^{\ast}E_{2,j}(y^{\ast})\in
T_{y^{\ast}}\mathcal{M}_{2}$, then following the same argument as in (24) and
(29), we have

$\displaystyle\nabla_{x}\xi_{1}(x^{\ast},y^{\ast})[\delta^{\ast}]$
$\displaystyle=\sum_{i=1}^{d_{1}}\left(\frac{\partial\bar{\xi}_{1,i}}{\partial
u_{1}}(u_{1}^{\ast},u_{2}^{\ast})\bar{\delta}^{\ast}\right)E_{1,i}(x^{\ast}),$
$\displaystyle D_{y}\xi_{1}(x^{\ast},y^{\ast})[\eta^{\ast}]$
$\displaystyle=\sum_{i=1}^{d_{1}}\left(\frac{\partial\bar{\xi}_{1,i}}{\partial
u_{2}}(u_{1}^{\ast},u_{2}^{\ast})\bar{\eta}^{\ast}\right)E_{1,i}(x^{\ast}),$
$\displaystyle D_{x}\xi_{2}(x^{\ast},y^{\ast})[\delta^{\ast}]$
$\displaystyle=\sum_{j=1}^{d_{2}}\left(\frac{\partial\bar{\xi}_{2,j}}{\partial
u_{1}}(u_{1}^{\ast},u_{2}^{\ast})\bar{\delta}^{\ast}\right)E_{2,j}(y^{\ast}),$
$\displaystyle\nabla_{y}\xi_{2}(x^{\ast},y^{\ast})[\eta^{\ast}]$
$\displaystyle=\sum_{j=1}^{d_{2}}\left(\frac{\partial\bar{\xi}_{2,j}}{\partial
u_{2}}(u_{1}^{\ast},u_{2}^{\ast})\bar{\eta}^{\ast}\right)E_{2,j}(y^{\ast}).$

From the definition of eigenvalue and eigenvector pairs, these equations
indicate that their eigenvalues are indeed the same.

##  Appendix G Proof of Theorem 3.2

Let
$\bar{\mathbf{T}}^{\prime}(u_{1}^{\ast},u_{2}^{\ast})=I_{d_{1}+d_{2}}+\gamma
M$. We show that the real part of each eigenvalue of $M$ is strictly smaller
than zero. From this, we can set
$\gamma^{\bullet}=-2\max_{k}\frac{\mbox{Re}(\lambda_{k}(M))}{|\lambda_{k}(M)|^{2}}>0$
and check that if $0<\gamma<\gamma^{\bullet}$, the spectral radius of
$\bar{\mathbf{T}}^{\prime}(u_{1}^{\ast},u_{2}^{\ast})$ is strictly smaller
than one.

From (38) and (31)-(33), the $M$ matrix is given by

$M=\left(\begin{array}[]{ll}-\bar{g}_{1}(u_{1}^{\ast})^{-1}\cdot
C&-\bar{g}_{1}(u_{1}^{\ast})^{-1}\cdot B\\\
\tau\bar{g}_{2}(u_{2}^{\ast})^{-1}\cdot
B^{\intercal}&-\tau\bar{g}_{2}(u_{2}^{\ast})^{-1}\cdot A\end{array}\right),$

where the matrices $A$,$B$,$C$ are defined in (26). To analyze the eigenvalue
of $M$, we could not apply [24, Lemma 5.2] to the matrix $M$ since the local
metric $\bar{g}_{1}(u_{1}^{\ast})$ and $\bar{g}_{2}(u_{2}^{\ast})$ are not
identity matrices. We consider a matrix $\mathbf{M}$ which is similar to $M$
so that they have the same eigenvalues,

$\mathbf{M}=\left(\begin{array}[]{ll}-\bar{\mathbf{C}}&-\bar{\mathbf{B}}\\\
\tau\bar{\mathbf{B}}^{\intercal}&-\tau\bar{\mathbf{A}}\end{array}\right)$ (39)

where

  * •

$\bar{\mathbf{C}}=\bar{g}_{1}(u_{1}^{\ast})^{-1/2}\cdot
C\cdot\bar{g}_{1}(u_{1}^{\ast})^{-1/2}$

  * •

$\bar{\mathbf{B}}=\bar{g}_{1}(u_{1}^{\ast})^{-1/2}\cdot
B\cdot\bar{g}_{2}(u_{2}^{\ast})^{-1/2}$

  * •

$\bar{\mathbf{A}}=\bar{g}_{2}(u_{2}^{\ast})^{-1/2}\cdot
A\cdot\bar{g}_{2}(u_{2}^{\ast})^{-1/2}$

We first check that the following conditions hold:

  * •

$\bar{\mathbf{A}}$ is p.d because $A$ is p.d.

  * •

$\bar{\mathbf{C}}+\bar{\mathbf{B}}\bar{\mathbf{A}}^{-1}\bar{\mathbf{B}}^{\intercal}$
is p.d because $C+BA^{-1}B^{\intercal}$ is p.d.

  * •

$\tau>\|\bar{\mathbf{C}}\|/\lambda_{min}(\bar{\mathbf{A}})$.

Indeed, we can relate these conditions to the following intrinsic quantities
on the manifold $\mathcal{M}_{1}\times\mathcal{M}_{2}$, by following the proof
of Theorem 2.1 (see Appendix B):

  * •

Let $\mathbf{A}=-\mbox{Hess}_{y}f(x^{\ast},y^{\ast})$, then $\mathbf{A}$ and
$\bar{\mathbf{A}}$ have the same eigenvalues.

  * •

Let $\mathbf{B}=\mbox{grad}^{2}_{yx}f(x^{\ast},y^{\ast})$,
$\mathbf{B}^{\intercal}=\mbox{grad}^{2}_{xy}f(x^{\ast},y^{\ast})$, then
$\mathbf{C}+\mathbf{B}\mathbf{A}^{-1}\mathbf{B}^{\intercal}$ and
$\bar{\mathbf{C}}+\bar{\mathbf{B}}\bar{\mathbf{A}}^{-1}\bar{\mathbf{B}}^{\intercal}$
have the same eigenvalues.

  * •

Let $\mathbf{C}=\mbox{Hess}_{x}f(x^{\ast},y^{\ast})$, then $\mathbf{C}$ and
$\bar{\mathbf{C}}$ have the same eigenvalues.

  * •

The condition
$\tau>\max_{k}|\lambda_{k}(\mathbf{C})|/\lambda_{min}(\mathbf{A})$ implies
that $\tau>\|\bar{\mathbf{C}}\|/\lambda_{min}(\bar{A})$.

###  G.1 Spectral analysis of $\mathbf{M}$

To analyze the eigenvalue of $\mathbf{M}$, we use the next proposition which
is adapted from [24, Lemma 5.2] with a refined range of $\tau$. For two real
symmetric matrices $A$ and $B$, we write $A\geq B$ if $A-B$ is semi-p.d. We
write $A>B$ if $A-B$ is p.d. For a symmetric matrix $A$, $\lambda_{max}(A)$
denotes the maximal eigenvalue of $A$, and $\lambda_{min}(A)$ denotes the
minimal eigenvalue of $A$.

######  Proposition G.1.

Let $A\in\mathbb{R}^{m\times m}$, $B\in\mathbb{R}^{n\times m}$,
$C\in\mathbb{R}^{n\times n}$ such that $A$ and $C+BA^{-1}B^{\intercal}$ are
positive definite. Assume that $\tau>\|C\|/\lambda_{min}(A)$, then any
eigenvalue $\lambda=\lambda_{0}+i\lambda_{1}$ (with
$\lambda_{0}\in\mathbb{R},\lambda_{1}\in\mathbb{R}$) of the following matrix,

$M=\left(\begin{array}[]{ll}-C&-B\\\ \tau B^{\intercal}&-\tau
A\end{array}\right)$

satisfies $\lambda_{0}<0$.

###### Proof.

We show that $\lambda_{0}\geq 0$ will lead to a contradiction, by following
the proof of [24, Lemma 5.2]. Assume that $\lambda$ is an eigenvalue of $M$,
i.e. $\mbox{det}(\lambda I-M)=0$. To compute $\mbox{det}(\lambda I-M)$, notice
that $\lambda I+\tau A$ is invertible since $\lambda_{0}I+\tau A$ is positive
definite ($A$ is positive definite and $\lambda_{0}\geq 0$). By the Schur
complement, we have

$\displaystyle\mbox{det}(\lambda I-M)$
$\displaystyle=\mbox{det}\left(\begin{array}[]{cc}C+\lambda I&B\\\ -\tau
B^{\intercal}&\tau A+\lambda I\end{array}\right)$ (42)
$\displaystyle=\mbox{det}(\lambda I+\tau A)\mbox{det}(H(\lambda))$ (43)

where $H(\lambda)=\lambda I+C+B(\lambda/\tau+A)^{-1}B^{\intercal}$. As
$\lambda I+\tau A$ is invertible, $\mbox{det}(\lambda I-M)=0$ implies that
$\mbox{det}(H(\lambda))=0$.

Let the spectral decomposition of $A$ be $U\Lambda_{A}U^{\intercal}$, with
orthogonal $U\in\mathbb{R}^{m\times m}$ and diagonal
$\Lambda_{A}=\mbox{diag}(\lambda_{1}(A),\cdots,\lambda_{m}(A))\in\mathbb{R}^{m\times
m}$. Then

$H(\lambda)=\lambda I+C+\tilde{B}D\tilde{B}^{\intercal},$

with $\tilde{B}=BU$ and $D=\mbox{diag}(d_{1},\cdots,d_{m})$ where

$d_{k}=\frac{1}{\lambda/\tau+\lambda_{k}(A)}=\frac{\lambda_{0}/\tau+\lambda_{k}(A)-i\lambda_{1}/\tau}{(\lambda_{0}/\tau+\lambda_{k}(A))^{2}+(\lambda_{1}/\tau)^{2}},\quad
1\leq k\leq m.$

It follows that if $\lambda_{0}>\|C\|$, then the real-part of $H(\lambda)$ is

$\mbox{Re}(H(\lambda))=\lambda_{0}I+C+\tilde{B}\mbox{Re}(D)\tilde{B}^{\intercal}\quad\mbox{p.d}.$
(44)

This is contradictory to the fact that $\mbox{det}(H(\lambda))=0$ [24,
Corollary 10.2]. On the other hand, if $0\leq\lambda_{0}\leq\|C\|$ and
$\lambda_{1}\neq 0$, then consider for $\beta\in\mathbb{R}$,

$\displaystyle\mbox{Re}(H(\lambda))+\frac{\tau\beta}{\lambda_{1}}\mbox{Im}(H(\lambda))$
$\displaystyle=\lambda_{0}I+C+\tilde{B}\mbox{Re}(D)\tilde{B}^{\intercal}+\frac{\tau\beta}{\lambda_{1}}(\lambda_{1}I+\tilde{B}\mbox{Im}(D)\tilde{B}^{\intercal})$
$\displaystyle=(\lambda_{0}+\tau\beta)I+C+\tilde{B}F\tilde{B}^{\intercal},$
(45)

where $F=\mbox{diag}(f_{1},\cdots,f_{m})$ with

$\displaystyle
f_{k}=\frac{\lambda_{0}/\tau+\lambda_{k}(A)-\beta}{(\lambda_{0}/\tau+\lambda_{k}(A))^{2}+(\lambda_{1}/\tau)^{2}},\quad
1\leq k\leq m.$ (46)

Take $\beta=\lambda_{min}(A)$, then $f_{k}\geq 0$ for each $k\leq m$. The
condition $\tau\lambda_{min}(A)>\|C\|$ implies $(\lambda_{0}+\tau\beta)I+C$ is
p.d. and together with (46), we have (45) is p.d., so it is contradictory to
the fact that $\mbox{det}(H(\lambda))=0$ [24, Lemma 10.1]. Lastly, if
$0\leq\lambda_{0}\leq\|C\|$ and $\lambda_{1}=0$, then from (44),

$H(\lambda)=\lambda_{0}I+C+\tilde{B}D\tilde{B}^{\intercal}$ (47)

with
$d_{k}=\frac{1}{\lambda_{0}/\tau+\lambda_{k}(A)}=\frac{1}{\lambda_{k}(A)}-\frac{\lambda_{0}/\tau}{(\lambda_{0}/\tau+\lambda_{k}(A))\lambda_{k}(A)}\geq\frac{1}{\lambda_{k}(A)}-\frac{\lambda_{0}/\tau}{\lambda_{min}(A)\lambda_{k}(A)}$.
It follows

$\displaystyle H(\lambda)$
$\displaystyle\geq\lambda_{0}I+C+\bigg{(}1-\frac{\lambda_{0}/\tau}{\lambda_{min}(A)}\bigg{)}BA^{-1}B^{\intercal}$
$\displaystyle=\bigg{(}1-\frac{\lambda_{0}/\tau}{\lambda_{min}(A)}\bigg{)}(C+BA^{-1}B^{\intercal})+\frac{\lambda_{0}/\tau}{\lambda_{min}(A)}C+\lambda_{0}I$

As $\tau>\|C\|/\lambda_{min}(A)$, $0\leq\lambda_{0}\leq\|C\|$ we have that
$0\leq\frac{\lambda_{0}/\tau}{\lambda_{min}(A)}<1$ and that
$I+\frac{1}{\tau\lambda_{min}(A)}C$ is p.d, i.e. we find that again
$H(\lambda)$ is p.d which is contradictory. In conclusion, $\lambda_{0}<0$. ∎

##  Appendix H Proof of Theorem 3.3

As in [10, Lemma 2.7], we apply the Ky Fan inequality [39, Lemma 4.20] to
analyze the symmetric part of the matrix $\mathbf{M}$ in (39): if
$\lambda_{k}(\mathbf{M})=\lambda_{0}+i\lambda_{1}$ is the $k$-th eigenvalue of
$\mathbf{M}$, then its real part
$\lambda_{0}\leq\lambda_{max}((\mathbf{M}+\mathbf{M}^{\intercal})/2)$. We show
that $\lambda_{max}((\mathbf{M}+\mathbf{M}^{\intercal})/2)<0$. This allows one
to set
$\gamma^{\bullet}=-2\max_{k}\frac{\mbox{Re}(\lambda_{k}(\mathbf{M}))}{|\lambda_{k}(\mathbf{M})|^{2}}$
as in the proof of Theorem 3.2 to conclude.

Based on the proof of [21, Proposition 26], we consider
($\bar{\mathbf{A}},\bar{\mathbf{B}},\bar{\mathbf{C}}$ are defined in (39))

$\displaystyle\mathbf{M}$
$\displaystyle=\left(\begin{array}[]{ll}0&\sqrt{\frac{1}{\tau}}I\\\
I&0\end{array}\right)\left(\begin{array}[]{ll}-\bar{\mathbf{C}}&-\bar{\mathbf{B}}\\\
\tau\bar{\mathbf{B}}^{\intercal}&-\tau\bar{\mathbf{A}}\end{array}\right)\left(\begin{array}[]{ll}0&\sqrt{\frac{1}{\tau}}I\\\
I&0\end{array}\right)^{-1}$
$\displaystyle=\left(\begin{array}[]{ll}\sqrt{\tau}\bar{\mathbf{B}}^{\intercal}&-\sqrt{\tau}\bar{\mathbf{A}}\\\
-\bar{\mathbf{C}}&-\bar{\mathbf{B}}\end{array}\right)\left(\begin{array}[]{ll}0&I\\\
\sqrt{\tau}I&0\end{array}\right)=\left(\begin{array}[]{ll}-\tau\bar{\mathbf{A}}&\sqrt{\tau}\bar{\mathbf{B}}^{\intercal}\\\
-\sqrt{\tau}\bar{\mathbf{B}}&-\bar{\mathbf{C}}\end{array}\right).$

From the Ky Fan inequality, the real part of each eigenvalue $\lambda$ of
$\mathbf{M}$ is upper bounded by the maximal eigenvalue of
$(\mathbf{M}+\mathbf{M}^{\intercal})/2$. We verify that indeed
$\lambda_{max}((\mathbf{M}+\mathbf{M}^{\intercal})/2)<0$ because
$\bar{\mathbf{A}}$ and $\bar{\mathbf{C}}$ are p.d., and $\tau>0$. Therefore,
the real part of $\lambda$ is strictly smaller than $0$.

##  Appendix I Derivation of deterministic $\tau$-SGA algorithm

The $\tau$-SGA algorithm modifies the vector field
$\xi(x,y)=(-\delta(x,y),\tau\eta(x,y))$ of $\tau$-GDA using the anti-symmetric
part of the Jacobian matrix of $\xi(x,y)$. As
$\delta(x,y)=\mbox{grad}_{x}f(x,y)$ and $\eta(x,y)=\mbox{grad}_{y}f(x,y)$, the
Jacobian matrix is the following linear transform on
$T_{x}\mathcal{M}_{1}\times T_{y}\mathcal{M}_{2}$ (which is a natural
extension from the Euclidean case),

$J(x,y)=\left(\begin{array}[]{ll}-\tilde{\mathbf{C}}(x,y)&-\tilde{\mathbf{B}}(x,y)\\\
\tau\tilde{\mathbf{B}}^{\intercal}(x,y)&-\tau\tilde{\mathbf{A}}(x,y)\end{array}\right)=\left(\begin{array}[]{ll}-\mbox{Hess}_{x}(x,y)&-\mbox{grad}^{2}_{yx}f(x,y)\\\
\tau\mbox{grad}^{2}_{xy}f(x,y)&\tau\mbox{Hess}_{y}(x,y)\end{array}\right).$
(48)

Here we introduce the notation
$\tilde{\mathbf{A}},\tilde{\mathbf{B}},\tilde{\mathbf{C}}$ in (48) to simplify
the equation. The $\tau$-SGA update rule is obtained from 555 Note that in the
original SGA rule [23, Proposition 5] the transpose of
$\frac{J(x,y)-J^{\intercal}(x,y)}{2}$ is considered since its definition of
$J$ has a sign difference compared to ours in (48).

$\displaystyle\xi(x,y)+\mu\left(\frac{J(x,y)-J^{\intercal}(x,y)}{2}\right)\xi(x,y)$
$\displaystyle=\xi(x,y)+\mu\frac{\tau+1}{2}\left(\begin{array}[]{ll}0&-\tilde{\mathbf{B}}(x,y)\\\
\tilde{\mathbf{B}}^{\intercal}(x,y)&0\end{array}\right)\xi(x,y)$
$\displaystyle=\left[\left(\begin{array}[]{l}-\delta\\\
\tau\eta\end{array}\right)+\mu\frac{\tau+1}{2}\left(\begin{array}[]{l}-\tau\tilde{\mathbf{B}}[\eta]\\\
-\tilde{\mathbf{B}}^{\intercal}[\delta]\end{array}\right)\right](x,y).$

We next show that the correction term, which is proportional to
$(\tau\tilde{\mathbf{B}}[\eta],\tilde{\mathbf{B}}^{\intercal}[\delta])$, is
orthogonal to the $\tau$-GDA direction $(-\delta,\tau\eta)$ at each $(x,y)$,
under the Riemannian metric on the tangent space $T_{x}\mathcal{M}_{1}\times
T_{y}\mathcal{M}_{2}$.

######  Proposition I.1.

For any $(x,y)\in\mathcal{M}_{1}\times\mathcal{M}_{2}$, we have

$\langle\tau\tilde{\mathbf{B}}(x,y)[\eta(x,y)],\delta(x,y)\rangle_{x}+\langle\tilde{\mathbf{B}}^{\intercal}(x,y)[\delta(x,y)],-\tau\eta(x,y)\rangle_{y}=0.$

###### Proof.

We follow the proof in Appendix B, by using a local coordinate chart
$(O_{1}\times O_{2},\varphi_{1}\times\varphi_{2})$ around the point $(x,y)$
rather than $(x^{\ast},y^{\ast})$. As in (17) and (18), this coordinate chart
maps $\delta(x,y)$ and $\eta(x,y)$ to their local coordinates
$\bar{\delta}(u_{1},u_{2})$ and $\bar{\eta}(u_{1},u_{2})$. It also induces a
canonical basis $\\{E_{1,i}(x)\\}_{i\leq d_{1}}$ on $T_{x}\mathcal{M}_{1}$ and
$\\{E_{2,j}(y)\\}_{j\leq d_{2}}$ on $T_{y}\mathcal{M}_{2}$ . By the definition
of cross-gradients,

$\displaystyle\tilde{\mathbf{B}}(x,y)[\eta(x,y)]$
$\displaystyle=D_{y}\mbox{grad}_{x}f(x,y)[\eta(x,y)]$
$\displaystyle=D_{y}\delta(x,y)[\eta(x,y)]$
$\displaystyle=\sum_{i}\left(\frac{\partial\bar{\delta}_{i}}{\partial
u_{2}}(u_{1},u_{2})\bar{\eta}(u_{1},u_{2})\right)E_{1,i}(x).$

Using the Riemannian metric $\bar{g}_{1}$ represented in the local coordinate
as in (19), it turns out that

$\langle\tilde{\mathbf{B}}(x,y)[\eta(x,y)],\delta(x,y)\rangle_{x}=\bar{\delta}(u_{1},u_{2})^{\intercal}\bar{g}_{1}(u_{1})\left(\frac{\partial\bar{\delta}}{\partial
u_{2}}(u_{1},u_{2})\bar{\eta}(u_{1},u_{2})\right).$ (49)

Similarly we have

$\langle\tilde{\mathbf{B}}^{\intercal}(x,y)[\delta(x,y)],\eta(x,y)\rangle_{y}=\bar{\eta}(u_{1},u_{2})^{\intercal}\bar{g}_{2}(u_{2})\left(\frac{\partial\bar{\eta}}{\partial
u_{1}}(u_{1},u_{2})\bar{\delta}(u_{1},u_{2})\right).$ (50)

We conclude that (49) equals to (50) because as in (32),(33)

$\bar{g}_{1}(u_{1})\frac{\partial\bar{\delta}}{\partial
u_{2}}(u_{1},u_{2})=\left(\bar{g}_{2}(u_{2})\frac{\partial\bar{\eta}}{\partial
u_{1}}(u_{1},u_{2})\right)^{\intercal}.$

∎

##  Appendix J Computational efficiency of $\tau$-SGA

We first discuss the computation of $\tilde{\mathbf{B}}[\eta]$ and
$\tilde{\mathbf{B}}^{\intercal}[\delta]$ when $\mathcal{M}_{1}$ (resp.
$\mathcal{M}_{2}$) is an embedded sub-manifold of
$\mathbb{R}^{d_{1}^{\prime}}$ (resp. $\mathbb{R}^{d_{2}^{\prime}}$). We then
propose a linear-time computational procedure using auto-differentiation when
$\mathcal{M}_{1}$ is Euclidean, which is applicable to Wasserstein GAN. When
$\mathcal{M}_{2}$ is also Euclidean, this procedure is equivalent to the one
proposed in [4].

To compute $\tilde{\mathbf{B}}[\eta]$ at $(x,y)$, we first fix
$\eta=\eta(x,y)\in T_{y}\mathcal{M}_{2}$. From the property of embedded sub-
manifolds, we use $\leftrightarrow$ to identify a tangent vector $\delta\in
T_{x}\mathcal{M}_{1}$ with a vector
$\delta^{\prime}=(\delta^{\prime}_{i})_{i\leq
d_{1}^{\prime}}\in\mathbb{R}^{d_{1}^{\prime}}$ (resp. $\eta\in
T_{y}\mathcal{M}_{2}$ with $\eta^{\prime}=(\eta^{\prime}_{i})_{i\leq
d_{2}^{\prime}}\in\mathbb{R}^{d_{2}^{\prime}}$).

This allows one to compute the cross-gradients in the embedded space, by

$\tilde{\mathbf{B}}[\eta](x,y)=D_{y}\mbox{grad}_{x}f(x,y)[\eta]\quad\leftrightarrow\quad
D_{y}\delta^{\prime}(x,y)[\eta^{\prime}]=(\langle(\mbox{grad}_{y}\delta^{\prime}_{i}(x,y))^{\prime},\eta^{\prime}\rangle_{y})_{i\leq
d_{1}^{\prime}}\in\mathbb{R}^{d_{1}^{\prime}}.$

Note that
$\langle(\mbox{grad}_{y}\delta^{\prime}_{i}(x,y))^{\prime},\eta^{\prime}\rangle_{y}$
is computed on $\mathbb{R}^{d_{2}^{\prime}}$ with a metric induced from
$T_{y}\mathcal{M}_{2}$. Assume that its computational time is
$O(d_{2}^{\prime})$ for each $i$. Then it takes
$O(d_{2}^{\prime}d_{1}^{\prime})$ to compute $\tilde{\mathbf{B}}[\eta](x,y)$
in the embedded space. We can obtain a similar cost for
$\tilde{\mathbf{B}}^{\intercal}[\delta](x,y)$.

When $\mathcal{M}_{1}=\mathbb{R}^{d_{1}}$, the computational complexity of
$\tilde{\mathbf{B}}[\eta](x,y)$ can be significantly reduced: $\forall i\leq
d_{1}$,

$D_{y}\partial_{x_{i}}f(x,y)[\eta]=\partial_{x_{i}}\langle\mbox{grad}_{y}f(x,y),\eta\rangle_{y}\quad\leftrightarrow\quad\partial_{x_{i}}\langle\eta^{\prime}(x,y),\eta^{\prime}\rangle_{y}.$
(51)

Importantly, one does not need to recompute
$\langle\eta^{\prime}(x,y),\eta^{\prime}\rangle_{y}$ for each $i$. Therefore,
the whole cost of $\tilde{\mathbf{B}}[\eta](x,y)$ is
$O(d_{1}+d_{2}^{\prime})$. Note that in (51), the term $\eta^{\prime}$ is
“detached”.

For $\tilde{\mathbf{B}}^{\intercal}[\delta](x,y)$, we “detach”
$\delta=\delta(x,y)\in T_{x}\mathcal{M}_{1}$ and compute

$\tilde{\mathbf{B}}^{\intercal}[\xi_{1}](x,y)=\sum_{i\leq
d_{1}}\partial_{x_{i}}\eta(x,y)\delta_{i}\quad\leftrightarrow\quad\partial_{\eta^{\prime\prime}}\left(\sum_{i\leq
d_{1}}\partial_{x_{i}}\langle\eta^{\prime}(x,y),\eta^{\prime\prime}\rangle\delta_{i}\right)_{|_{\eta^{\prime\prime}=0}}.$

This implies that we can first compute
$\partial_{x_{i}}\langle\eta^{\prime}(x,y),\eta^{\prime\prime}\rangle$ for
each $i$ as in (51). We then compute its sum with $\delta_{i}$ which takes
$O(d_{1})$. Finally an extra auto-differentiation is taken with respect to
$\eta^{\prime\prime}$ which costs $O(d_{2}^{\prime})$. The whole cost of
$\tilde{\mathbf{B}}^{\intercal}[\xi_{1}](x,y)$ is therefore
$O(d_{1}+d_{2}^{\prime})$. Note that in this case, we use the Euclidean metric
on $\mathbb{R}^{d_{2}^{\prime}}$ to evaluate
$\langle\eta^{\prime}(x,y),\eta^{\prime\prime}\rangle$ rather than the induced
Riemannian metric in (51).

In the Example 2 in Figure 2 and the MNIST (digit 0) in Table 2, we find that
the computational time per iteration $t$ of $\tau$-SGA is roughly 3-4 times
that of $\tau$-GDA.

###  J.1 Extension to stochastic $\tau$-SGA

We construct stochastic $\tau$-SGA through an unbiased estimation of terms in
the update rule of deterministic $\tau$-SGA. To achieve this, we compute
$\partial_{x_{i}}\langle\eta^{\prime}(x,y),\eta^{\prime}\rangle_{y}$ in (51)
using two mini-batches independently sampled from a training set, one to
estimate $\eta^{\prime}(x,y)$, the other to estimate $\eta^{\prime}$.
Similarly, we use the same two mini-batches to estimate $\eta^{\prime}(x,y)$
and $\delta$ in $\tilde{\mathbf{B}}^{\intercal}[\xi_{1}](x,y)$.

##  Appendix K Proof of Theorem 3.4

Since $f$ is twice continuously differentiable, we can apply Theorem 3.1 to
analyze the induced dynamics $\bar{\mathbf{T}}$ ($\mathbf{T}$ defined by
(13),(14)). Following (37), we obtain

$\bar{\mathbf{T}}(u_{1},u_{2})=\left(\begin{array}[]{l}\bar{\mathcal{R}}_{1}\left(u_{1},-\gamma\left[\bar{\delta}+\mu\frac{(\tau+1)\tau}{2}\frac{\partial\bar{\delta}}{\partial
u_{2}}\bar{\eta}\right](u_{1},u_{2})\right)\\\
\bar{\mathcal{R}}_{2}\left(u_{2},\gamma\tau\bar{\eta}(u_{1},u_{2})\right)\end{array}\right),$
(52)

Firstly, we remark that if we want to analyze the $\tau$-SGA method beyond the
asymptotic regime, one needs to consider this induced dynamics
$\bar{\mathcal{R}}_{2}\left(u_{2},\gamma\left[\tau\bar{\eta}-\mu\frac{\tau+1}{2}\frac{\partial\bar{\eta}}{\partial
u_{1}}\bar{\delta}\right](u_{1},u_{2})\right)$ for the variable $u_{2}$. The
analysis of the spectral radius of the Jacobian matrix
$\bar{\mathbf{T}}^{\prime}(u_{1}^{\ast},u_{2}^{\ast})$ is harder as one could
not easily adapt the proof of Proposition G.1 to this case. This remains an
interesting open question.

From (52), we compute the Jacobian matrix
$\bar{\mathbf{T}}^{\prime}(u_{1}^{\ast},u_{2}^{\ast})=I+\gamma M$, where

$M=\left(\begin{array}[]{ll}-\frac{\partial\bar{\delta}}{\partial
u_{1}}(u_{1}^{\ast},u_{2}^{\ast})&-\frac{\partial\bar{\delta}}{\partial
u_{2}}(u_{1}^{\ast},u_{2}^{\ast})\\\ \tau\frac{\partial\bar{\eta}}{\partial
u_{1}}(u_{1}^{\ast},u_{2}^{\ast})&\tau\frac{\partial\bar{\eta}}{\partial
u_{2}}(u_{1}^{\ast},u_{2}^{\ast})\end{array}\right)-\mu\frac{(\tau+1)\tau}{2}\left(\begin{array}[]{ll}[\frac{\partial\bar{\delta}}{\partial
u_{2}}\frac{\partial\bar{\eta}}{\partial
u_{1}}](u_{1}^{\ast},u_{2}^{\ast})&[\frac{\partial\bar{\delta}}{\partial
u_{2}}\frac{\partial\bar{\eta}}{\partial u_{2}}](u_{1}^{\ast},u_{2}^{\ast})\\\
0&0\end{array}\right)$

This can be reduced to the analysis of the eigenvalues of a similar matrix
$\mathbf{M}$ as in (39),

$\mathbf{M}=\left(\begin{array}[]{ll}-\bar{\mathbf{C}}&-\bar{\mathbf{B}}\\\
\tau\bar{\mathbf{B}}^{\intercal}&-\tau\bar{\mathbf{A}}\end{array}\right)+\mu\frac{(\tau+1)\tau}{2}\left(\begin{array}[]{ll}-\bar{\mathbf{B}}\bar{\mathbf{B}}^{\intercal}&\bar{\mathbf{B}}\bar{\mathbf{A}}\\\
0&0\end{array}\right)$

As in the proof of Theorem 3.2, we verify that
$\tau>\min(\|\bar{\mathbf{C}}\|,\|\bar{\mathbf{C}}+\theta\bar{\mathbf{B}}\bar{\mathbf{B}}^{\intercal}\|)/\lambda_{min}(\bar{\mathbf{A}})$,
and $\mu=\theta\frac{2}{\tau(\tau+1)}$ with $0\leq\theta\leq
1/\lambda_{max}(\bar{\mathbf{A}})$. From the Proposition K.1 (see next), we
conclude that the real-part of each eigenvalue of $\mathbf{M}$ is strictly
negative. Therefore as in the proof of Theorem 3.3, we can find
$\gamma^{\bullet}>0$ to guarantee the local convergence for
$0<\gamma<\gamma^{\bullet}$.

###  K.1 Spectral analysis of $\mathbf{M}$

The following proposition is adapted from Proposition G.1 to analyze the
eigenvalues of $\mathbf{M}$.

######  Proposition K.1.

Let $A\in\mathbb{R}^{m\times m}$, $B\in\mathbb{R}^{n\times m}$,
$C\in\mathbb{R}^{n\times n}$ such that $A$ and $C+BA^{-1}B^{\intercal}$ are
positive definite. Assume that $\tau>\min(\|C\|,\|C+\theta
BB^{\intercal}\|)/\lambda_{min}(A)$, and $\mu=\theta\frac{2}{\tau(\tau+1)}$
with $0\leq\theta\leq 1/\lambda_{max}(A)$, then any eigenvalue
$\lambda=\lambda_{0}+i\lambda_{1}$ (where
$\lambda_{0}\in\mathbb{R},\lambda_{1}\in\mathbb{R}$) of the following matrix,

$M=\left(\begin{array}[]{ll}-C&-B\\\ \tau B^{\intercal}&-\tau
A\end{array}\right)+\mu\frac{(\tau+1)\tau}{2}\left(\begin{array}[]{ll}-BB^{\intercal}&BA^{\intercal}\\\
0&0\end{array}\right)$

satisfies $\lambda_{0}<0$.

###### Proof.

We show that $\lambda_{0}\geq 0$ will lead to a contradiction, by following
the proof of Proposition G.1. Assume that $\lambda$ is an eigenvalue of $M$,
i.e. $\mbox{det}(\lambda I-M)=0$. To compute $\mbox{det}(\lambda I-M)$, notice
that $\lambda I+\tau A$ is invertible since $\lambda_{0}I+\tau A$ is positive
definite ($A$ is positive definite and $\lambda_{0}\geq 0$). By the Schur
complement, we have

$\displaystyle\mbox{det}(\lambda I-M)$
$\displaystyle=\mbox{det}\left(\begin{array}[]{cc}C+\theta
BB^{\intercal}+\lambda I&B-\theta BA\\\ -\tau B^{\intercal}&\tau A+\lambda
I\end{array}\right)$ (55) $\displaystyle=\mbox{det}(\lambda I+\tau
A)\mbox{det}(H(\lambda))$ (56)

where $H(\lambda)=\lambda I+C+\theta BB^{\intercal}+(B-\theta
BA)(\lambda/\tau+A)^{-1}B^{\intercal}$. As $\lambda I+\tau A$ is invertible,
$\mbox{det}(\lambda I-M)=0$ implies that $\mbox{det}(H(\lambda))=0$.

Let the spectral decomposition of $A$ be $U\Lambda_{A}U^{\intercal}$, with
orthogonal $U\in\mathbb{R}^{m\times m}$ and diagonal
$\Lambda_{A}=\mbox{diag}(\lambda_{1}(A),\cdots,\lambda_{m}(A))\in\mathbb{R}^{m\times
m}$. Then

$H(\lambda)=\lambda I+C+\theta
BB^{\intercal}+\tilde{B}D\tilde{B}^{\intercal},$

with $\tilde{B}=BU$ and $D=\mbox{diag}(d_{1},\cdots,d_{m})$ where

$d_{k}=\frac{1-\theta\lambda_{k}(A)}{\lambda/\tau+\lambda_{k}(A)}=(1-\theta\lambda_{k}(A))\frac{\lambda_{0}/\tau+\lambda_{k}(A)-i\lambda_{1}/\tau}{(\lambda_{0}/\tau+\lambda_{k}(A))^{2}+(\lambda_{1}/\tau)^{2}},\quad
1\leq k\leq m.$

It follows that if $\lambda_{0}>\min(\|C\|,\|C+\theta BB^{\intercal}\|)$, then
the real-part of $H(\lambda)$ is

$\mbox{Re}(H(\lambda))=\lambda_{0}I+C+\theta
BB^{\intercal}+\tilde{B}\mbox{Re}(D)\tilde{B}^{\intercal}\geq\lambda_{0}I+C+\theta
BB^{\intercal}\quad\mbox{p.d}$ (57)

This is contradictory to the fact that $\mbox{det}(H(\lambda))=0$ [24,
Corollary 10.2]. On the other hand, if
$0\leq\lambda_{0}\leq\min(\|C\|,\|C+\theta BB^{\intercal}\|)$ and
$\lambda_{1}\neq 0$, then consider for $\beta\in\mathbb{R}$,

$\displaystyle\mbox{Re}(H(\lambda))+\frac{\tau\beta}{\lambda_{1}}\mbox{Im}(H(\lambda))$
$\displaystyle=\lambda_{0}I+C+\theta
BB^{\intercal}+\tilde{B}\mbox{Re}(D)\tilde{B}^{\intercal}+\frac{\tau\beta}{\lambda_{1}}(\lambda_{1}I+\tilde{B}\mbox{Im}(D)\tilde{B}^{\intercal})$
$\displaystyle=(\lambda_{0}+\tau\beta)I+C+\theta
BB^{\intercal}+\tilde{B}F\tilde{B}^{\intercal},$ (58)

where $F=\mbox{diag}(f_{1},\cdots,f_{m})$ with

$\displaystyle
f_{k}=(1-\theta\lambda_{k}(A))\frac{\lambda_{0}/\tau+\lambda_{k}(A)-\beta}{(\lambda_{0}/\tau+\lambda_{k}(A))^{2}+(\lambda_{1}/\tau)^{2}},\quad
1\leq k\leq m.$ (59)

Take $\beta=\lambda_{min}(A)$, then $f_{k}\geq 0$ for each $k\leq m$. The
condition $\tau\lambda_{min}(A)>\min(\|C\|,\|C+\theta BB^{\intercal}\|)$
implies $(\lambda_{0}+\tau\beta)I+C+\theta BB^{\intercal}$ is p.d. and
together with (59), we have (58) is p.d., so it is contradictory to the fact
that $\mbox{det}(H(\lambda))=0$ [24, Lemma 10.1]. Lastly, if
$0\leq\lambda_{0}\leq\min(\|C\|,\|C+\theta BB^{\intercal}\|)$ and
$\lambda_{1}=0$, then from (57),

$H(\lambda)=\lambda_{0}I+C+\theta
BB^{\intercal}+\tilde{B}D\tilde{B}^{\intercal}$

with
$d_{k}=\frac{1-\theta\lambda_{k}(A)}{\lambda_{0}/\tau+\lambda_{k}(A)}=\frac{1-\theta\lambda_{k}(A)}{\lambda_{k}(A)}-\frac{(1-\theta\lambda_{k}(A))\lambda_{0}/\tau}{(\lambda_{0}/\tau+\lambda_{k}(A))\lambda_{k}(A)}\geq(1-\theta\lambda_{k}(A))(\frac{1}{\lambda_{k}(A)}-\frac{\lambda_{0}/\tau}{\lambda_{min}(A)\lambda_{k}(A)})$.
As $\tau>\min(\|C\|,\|C+\theta BB^{\intercal}\|)/\lambda_{min}(A)$ and
$0\leq\lambda_{0}\leq\min(\|C\|,\|C+\theta BB^{\intercal}\|)$ we have
$0\leq\frac{\lambda_{0}/\tau}{\lambda_{min}(A)}<1$ and it follows

$\displaystyle H(\lambda)$ $\displaystyle\geq\lambda_{0}I+C+\theta
BB^{\intercal}+\left(1-\frac{\lambda_{0}/\tau}{\lambda_{min}(A)}\right)BA^{-1}B^{\intercal}-\theta(1-\frac{\lambda_{0}/\tau}{\lambda_{min}(A)})\tilde{B}\tilde{B}^{\intercal}$
$\displaystyle=\left(1-\lambda_{0}c_{0}\right)(C+BA^{-1}B^{\intercal})+\lambda_{0}(c_{0}C+c_{0}\theta
BB^{\intercal}+I)$
$\displaystyle\geq\left(1-\lambda_{0}c_{0}\right)(C+BA^{-1}B^{\intercal})$
(60)

where $c_{0}=\frac{1}{\tau\lambda_{min}(A)}$. The last inequity (60) is due to
$\lambda_{0}\geq 0$ and the condition $1>c_{0}\min(\|C\|,\|C+\theta
BB^{\intercal}\|)$. They imply that $\lambda_{0}(c_{0}C+c_{0}\theta
BB^{\intercal}+I)\geq 0$. As $1-\lambda_{0}c_{0}\in(0,1]$, (60) implies that
$H(\lambda)$ is p.d which is contradictory. In conclusion, $\lambda_{0}<0$. ∎

##  Appendix L Details of numerical experiments

###  L.1 Choice of smooth non-linearity

We aim to build GAN models whose value function $f$ is twice continuously
differentiable, based on smooth non-linearities studied in [6]. For the
discriminator of Gaussian distribution, $\rho$ is a smooth approximation of
the absolute value non-linearity of the form

$\rho(a)=\sqrt{a^{2}+\epsilon^{2}},\quad\epsilon=10^{-6}.$

This function is twice continuously differentiable since
$\rho^{\prime\prime}(a)=\frac{\epsilon^{2}}{2(a^{2}+\epsilon^{2})^{3/2}}$.
Similarly, for the discriminator of MNIST digits, $\rho$ is a smooth ReLu non-
linearity,

$\rho(a)=(a+\sqrt{a^{2}+\epsilon^{2}})/2.$

###  L.2 DCGAN generator for MNIST

We use the default DCGAN generator to model MNIST (digit 0), implemented in
WGAN-GP666https://github.com/caogang/wgan-gp. From the 60k training samples of
MNIST, we take $4932$ training samples, $991$ validation samples of
$\phi_{data}$ to build the dataset MNIST (digit 0). We also make a slight
modification of the generator so that the function $x\mapsto G_{x}(Z)$ is
twice continuously differentiable at any $x\in\mathcal{M}_{1}$ for a given
$Z$. For this, each ReLu non-linearity in $Z\mapsto G_{x}(Z)$ is replaced by
the smooth ReLu $\rho$ in Appendix L.1.

###  L.3 Scattering CNN discriminator for MNIST

We construct a smooth Lipschitz-continuous discriminator for MNIST digits,

$D_{y}(\phi)=\langle v_{y},\rho(w_{y}\star P(\phi)+b_{y})\rangle,$

where $\rho$ is the smooth ReLu defined in Appendix L.1. Given $\phi$, this
makes the function $y\mapsto D_{y}(\phi)$ twice continuously differentiable at
any $y\in\mathcal{M}_{2}$. However, the function $\phi\mapsto D_{y}(\phi)$ is
not everywhere twice continuously differentiable due to the modulus non-
linearity in the scattering transform. We therefore replace this modulus non-
linearity $z\mapsto|z|=|z_{re}+iz_{im}|$ by
$z\mapsto\sqrt{z_{re}^{2}+z_{im}^{2}+\epsilon^{2}}$ for each complex number
input $z=z_{re}+iz_{im}$. This makes the function $\phi\mapsto D_{y}(\phi)$
twice continuously differentiable. As a consequence, the function
$(x,y)\mapsto f(x,y)$ is also twice continuously differentiable, which is
induced from the smoothness of $(x,y)\mapsto D_{y}(G_{x}(Z))$ and $y\mapsto
D_{y}(\phi_{data})$.

##### Scattering transform

The input $\phi$ of dimension $d=784$ is represented as an image of size
$28\times 28$. It is pre-processed by the wavelet scattering transform
$P(\phi)$ to extract stable edge-like information using Morlet wavelet of
different orientations and scales. We use the second-order scattering
transform with four wavelet orientations (between $[0,\pi)$) and two wavelet
scales. It first computes the convolution of $\phi$ with each wavelet filter,
then a smooth modulus non-linearity is applied to each feature map. This
computation is repeated one more time on each obtained feature maps and then a
low-pass filter is applied to each of the channels. The obtained scattering
features $P(\phi)$ is an image of size $9\times 9$ with $25$ channels
($I=25$,$n=9$).

##### Orthogonal CNN layer

The orthogonal CNN layer is parameterized by the kernel $w_{y}$ and bias
$b_{y}$. The kernel $w_{y}$ has $5\times 5$ spatial size ($k=5$). With a
suitable padding and stride (two by two), we obtain an output image of size
$5\times 5$ ($N=5$) with $J=256$ channels. Therefore the embedding space
dimension of $v_{y}$ is $JN^{2}=6400$.

###  L.4 Stiefel manifold geometry

For the discriminators of the Gaussian distribution and MNIST (digit 0), part
of the parameters in $y$ belong to Stiefel manifolds. To choose a Riemannian
metric on a Stiefel manifold (which is non-Euclidean), we use the one in [28,
equation 20]. We also use the SVD projection in [28, Proposition 12] as the
retraction $\mathcal{R}_{2}$ on each Stiefel manifold.

###  L.5 Initialization for local convergence

##### Simultaneous $\tau$-GDA initialization for Gaussian distribution

Starting from a random initialization of $(x,y)$, we apply the $\tau$-GDA to
build a pre-trained model. It is pre-trained with batch size $1000$, learning
rate $\gamma=0.0002$ and $\tau=100$ for $T=50000$ iterations.

##### Alternating $\tau$-GDA initialization for MNIST (digit 0)

The alternating $\tau$-GDA [40] is often used in the training of WGAN-GP [16],
and it can be extended to Riemannian manifold using the retraction. Starting
from a random initialization of $(x,y)$, we apply the alternating $\tau$-GDA
to build a pre-trained model. Each iteration amounts to perform $\tau=5$
gradient updates of $y$ (with learning rate $0.1$) and one gradient update of
$x$ (with learning rate $0.1$). It is pre-trained with batch size $128$ for a
total number of $10000$ iterations.