# Saddle Point Optimization with Approximate Minimization Oracle

Youhei Akimoto  University of Tsukuba & RIKEN AIP1-1-1 TennodaiTsukubaIbaraki,
Japan305-8573<EMAIL_ADDRESS>
(The DOI of the original article is as follows)

###### Abstract.

A major approach to saddle point optimization $\min_{x}\max_{y}f(x,y)$ is a
gradient based approach as is popularized by generative adversarial networks
(GANs). In contrast, we analyze an alternative approach relying only on an
oracle that solves a minimization problem approximately. Our approach locates
approximate solutions $x^{\prime}$ and $y^{\prime}$ to
$\min_{x^{\prime}}f(x^{\prime},y)$ and $\max_{y^{\prime}}f(x,y^{\prime})$ at a
given point $(x,y)$ and updates $(x,y)$ toward these approximate solutions
$(x^{\prime},y^{\prime})$ with a learning rate $\eta$. On locally strong
convex–concave smooth functions, we derive conditions on $\eta$ to exhibit
linear convergence to a local saddle point, which reveals a possible
shortcoming of recently developed robust adversarial reinforcement learning
algorithms. We develop a heuristic approach to adapt $\eta$ derivative-free
and implement zero-order and first-order minimization algorithms. Numerical
experiments are conducted to show the tightness of the theoretical results as
well as the usefulness of the $\eta$ adaptation mechanism.

Minimax Optimization, Saddle Point Optimization, Robust Optimization,
Convergence

††copyright: none††journalyear: 2021††conference: 2021 Genetic and
Evolutionary Computation Conference; July 10–14, 2021; Lille, France††doi:
10.1145/3449639.3459266††ccs: Mathematics of computing Continuous
optimization††ccs: Theory of computation Convergence and learning in
games††ccs: Theory of computation Theory of randomized search heuristics

##  1\. Introduction

We consider the following min–max optimization problem

(1) $\min_{x\in\mathbb{R}^{m}}\max_{y\in\mathbb{R}^{n}}f(x,y)\enspace,$

where $f:\mathbb{R}^{m}\times\mathbb{R}^{n}\to\mathbb{R}$ is the objective
function. This arises in many fields, including constrained optimization
(Cherukuri et al., 2017), robust optimization (Conn and Vicente, 2012; Qiu et
al., 2018), robust reinforcement learning (RL) (Pinto et al., 2017; Shioya et
al., 2018), and generative adversarial networks (GANs) (Goodfellow et al.,
2014; Salimans et al., 2016).

Unless $f$ is convex–concave, i.e., convex in $x$ and concave in $y$, it is
often computationally intractable to locate the global optimal solution. A
more realistic target is to locate a local min–max saddle point
$(x^{*},y^{*})$, which is a local minimum of $f$ in $x$ and a local maximum of
$f$ in $y$. In this study, we focus on locating a local min–max saddle point
of (1). Unless otherwise specified, by “a saddle point,” we mean a min–max
saddle point in this paper.

First-order approaches, as popularized by the success of GANs, are often
employed for this purpose. A simultaneous gradient method

(2)
$(x_{t+1},y_{t+1})=(x_{t},y_{t})+\eta(-\nabla_{x}f(x_{t},y_{t}),\nabla_{y}f(x_{t},y_{t})),$

has been analyzed for its local and global convergence properties. On $C^{2}$
functions, the continuous time dynamics obtained for the limit $\eta\to 0$ are
analyzed and their local asymptotic stability around strict local saddle
points is shown (Nagarajan and Kolter, 2017; Cherukuri et al., 2017). The
analysis is extended to (2) with a finite $\eta>0$, and the condition on
$\eta$ for a strict local saddle point to be asymptotically stable is derived
in (Mescheder et al., 2017). Adolphs et al. (2019) has shown that there are
locally asymptotically stable points that are not saddle points and proposed a
modification to avoid these undesired stable points. Liang and Stokes (2019)
have shown that (2) globally converges toward the strict saddle point if the
objective function is globally and strongly convex–concave and that some
improved gradient-based approaches (Daskalakis et al., 2018; Yadav et al.,
2018; Mescheder et al., 2017) can converge toward a nonstrict saddle point on
a bilinear function. For constrained cases, the problems are often treated as
variational inequalities, and first-order approaches are employed (Gidel et
al., 2017; Nouiehed et al., 2019).

Zero-order approaches for (1) are categorized as coevolutionary approaches
(Al-Dujaili et al., 2019; Qiu et al., 2018; Jensen, 2004; Branke and
Rosenbusch, 2008), Bayesian optimization approaches (Picheny et al., 2019;
Bogunovic et al., 2018), trust-region approaches (Conn and Vicente, 2012), and
gradient approximation approaches (Liu et al., 2020). They are often designed
heuristically, and little attention has been paid to convergence guarantees
and convergence rate analysis of these approaches. It is known that
coevolutionary approaches suffer from nonconvergent behavior (Al-Dujaili et
al., 2019; Qiu et al., 2018). Recently, (Bogunovic et al., 2018) showed regret
bounds for a Bayesian optimization approach and (Liu et al., 2020) showed an
error bound for a gradient approximation approach, where the error is measured
by the square norm of the gradient. Both analyses show sublinear rates under
possibly stochastic (i.e., noisy) versions of (1).

In contrast to previous works, we analyze an alternative approach that relies
solely on approximate minimization oracles, which can be zero-order, first-
order, or higher-order approaches. Given $(x,y)$, this approach first locates
approximate solutions $x^{\prime}$ and $y^{\prime}$ to the minimization
problems
$\operatornamewithlimits{argmin}_{x^{\prime}\in\mathbb{R}^{m}}f(x^{\prime},y)$
and
$\operatornamewithlimits{argmin}_{y^{\prime}\in\mathbb{R}^{n}}-f(x,y^{\prime})$
by using approximate minimization oracles. Then, it updates $(x,y)$ toward
$(x^{\prime},y^{\prime})$ with a learning rate $\eta>0$, i.e.,
$(x,y)\leftarrow(1-\eta)(x,y)+\eta(x^{\prime},y^{\prime})$. One can choose
domain-specific methods to minimize $f(\cdot,y)$ and $-f(x,\cdot)$. This work
is motivated by a recent advance in robust RL (Pinto et al., 2017; Shioya et
al., 2018), where the objective is to locate a robust policy $x$ under
adversarial disturbance $y$, and these parameters are updated in the
abovementioned manner; the approximate minimization oracles are some standard
RL approaches. The applications of our analysis are of course not limited to
robust RL. In particular, it is suitable for problems involving numerical
simulations and requiring zero-order approaches not necessarily connected to
gradient descent.

The contributions of this paper are twofold. We analyze the local and global
convergence properties of an oracle-based saddle point optimization on locally
and globally strong convex–concave smooth functions. We derive a sufficient
condition on the learning rate $\eta$ to guarantee linear (i.e., geometric)
convergence and derive an upper bound on the convergence rate. In contrast to
the abovementioned analysis for zero-order approaches showing a sublinear
decrease (Bogunovic et al., 2018; Liu et al., 2020), our analysis is for
linear convergence; hence, the result is more related to the one obtained for
a simultaneous gradient update (Mescheder et al., 2017; Liang and Stokes,
2019). We show the condition on $\eta$ to be not only sufficient but also
necessary for convergence on a convex–concave quadratic function. This reveals
a possible shortcoming of approaches with $\eta=1$, which do not guarantee
convergence but are employed in e.g. robust adversarial RL (Pinto et al.,
2017; Shioya et al., 2018) or coevolutionary approaches. The tightness and
possible room for improvement are demonstrated in numerical experiments.

We propose a heuristic method to adapt the learning rate $\eta$. Our approach
is derivative free and black box, that is, with no gradient information and no
characteristic constant such as a smoothness parameter. Because $f$ is often
black-box and no sufficient information is available to select the right
$\eta$ in advance, especially when zero-order approaches are desired, $\eta$
adaptation is essential for practical approaches. We instantiate the whole
algorithm using a zero-order randomized hill-climbing approach, namely, the
(1+1)-ES with 1/5 success rule (Devroye, 1972; Schumer and Steiglitz, 1968;
Rechenberg, 1973) as well as sequential least squares programming (SLSQP)
(Kraft, 1988). We demonstrate that our learning rate adaptation method locates
a nearly optimal saddle point by iterating up to three times more $f$-calls
compared with the best fixed $\eta$.

##  2\. Formulation

#### Notation

Suppose that $f:\mathbb{R}^{m}\times\mathbb{R}^{n}\to\mathbb{R}$ is a twice
continuously differentiable function, i.e.,
$f\in\mathcal{C}^{2}(\mathbb{R}^{m}\times\mathbb{R}^{n},\mathbb{R})$. Let
$H_{x,x}(x,y)$, $H_{x,y}(x,y)$, $H_{y,x}(x,y)$, and $H_{y,y}(x,y)$ be the
blocks of the Hessian matrix of $f$, whose $(i,j)$-th component is
$\partial^{2}f/\partial x_{i}\partial x_{j}$, $\partial^{2}f/\partial
x_{i}\partial y_{j}$, $\partial^{2}f/\partial y_{i}\partial x_{j}$, and
$\partial^{2}f/\partial y_{i}\partial y_{j}$, respectively, evaluated at a
given point $(x,y)$. Let $J_{g}(z)$ be the Jacobian of a differentiable
$g=(g_{1},\dots,g_{k}):\mathbb{R}^{\ell}\to\mathbb{R}^{k}$, where the ($i$,
$j$)-th element is $\partial g_{i}/\partial z_{j}$ evaluated at
$z=(z_{1},\dots,z_{\ell})\in\mathbb{R}^{\ell}$. If $k=1$, we write
$J_{g}(z)=\nabla g(z)^{\mathrm{T}}$. For a positive definite symmetric matrix
$A$, let $\sqrt{A}$ denote the matrix square root.

#### Saddle Point

A point $(x^{*},y^{*})\in\mathbb{R}^{m}\times\mathbb{R}^{n}$ is a (strict)
local saddle point of a function
$f:\mathbb{R}^{m}\times\mathbb{R}^{n}\to\mathbb{R}$ if there exists a
neighborhood
$\mathcal{E}_{x}\times\mathcal{E}_{y}\subseteq\mathbb{R}^{m}\times\mathbb{R}^{n}$
including $(x^{*},y^{*})$ such that for any
$(x,y)\in\mathcal{E}_{x}\times\mathcal{E}_{y}$, the condition
$f(x,y^{*})\geqslant(>)f(x^{*},y^{*})\geqslant(>)f(x^{*},y)$ holds. If
$\mathcal{E}_{x}=\mathbb{R}^{m}$ and $\mathcal{E}_{y}=\mathbb{R}^{n}$,
$(x^{*},y^{*})$ is called the (strict) global saddle point. For
$f\in\mathcal{C}^{2}(\mathbb{R}^{m}\times\mathbb{R}^{n},\mathbb{R})$, a point
$(x^{*},y^{*})$ is a strict saddle point if it is a critical point
($\nabla_{x}f(x^{*},y^{*})=0$ and $\nabla_{y}f(x^{*},y^{*})=0$) and
$H_{x,x}(x^{*},y^{*})\succ 0$ and $H_{y,y}(x^{*},y^{*})\prec 0$ both hold.

#### Convex–concave Function

A function $f:\mathbb{R}^{m}\times\mathbb{R}^{n}\to\mathbb{R}$ is a (strict)
convex–concave function if $f$ is (strictly) convex in $x$ and (strictly)
concave in $y$. It is locally (and strictly) convex–concave in open sets
$\mathcal{E}_{x}\subseteq\mathbb{R}^{m}$ and
$\mathcal{E}_{y}\subseteq\mathbb{R}^{n}$ if the restriction of $f$ to
$\mathcal{E}_{x}\times\mathcal{E}_{y}$ is a (strict) convex–concave function.
Moreover, $f$ is called a (locally) strong convex–concave function if (the
restriction to $\mathcal{E}_{x}\times\mathcal{E}_{y}$ of) $f$ is strongly
convex in $x$ and strongly concave in $y$. For
$f\in\mathcal{C}^{2}(\mathbb{R}^{m}\times\mathbb{R}^{n},\mathbb{R})$, $f$ is
(locally) strong convex–concave if and only if there exists $\mu>0$ such that
$H_{x,x}(x,y)\succcurlyeq\mu$ and $H_{y,y}(x,y)\preccurlyeq-\mu$ for all
$(x,y)\in\mathbb{R}^{m}\times\mathbb{R}^{n}$
($(x,y)\in\mathcal{E}_{x}\times\mathcal{E}_{y}$).

#### Suboptimality Error

We define a quantity to measure the progress of the optimization towards the
saddle point. For a strict convex–concave function, there exists a unique
saddle point $(x^{*},y^{*})$. The suboptimality error (Gidel et al., 2017) is
then defined as

(3) $\displaystyle G(x,y)$
$\displaystyle:=\max_{y^{\prime}\in\mathbb{R}^{n}}f(x,y^{\prime})-\min_{x^{\prime}\in\mathbb{R}^{m}}f(x^{\prime},y)\enspace.$

It is easy to see that the suboptimality error is non-negative and is zero
only at the saddle point $(x^{*},y^{*})$. If $f$ is locally strictly
convex–concave in a neighborhood $\mathcal{E}_{x}\times\mathcal{E}_{y}$ of a
local saddle point $(x^{*},y^{*})$, the suboptimality error can be extended to
a nonconvex–concave scenario as

(4) $\displaystyle G(x,y)$
$\displaystyle:=\max_{y^{\prime}\in\mathcal{E}_{y}}f(x,y^{\prime})-\min_{x^{\prime}\in\mathcal{E}_{x}}f(x^{\prime},y)\enspace.$

In this study, we use a quadratic approximation of (3) and (4) to measure the
progress toward a saddle point $(x^{*},y^{*})$. If
$f\in\mathcal{C}^{2}(\mathbb{R}^{m}\times\mathbb{R}^{n},\mathbb{R})$ and
(locally) strong convex–concave around a (local) saddle point $(x^{*},y^{*})$,
the (local) suboptimality error can be approximated around $(x^{*},y^{*})$ by

(5) $\tilde{G}(x,y)=\lVert x-x^{*}\rVert_{G_{x,x}^{*}}^{2}/2+\lVert
y-y^{*}\rVert_{G_{y,y}^{*}}^{2}/2\enspace,$

where (dropping $(x,y)$ from $H_{x,x}$, $H_{y,y}$, $H_{x,y}$, and $H_{y,x}$)

(6) $\begin{split}G_{x,x}(x,y)=&H_{x,x}+H_{x,y}(-H_{y,y})^{-1}H_{y,x}\\\
G_{y,y}(x,y)=&(-H_{y,y})+H_{y,x}H_{x,x}^{-1}H_{x,y}\end{split}$

and

(7) $\begin{split}\lVert
x-x^{*}\rVert_{G_{x,x}^{*}}=&[(x-x^{*})^{\mathrm{T}}G_{x,x}(x^{*},y^{*})(x-x^{*})]^{1/2},\\\
\lVert
y-y^{*}\rVert_{G_{y,y}^{*}}=&[(y-y^{*})^{\mathrm{T}}G_{y,y}(x^{*},y^{*})(y-y^{*})]^{1/2}.\end{split}$

It is non-negative for all $(x,y)\in\mathbb{R}^{m}\times\mathbb{R}^{n}$ and
zero only at $(x^{*},y^{*})$. If $f$ is a convex–concave quadratic function,
we have $G(x,y)=\tilde{G}(x,y)$ for all
$(x,y)\in\mathbb{R}^{m}\times\mathbb{R}^{n}$. In the following we use
$\tilde{G}(x,y)$ to measure progress.

##  3\. Algorithm

#### Approximate Minimization Oracle

Consider a minimization of $g:\mathbb{R}^{\ell}\to\mathbb{R}$. Given a
positive definite symmetric matrix $A\in\mathbb{R}^{\ell\times\ell}$, consider
a local minimization of $g$ in the neighborhood
$U_{A}(z,r):=\\{z^{\prime}\in\mathbb{R}^{\ell}:\lVert
z^{\prime}-z\rVert_{A}\leqslant r\\}$ of an initial search point $z$ with
radius $r>0$ under $\lVert z\rVert_{A}=[z^{\mathrm{T}}Az]^{1/2}$,

(8) $\displaystyle z^{*}=\operatornamewithlimits{argmin}_{z^{\prime}\in
U_{A}(z,r)}g(z^{\prime})\enspace.$

A local minimization oracle $M_{\epsilon,r}^{A}(g,z)$ takes the function $g$
to be minimized, and the initial solution $z$ as inputs and outputs an
approximate solution $z^{\prime}$ to the above minimization problem (8)
satisfying the following condition:

(9) $\displaystyle\lVert
M_{\epsilon,r}^{A}(g,z)-z^{*}\rVert_{A}^{2}\leqslant\epsilon\cdot\lVert
z-z^{*}\rVert_{A}^{2}\enspace.$

In other words, it locates a point that decreases the squared distance from
the exact local minimum under $A$ by the factor $\epsilon\in(0,1)$ compared to
the initial point $z$.

Important examples are algorithms that exhibit linear convergence, where the
runtime (number of $f$-calls and/or $\nabla f$-calls) to shrink the distance
to the optimum by the factor $\epsilon$ is estimated by
$\Theta(\log(1/\epsilon))$. Therefore, by running $\Theta(\log(1/\epsilon))$
iterations of such algorithms, one can implement the above oracle. However, it
is not limited to linearly convergent algorithms. The oracle requirement (9)
can be satisfied with algorithms that converge sublinearly. In such cases, the
number of iterations required to satisfy (9) becomes greater as the candidate
solution approaches the optimum. Therefore, the stopping condition ($\tau$ in
our proposed algorithm, see Section 5) for the search algorithm inside the
oracle needs to be tuned more carefully.

#### Oracle-based Saddle Point Optimization

We consider an approach to the saddle point optimization (1) based solely on
approximate local minimization oracles satisfying (9). We first find
approximate local solutions to
$\operatornamewithlimits{argmin}_{x^{\prime}\in\mathbb{R}^{m}}f(x^{\prime},y_{t})$
and
$\operatornamewithlimits{argmin}_{y^{\prime}\in\mathbb{R}^{n}}-f(x_{t},y^{\prime})$
by approximate minimization oracles $M_{\epsilon_{x},r_{x}}^{A_{x}}$ and
$M_{\epsilon_{y},r_{y}}^{A_{y}}$, respectively.111We remark that if the right-
hand side of condition (9) is replaced with a constant, the whole minimax
algorithm cannot be guaranteed to converge to a local saddle point. Because
our objective is to derive the linear convergence rate, such a situation is
beyond the scope of this paper. Then, $x_{t}$ and $y_{t}$ are updated with the
learning rate $\eta>0$ as

(10) $\begin{split}x_{t+1}&=(1-\eta)x_{t}+\eta\cdot
M_{\epsilon_{x},r_{x}}^{A_{x}}(f(\cdot,y_{t}),x_{t})\enspace,\\\
y_{t+1}&=(1-\eta)y_{t}+\eta\cdot
M_{\epsilon_{y},r_{y}}^{A_{y}}(-f(x_{t},\cdot),y_{t})\enspace.\end{split}$

##  4\. Convergence Analysis

We investigate the global and local convergence properties of the proposed
approach (10) on
$f\in\mathcal{C}^{2}(\mathbb{R}^{m}\times\mathbb{R}^{n},\mathbb{R})$. We are
especially interested in knowing how small the learning rate $\eta$ needs to
be and how quickly it converges.

In the following, let $(x^{*},y^{*})$ be a strict saddle point of $f$ and $f$
be locally strong convex–concave around $(x^{*},y^{*})$. For notational
simplicity, let $H_{x,x}^{*}=H_{x,x}(x^{*},y^{*})$,
$H_{y,y}^{*}=H_{y,y}(x^{*},y^{*})$, $H_{x,y}^{*}=H_{x,y}(x^{*},y^{*})$,
$H_{y,x}^{*}=H_{y,x}(x^{*},y^{*})$, $G_{x,x}^{*}=G_{x,x}(x^{*},y^{*})$, and
$G_{y,y}^{*}=G_{y,y}(x^{*},y^{*})$.

Our analysis is based on the implicit function theorem (e.g., Theorem 5 of (de
Oliveira, 2013)), which shows the existence and uniqueness of the solutions
$\hat{x}(y)$ and $\hat{y}(x)$ to $\nabla_{x}f(\hat{x}(y),y)=0$ and
$\nabla_{y}f(x,\hat{y}(x))=0$. The proofs of the following results are
provided in Appendix A.

###  4.1. Global Linear Convergence

Suppose that
$f\in\mathcal{C}^{2}(\mathbb{R}^{m}\times\mathbb{R}^{n},\mathbb{R})$ is a
strong convex–concave function. Then, there exists a unique global saddle
point $(x^{*},y^{*})$. We assume that for

(11)
$\begin{split}\Delta_{x,y}(x,y)&=(H_{x,x}(x,y))^{-1}H_{x,y}(x,y)-(H_{x,x}^{*})^{-1}H_{x,y}^{*}\\\
\Delta_{x,y}(x,y)&=(H_{y,y}(x,y))^{-1}H_{y,x}(x,y)-(H_{y,y}^{*})^{-1}H_{y,x}^{*}\enspace,\end{split}$

there exists $\delta\in[0,1)$ such that

(12)
$\displaystyle\sigma(\sqrt{\smash[b]{G_{x,x}^{*}}}\Delta_{x,y}(x,y)\sqrt{\smash[b]{G_{y,y}^{*}}}^{-1})\leqslant\delta$
$\displaystyle\sigma(\sqrt{\smash[b]{G_{y,y}^{*}}}\Delta_{y,x}(x,y)\sqrt{\smash[b]{G_{x,x}^{*}}}^{-1})\leqslant\delta$

hold for all $(x,y)\in\mathbb{R}^{m}\times\mathbb{R}^{n}$, where
$\sigma(\cdot)$ denotes the greatest singular value of the argument matrix. We
let the class of such functions be denoted by $\mathcal{F}$. An instance in
$\mathcal{F}$ is a convex–concave quadratic function, where $\delta=0$.

The following theorem states the conditions on the approximation error
$\epsilon_{x}$ and $\epsilon_{y}$ of the approximate minimization oracle (9)
and the learning rate $\eta$ to guarantee the global linear convergence and
shows the upper bound of the convergence rate.

###### Theorem 4.1.

Suppose $f\in\mathcal{F}$ and $(x^{*},y^{*})$ is the global saddle point of
$f$. Let $\tilde{G}(x,y)$ be defined as (5). Let $\bar{\sigma}$ be the
greatest singular value of
$\sqrt{\smash[b]{G_{x,x}^{*}}}(H_{x,x}^{*})^{-1}H_{x,y}^{*}\sqrt{\smash[b]{G_{y,y}^{*}}}^{-1}$,
which is equal to
$\sqrt{\smash[b]{G_{x,x}^{*}}}^{-1}H_{x,y}^{*}(-H_{y,y}^{*})^{-1}\sqrt{\smash[b]{G_{y,y}^{*}}}$.
Suppose that there are approximate minimization oracles
$M_{\epsilon_{x},r_{x}}^{A_{x}}$ and $M_{\epsilon_{y},r_{y}}^{A_{y}}$
satisfying oracle condition (9) with $A_{x}=G_{x,x}^{*}$ and
$A_{y}=G_{y,y}^{*}$. Consider algorithm (10) with the neighborhood parameter
$r_{x}=r_{y}=\infty$. Let $\bar{\epsilon}=\max[\epsilon_{x},\epsilon_{y}]$. If

(13) $\displaystyle\bar{\epsilon}<1-\delta$ $\displaystyle\mathrm{and}$
$\displaystyle\eta<\bar{\eta}:=\frac{2(1-(\bar{\epsilon}+\delta))}{1+\bar{\sigma}^{2}-(\bar{\epsilon}+\delta)^{2}}\enspace,$

then for any $(x_{0},y_{0})\in\mathbb{R}^{m}\times\mathbb{R}^{n}$,

(14) $\tilde{G}(x_{t},y_{t})\leqslant\gamma^{2t}\tilde{G}(x_{0},y_{0})$

holds for
$\gamma=((1-\eta)^{2}+\eta^{2}\bar{\sigma}^{2})^{1/2}+\eta(\bar{\epsilon}+\delta)<1$.

###  4.2. Local Linear Convergence

If $f\in\mathcal{C}^{2}(\mathbb{R}^{m}\times\mathbb{R}^{n},\mathbb{R})$, and
it is locally strong convex–concave in a neighborhood around a local saddle
point $(x^{*},y^{*})$, we can derive a local convergence condition, as shown
in Theorem 4.2. This provides a tighter estimate of the asymptotic rate
$\bar{\gamma}$ of convergence than that in Theorem 4.1.

###### Theorem 4.2.

Let $f\in\mathcal{C}^{2}(\mathbb{R}^{m}\times\mathbb{R}^{n},\mathbb{R})$ and
locally strong convex–concave in a neighborhood around a local saddle point
$(x^{*},y^{*})$. Let $\tilde{G}(x,y)$ be defined as (5). Suppose that there
are approximate minimization oracles $M_{\epsilon_{x},r_{x}}^{A_{x}}$ and
$M_{\epsilon_{y},r_{y}}^{A_{y}}$ satisfying oracle condition (9) with
$A_{x}=G_{x,x}^{*}$ and $A_{y}=G_{y,y}^{*}$. Consider algorithm (10) with
neighborhood parameters $r_{x}$ and $r_{y}$. Let $\bar{\sigma}$ and
$\bar{\epsilon}$ be defined in Theorem 4.1. If

(15) $\displaystyle\bar{\epsilon}<1$ and
$\displaystyle\eta<\bar{\eta}:=\frac{2(1-\bar{\epsilon})}{1+\bar{\sigma}^{2}-\bar{\epsilon}^{2}}\enspace,$

then for any $\gamma\in(\bar{\gamma},1)$ with
$\bar{\gamma}=((1-\eta)^{2}+\eta^{2}\bar{\sigma}^{2})^{1/2}+\eta\bar{\epsilon}$,
there exist constants $\bar{r}_{x}>0$ and $\bar{r}_{y}>0$ such that for any
$r_{x}\in(0,\bar{r}_{x})$ and $r_{y}\in(0,\bar{r}_{y})$, there exists a
neighborhood $U\subseteq\mathbb{R}^{m}\times\mathbb{R}^{n}$ of $(x^{*},y^{*})$
satisfying

(16)
$\tilde{G}(x_{t},y_{t})\leqslant\gamma^{2t}\tilde{G}(x_{0},y_{0}),\quad\forall(x_{0},y_{0})\in
U.$

###  4.3. Discussion

Theorems 4.1 and 4.2 show that the approximate suboptimality error
$\tilde{G}(x_{t},x_{t})$ converges linearly in terms of the number of queries
to the approximate minimization oracle as long as the conditions are
satisfied. If the approximate minimization oracle requires $\mathcal{O}(1)$
$f$-calls (or $\nabla f$-calls) on average, this implies that
$\tilde{G}(x_{t},y_{t})$ converges linearly in terms of $\\#f$-calls (or
$\nabla f$-calls) as well. As mentioned in Section 3, a linearly converging
minimization approach requires $\mathcal{O}(\log(1/\epsilon))$ $f$-calls, and
hence the resulting approach (10) converges linearly in terms of queries. This
result is distinguished from the results of Bogunovic et al. (2018); Liu et
al. (2020), where a sublinear decrease has been analyzed in a stochastic
(noisy) setting, rather than a linear decrease in a deterministic one.

Theorem 4.2 provides the upper bound $\bar{\gamma}$ of the convergence rate
$\lim_{t\to\infty}\sqrt{\smash[b]{\tilde{G}(x_{t+1},y_{t+1})/\tilde{G}(x_{t},y_{t})}}$.
It is minimized at $\eta=\eta^{*}$, and its optimal value is
$\bar{\gamma}^{*}$, where

(17)
$\begin{split}\eta^{*}&:=\frac{1-\sqrt{\bar{\sigma}^{2}\bar{\epsilon}^{2}/(1-\bar{\epsilon}^{2}+\bar{\sigma}^{2})}}{1+\bar{\sigma}^{2}}\enspace,\\\
\bar{\gamma}^{*}&:=\frac{\bar{\sigma}(1+\bar{\sigma}^{2})+(\bar{\epsilon}+\bar{\epsilon}^{2}\bar{\sigma})\sqrt{1-\bar{\epsilon}^{2}+\bar{\sigma}^{2}}}{(1+\bar{\sigma}^{2})\sqrt{1-\bar{\epsilon}^{2}+\bar{\sigma}^{2}}}\enspace.\end{split}$

We have $\eta^{*}=1/(1+\bar{\sigma}^{2})$ and
$\bar{\gamma}^{*}=\bar{\sigma}/\sqrt{1+\bar{\sigma}^{2}}$ for
$\bar{\epsilon}=0$ (an exact minimization oracle) and $\eta^{*}\downarrow 0$
and $\bar{\gamma}^{*}\uparrow 1$ as $\bar{\epsilon}\uparrow 1$. The number of
iterations to halve $\tilde{G}$ is
$T\leqslant\log(2)/\log(1/\bar{\gamma}^{2})$.

Conditions (13) and (15) are shown to be sufficient in Theorems 4.1 and 4.2.
The question naturally arises as to whether it is necessary as well. For
example, consider $f(x,y)=x^{2}/2+\bar{\sigma}xy-y^{2}/2$. Then,
$G_{x,x}^{*}=G_{y,y}^{*}=1+\bar{\sigma}^{2}$, where $\bar{\sigma}$ is
equivalent to $\bar{\sigma}$ in Theorem 4.1. We have
$\hat{y}(x)=\bar{\sigma}x$ and $\hat{x}(y)=-\bar{\sigma}y$. For the exact
minimization oracle ($\bar{\epsilon}=0$), we have
$x_{t+2}=[(1-\eta)^{2}+\eta^{2}\bar{\sigma}^{2}]x_{t}$ and
$y_{t+2}=[(1-\eta)^{2}+\eta^{2}\bar{\sigma}^{2}]y_{t}$, and it is easy to see
that $(x_{t},y_{t})$ converges if and only if
$(1-\eta)^{2}+\eta^{2}\bar{\sigma}^{2}<1$, i.e.,
$\eta<2/(1+\bar{\sigma}^{2})$, the right-hand side of which is $\bar{\eta}$ in
(13) and (15). Therefore, conditions (13) and (15) are also necessary for such
a situation. Robust adversarial reinforcement learning (RARL) (Pinto et al.,
2017) and its extension (Shioya et al., 2018) fall into our framework with
$\eta=1$. Our theoretical investigation reveals a possible limitation in these
two and a need for the introduction of $\eta<1$.

A local convergence analysis of a simultaneous gradient update (2) in
(Mescheder et al., 2017) shows that $(x_{t},y_{t})$ converges toward
$(x^{*},y^{*})$ if and only if all the eigenvalues of the matrix
$\Big{[}\begin{smallmatrix}I-\eta H_{x,x}^{*}&-\eta H_{x,y}^{*}\\\ \eta
H_{y,x}^{*}&I+\eta H_{y,y}^{*}\end{smallmatrix}\Big{]}$ live in the unit disk
on a complex space. Again, considering the abovementioned case of
$f(x,y)=x^{2}/2+\bar{\sigma}xy-y^{2}/2$, the absolute values of the
eigenvalues in a complex space of the above matrix are all
$((1-\eta)^{2}+\eta^{2}\bar{\sigma}^{2})^{1/2}$. That is, the necessary and
sufficient condition $\eta<2/(1+\bar{\sigma}^{2})$ for the simultaneous
gradient update (2) is equivalent to the condition on $\eta$ for (10) with
$\bar{\epsilon}=0$.

For the global convergence in Theorem 4.1, we assume that there exists
$\delta\in[0,1)$ satisfying (12) for all
$(x,y)\in\mathbb{R}^{m}\times\mathbb{R}^{n}$. In contrast, Theorem 1 of (Liang
and Stokes, 2019) shows the global linear convergence of the simultaneous
gradient approach with a sufficiently small $\eta$ on a globally and strongly
convex–concave $f$. Our condition $\delta<1$ is stronger. We suspect that this
is not a potential difference between the simultaneous gradient approach and
the zero-order saddle point optimization (10) and that it is not necessary for
the global convergence itself. The looseness of the result will be
demonstrated in Section 6.

##  5\. Learning Rate Adaptation

Often, $\bar{\sigma}$ is unknown, especially when zero-order optimization is
required. Here we propose a heuristic approach to adapt the learning rate
$\eta$. We implemented zero- and first-order saddle point optimization
algorithms using a randomized hill-climbing algorithm, namely, the (1+1)-ES
with 1/5 success rule (Devroye, 1972; Schumer and Steiglitz, 1968; Rechenberg,
1973), and a sequential least squares programming (SLSQP) subroutine (Kraft,
1988). Our implementation is publicly
available222https://gist.github.com/youheiakimoto/15212dbf46dc546af20af38b0b48ff17.

###  5.1. Adaptation Mechanism

The suboptimality error $G(x_{t},y_{t})$ is approximated by using the oracle
outputs $x^{\prime}=M^{A_{x}}_{\epsilon_{x},r_{x}}(f(\cdot,y_{t}),x_{t})$ and
$y^{\prime}=M^{A_{y}}_{\epsilon_{y},r_{y}}(-f(x_{t},\cdot),y_{t})$ as
$G(x_{t},y_{t})\approx
f(x_{t},y^{\prime})-f(x^{\prime},y_{t})=:F(x_{t},y_{t})$. More precisely, we
have $(1-2\bar{\epsilon})G(x_{t},y_{t})\leqslant F(x_{t},y_{t})\leqslant
G(x_{t},y_{t})$. This implies that if there exists $\tilde{\gamma}<1$ such
that $F(x_{t},y_{t})\leqslant\tilde{\gamma}^{2t}F(x_{0},y_{0})$, we can
guarantee that
$G(x_{t},y_{t})\leqslant(1-2\bar{\epsilon})^{-1}\tilde{\gamma}^{2t}G(x_{0},y_{0})$.
We note also that $G\approx\tilde{G}$ in a neighborhood of $(x^{*},y^{*})$ for
$\mathcal{C}^{2}$ functions.

Algorithm 1 Learning Rate Adaptation

1:$x\in\mathbb{R}^{m}$, $y\in\mathbb{R}^{n}$, $a_{\eta}>0$, $b_{\eta}\geqslant
0$, $c_{\eta}>1$

2:$\eta\leftarrow 1$, $\tilde{\gamma}\leftarrow 0$

3:for $t=1,\cdots,T$ do

4: $x_{t}\leftarrow x$, $y_{t}\leftarrow y$

5: $\eta_{c}\leftarrow\\{\min(\eta\cdot c_{\eta},1),\eta,\eta/c_{\eta}\\}$
w.p. $1/3$ for each

6: $N_{\text{step}}\leftarrow\lfloor b_{\eta}+a_{\eta}/\eta_{c}\rfloor$

7: for $s=1,\cdots,N_{\text{step}}$ do

8: $\hat{x}\leftarrow M^{A_{x}}_{\epsilon_{x},r_{x}}(f(\cdot,y),x)$

9: $\hat{y}\leftarrow M^{A_{y}}_{\epsilon_{y},r_{y}}(-f(x,\cdot),y)$

10: $F_{s}\leftarrow f(x,\hat{y})-f(\hat{x},y)$

11: $x\leftarrow x+\eta_{c}(\hat{x}-x)$

12: $y\leftarrow y+\eta_{c}(\hat{y}-y)$

13: break if $s\geqslant b_{\eta}$ and $F_{s}>\cdots>F_{s-b_{\eta}+1}$

14: end for

15: $\tilde{\gamma}_{c},\
\sigma_{\tilde{\gamma}_{c}}\leftarrow\textsc{slope}(\log(F_{1}),\dots,\log(F_{s}))$

16: if $\tilde{\gamma}\geqslant 0$ and $\tilde{\gamma}_{c}\geqslant 0$ then

17: $\eta\leftarrow\eta/c_{\eta}^{3}$

18: else if $\tilde{\gamma}_{c}\leqslant\tilde{\gamma}$ or $\eta=\eta_{c}$
then

19: $\eta\leftarrow\eta_{c}$, $\tilde{\gamma}\leftarrow\tilde{\gamma}_{c}$

20: end if

21: $x\leftarrow x_{t}$, $y\leftarrow y_{t}$ if
$\tilde{\gamma}-2\sigma_{\tilde{\gamma}_{c}}>0$

22:end for

From the above observation, we propose to adapt the learning rate $\eta>0$ by
Algorithm 1. The object is to find the value of $\eta$ that minimizes the
convergence rate of $F(x_{t},y_{t})$. We perform a random search for the
optimal $\eta$. The learning rate $\eta_{c}$ to be tested is one of $\eta$,
$\eta\cdot c_{\eta}$, or $\eta/c_{\eta}$ chosen at random with a probability
of $1/3$ each. The logarithm of the convergence rate, $\tilde{\gamma}_{c}$, is
estimated by running the algorithm (10) with $\eta_{c}$ for $N_{\text{step}}$
iterations (Lines 7 to 15). Based on the estimated log convergence rate
$\tilde{\gamma}_{c}$, we search for the best $\eta$.

To estimate the convergence rate, we run $N_{\text{step}}\propto 1/\eta$
iterations. The rationale behind the choice of $N_{\text{step}}\propto 1/\eta$
is as follows. From Theorem 4.1, we know that the logarithm of the convergence
rate, $\log(\gamma)$ in (14), is $-\mathcal{O}(\eta)$. On the other hand,
because of the approximation error, i.e.,
$G(x_{t},y_{t})\leqslant(1-2\bar{\epsilon})^{-1}\tilde{\gamma}^{2t}G(x_{0},y_{0})$,
the logarithm of the convergence rate of $F$ estimated over $N_{\text{step}}$
iterations can deviate from that of $G$ by
$\pm\lvert\log(1-2\bar{\epsilon})\rvert/N_{\text{step}}$. To alleviate the
effect of the approximation error, we need
$\epsilon/N_{\text{step}}\in\mathcal{O}(\eta)$. In words, $N_{\text{step}}$
needs to be no less than proportional to $1/\eta$.

Algorithm 1 summarizes the overall framework of the proposed approach with the
adaptation mechanism.333In practice, it is advised to try out a fresh $y$
after every $M^{A_{y}}_{\epsilon_{y},r_{y}}(-f(x_{t},\cdot),y_{t})$ calls if
the minimization oracle is for a local minimization but it is desired to avoid
(rather) suboptimal $y$. After Line 9, one may add the following lines:
1:$\tilde{y}\sim p_{y}$ 2:if $f(x,\tilde{y})>f(x,\hat{y})$ then 3:
$\hat{y}\leftarrow\tilde{y}$ 4:end if A possible choice of $p_{y}$ is the
distribution from which the initial $y$ is drawn. If the search space of $y$
is bounded, the uniform distribution over the search domain of $y$ is a
candidate for $p_{y}$. One can also perform analogous steps for $x$. We didn’t
implement it in our experiments as the objective of the experiments are to
understand the proposed approach and the adaptation mechanism.  We introduce
three hyperparameters $a_{\eta}$, $b_{\eta}$, and $c_{\eta}$ to adapt one
parameter $\eta$. Arguably, tuning them is easier than tuning $\eta$ itself. A
greater $a_{\eta}$ and $b_{\eta}$ will lead to a more accurate estimation of
$\tilde{\gamma}_{c}$, while spending more $f$-calls and requiring a longer
adaptation time for $\eta$. The granularity of the randomized line search for
$\eta$ is controlled by $c_{\eta}$. A smaller $c_{\eta}$ will lead to a
smoother change in $\eta$, while requiring a longer adaptation time.

Algorithm 2 (1+1)-ES

1:$h:\mathbb{R}^{\ell}\to\mathbb{R}$, $z\in\mathbb{R}^{\ell}$,
$\sigma\in(0,\sigma_{\max}]$, $h_{z}=h(z)$, $\tau_{\text{es}}\in\mathbb{N}$

2:$c=e^{1/\sqrt{2\ell}}$

3:$n_{\text{succ}}=0$

4:while $n_{\text{succ}}<\tau_{\text{es}}\cdot\ell$ do

5: $z^{\prime}\leftarrow z+\sigma\mathcal{N}(0,I)$

6: $h_{z^{\prime}}\leftarrow h(z^{\prime})$

7: if $h_{z^{\prime}}\leqslant h_{z}$ then

8: $\sigma\leftarrow\min(\sigma\cdot c,\sigma_{\max})$, $z\leftarrow
z^{\prime}$, $h_{z}\leftarrow h_{z^{\prime}}$

9: $n_{\text{succ}}\leftarrow n_{\text{succ}}+1$

10: else

11: $\sigma\leftarrow\sigma\cdot c^{-1/4}$

12: end if

13:end while

14:return $z$, $h_{z}$, $\sigma$

(a) AES ($n=10$, $b=1$)

(b) AES (varying $n$, $b=1$)

(c) AES ($n=10$, varying $b$)

(d) ASLSQP ($n=10$, $b=1$)

(e) ASLSQP (varying $n$, $b=1$)

(f) ASLSQP ($n=10$, varying $b$)

Figure 1. The number of $f$-calls until $G(x_{t},y_{t})\leqslant 10^{-5}$ is
obtained on $f_{1}$. (a–c) AES. (d–f) ASLSQP. (a, d) Boxplots of the results
with and without $\eta$ adaptation. Missing data imply that no run has
succeeded in achieving $G(x_{t},y_{t})\leqslant 10^{-5}$. (b, e) Median (solid
line) and inter-quartile range (band) for dimensions $n=m$. (c, f) Median
(solid line) and inter-quartile range (band) for interaction term $b$.

###  5.2. Adversarial Evolution Strategy

We approximate $M_{\epsilon,r}^{A}$ with the so-called (1+1)-ES with a
$1/5$-success rule. The (1+1)-ES is a randomized hill-climbing algorithm with
an adaptive step-size mechanism. We adopt the simplified update of the step-
size proposed by (Kern et al., 2004). Algorithm 2 summarizes the (1+1)-ES on
$h:\mathbb{R}^{\ell}\to\mathbb{R}$. In our case, $h$ corresponds to
$f(\cdot,y_{t})$ for $x$ updates and $-f(x_{t},\cdot)$ for $y$ updates. It
generates a candidate solution $z^{\prime}$ from a normal distribution
centered at $z$ with standard deviation $\sigma$. If the candidate solution
has a function value smaller than or equal to that of the current solution, we
replace the current solution with the candidate and call each such iteration
successful. The step size is increased by multiplying $c$ if the iteration is
successful. Otherwise, it is decreased by multiplying $c^{-1/4}$. Such a step-
size adaptation mechanism is called a $1/5$-success rule because it maintains
the step size to keep the probability of success to $1/5$.

The runtime time analysis of the (1+1)-ES algorithm reveals that the expected
number $\mathbb{E}[T]$ of iterations $T$ to achieve $\lVert
z-z^{*}\rVert\leqslant\epsilon R_{0}$ from the initial condition $R_{0}=\lVert
z-z^{*}\rVert$ is $\Theta(\log(1/\epsilon))$ on strongly convex functions with
Lipschitz continuous gradients and their strictly increasing transformations
(Morinaga and Akimoto, 2019). Moreover, its scaling w.r.t. dimension $\ell$ is
proven to be $\Theta(\ell)$ on convex quadratic functions (Morinaga et al.,
2021). Therefore, we expect that Algorithm 2 approximates $M_{\epsilon,r}^{A}$
with $\epsilon\in\exp(-\Theta(T/\ell))$.

We propose a zero-order saddle point optimization algorithm using the
(1+1)-ES, called _Adversarial Evolution Strategy (AES)_. AES replaces
$M^{A_{x}}_{\epsilon_{x},r_{x}}(f(\cdot,y_{i}),x_{i})$ and
$M^{A_{y}}_{\epsilon_{y},r_{y}}(-f(x_{i},\cdot),y_{i})$ in Algorithm 1 with
Algorithm 2. In Algorithm 1, in addition to $x$ and $y$, AES maintains two
step-size parameters $\sigma_{x}$ for $x$ update and $\sigma_{y}$ for $y$
update. On Lines 8 and 9, we input $\sigma_{x}$ and $\sigma_{y}$ to Algorithm
2 to obtain the updated values. On Lines 4 and 21, we maintain $\sigma_{x}$
and $\sigma_{y}$ as well as $x$ and $y$.

Importantly, the number $T$ of iterations of Algorithm 2 is not fixed for each
oracle call. Instead, we count the number of successful iterations, where a
better or equally accurate objective function value is generated. The reason
for this design choice is as follows. As mentioned above, we expect
$\epsilon\in\exp(-\Theta(T/\ell))$. However, if, for example, $\sigma$ is too
large to produce a successful candidate solution444It is reported in (Morinaga
et al., 2021) that $\sigma$ needs to be proportional to the norm of the
gradient to produce a successful candidate solution if the objective function
is convex quadratic and that the (1+1)-ES with 1/5 success rule maintains
$\sigma$ to be proportional to the norm of the gradient., Algorithm 2 has
little chance to produce any improvement until $\sigma$ is well adapted.
Because the correct $\sigma$ approaches zero as the candidate solution
approaches the optimum, it requires more iterations to adapt $\sigma$ if it is
initialized to a fixed value for each oracle call. Then, the output of
Algorithm 2 with a fixed number of iterations will not satisfy the oracle
condition (9) for some $\epsilon$. To reduce the effect of an excessively
large initial $\sigma$, we count the number of successful iterations rather
than total iterations.

To further relax the effect of an unsuitable initial $\sigma$ value, AES
shares $\sigma$ for successive oracle calls. With this parameter being shared
between oracle calls, we expect the following. Compared to the case where the
initial $\sigma$ is too large, we expect this to reduce the number of
iterations (objective function calls) to satisfy the stopping condition
because the adaptation time for $\sigma$ is reduced. Compared to the case
where the initial $\sigma$ is too small, where the iteration can be successful
with high probability despite the improvement being rather small, we expect to
avoid excessively early stopping of Algorithm 2.555Our preliminary experiments
revealed that the $\sigma$ sharing effect is statistically significant;
however, it may not be practically important. The ratio between the numbers of
$f$-calls to reach the same target threshold with AES and AES without $\sigma$
sharing was less than 1.1 on $f_{1}$ with $n=m=10$ and $b=2$. See Section 6
for the detailed experimental setting. The reason is simply that $\sigma$
adaptation in the (1+1)-ES is rather quick. From (Akimoto et al., 2018), we
know that the adaptation to forget an inaccurate initial step size is
proportional to $\log(\sigma_{0}/\sigma^{*})$, where $\sigma_{0}$ is the
initial step size and $\sigma^{*}$ denotes the optimal step size, which is
proportional to $\lVert x_{0}-x^{*}\rVert$ on the spherical function with
optimum at $x^{*}$.

###  5.3. Adversarial SLSQP

To demonstrate the applicability of the proposed $\eta$ adaptation mechanism,
we also implement it with SLSQP as $M_{\epsilon,r}^{A}$. This is a sequential
quadratic programming approach. We allow access to the gradient of $f$ to the
SLSQP procedure and limit the maximum number of iterations to
$\tau_{\mathrm{slsqp}}$. In contrast to AES, we share nothing but the
solutions $x$ and $y$ between different oracle calls. We call this first-order
approach _Adversarial SLSQP (ASLSQP)_.

##  6\. Experiments

We conducted numerical experiments with AES and ASLSQP with and without $\eta$
adaptation. Our objective was to confirm the tightness of the bounds obtained
in Theorems 4.1 and 4.2 and the applicability of the proposed $\eta$
adaptation mechanism. AES and ASLSQP without $\eta$ adaptation repeat Lines 8
to 12 with a fixed learning rate. We set $a_{\eta}=1$, $b_{\eta}=5$,
$c_{\eta}=1.1$, and $\tau_{\text{es}}=\tau_{\text{slsqp}}=5$ for the following
experiments, unless otherwise specified.

###  6.1. Ex 1: Convex–concave Quadratic Case

In this experiment, we confirm that (A) $\bar{\eta}$ in (13) is a tight upper
bound; (B) $\bar{\eta}^{*}$ in (17) is a good estimate of the optimal $\eta$;
(C) $\bar{\eta}$ in (13) and $\eta^{*}$ and $\bar{\gamma}^{*}$ in (17)
accurately reflect the dependency of $\bar{\sigma}$; and (D) the proposed
$\eta$ adaptation demonstrates a reasonable performance compared to the best
fixed $\eta$.

We consider the following convex–concave quadratic function
$f:\mathbb{R}^{m}\times\mathbb{R}^{n}\to\mathbb{R}$:

$f_{1}(x,y)=\lVert x\rVert^{2}/2+bx^{\mathrm{T}}y-\lVert
y\rVert^{2}/2\enspace.$

For this function, $\bar{\sigma}=b$ and $\delta=0$ in Theorem 4.1. The upper
bound on $\eta$ for linear convergence is $\bar{\eta}=2/(1+b^{2})$. Moreover,
$\bar{\gamma}^{*}=\lvert b\rvert/(1+b^{2})^{1/2}$ and $\eta^{*}=1/(1+b^{2})$
if $\bar{\epsilon}=0$666The reason that we approximated $\bar{\gamma}^{*}$ and
$\eta^{*}$ with $\epsilon=0$ in the experiments is because the rigorous
formula (17) is less understandable. In the experiments, we observed that the
oracle condition (9) is satisfied with $\epsilon<5\times 10^{-4}$ for 50% of
iterative checks and $\epsilon<5\times 10^{-3}$ for 99% of iterative checks
for AES with $\tau=5$ on $n=20$, and even smaller for ASLSQP. Such small
values (17) can be well approximated by those computed with
$\bar{\epsilon}=0$.. We measured the performance of the algorithms by the
number of $f$-calls until they reached $G(x_{t},y_{t})\leqslant 10^{-5}$. For
each setting, we ran $50$ independent trials with different initial search
points generated from $\mathcal{N}(0,I)$. The step-size $\sigma_{x}$ and
$\sigma_{y}$ for AES are initialized to $2$ and their maximal values, denoted
by $\sigma_{\max}$ in Algorithm 2, are also set to $2$.

To confirm (A–C), we ran algorithms without $\eta$ adaptation with
$\eta=\bar{\eta}\times 10^{-k/10}$ for $k=0,1,\dots,15$ on $f_{1}$ with (i)
$n=m=10$ and $b=0.5$, $1$, $2$, $4$, $8$, $16$ and (ii) $n=m=5$, $10$, $20$,
$40$, $80$ and $b=1$. To confirm (D), we ran algorithms with $\eta$ adaptation
on the same problems. The results are summarized in Figure 1.

(A) No run succeeded with $\eta=\bar{\eta}$ for any $n$ and $b$, as
demonstrated in Figure 1-(a, d) for the case of $n=m=10$ and $b=1$, whereas
all runs succeeded with $\eta=\bar{\eta}\times 10^{-1/10}$. This demonstrates
the tightness of $\bar{\eta}$.

(B) For all tested cases, $\eta=\bar{\eta}\times 10^{-3/10}\approx\eta^{*}$
showed the best performance, as demonstrated in Figure 1-(a–f). That is,
$\eta^{*}$ in (17) is a good estimate of the best $\eta$.

(C) In Figure 1-(c, f), we observe that the best $\eta$ was approximately
$\bar{\eta}\times 10^{-3/10}\approx\eta^{*}$ for all $b$ values. The number of
$f$-calls was proportional to
$1/\log(1+b^{-2})\in\mathcal{O}(1/\log(1/\bar{\gamma}^{*}))$, implying that
our expectation of the dependency of $\bar{\sigma}=b$ on $\gamma^{*}$ is
accurate. Moreover, as mentioned above, we observe that all runs succeeded
with $\eta=\bar{\eta}\times 10^{-1/10}$, whereas no run was solved with
$\eta=\bar{\eta}$. We conclude that $\bar{\eta}$, $\eta^{*}$, and
$\bar{\gamma}^{*}$ accurately reflect their dependent value $\bar{\sigma}$ as
expected.

(D) Figure 1-(a–f) shows that AES and ASLSQP with $\eta$ adaptation solved the
problems with up to three times more $\\#f$-calls than AES and ASLSQP with the
best fixed $\eta=\bar{\eta}\times 10^{-3/10}$ for upper quartile cases. We
emphasize that for AES without $\eta$ adaptation to achieve performance
comparable to AES with $\eta$ adaptation, one needs to find
$\eta\in[\bar{\eta}/10,\bar{\eta})$. It is rather difficult to find such
intervals within a few trials when $\bar{\sigma}$ is unknown.

Theorems 4.1 and 4.2 indicate no dependency of $n$ and $m$ on $\gamma$. This
implies that the number of oracle calls does not scale with $n$ and $m$,
whereas the number of $f$-calls may depend on them. In ASLSQP (Figure 1-(e)),
$\\#f$-calls is nearly constant for varying dimensions as its runtime does not
scale with dimension. (note that ASLSQP also requires a similar amount of
$\\#\nabla f$-calls), whereas AES (Figure 1-(b)) requires $\\#f$-calls
proportional to the dimension, as its runtime is proportional.

###  6.2. Ex 2: Convex–concave Case with $\delta>1$

Next, we demonstrate that $\delta<1$ in Theorem 4.1 is not necessary for
convergence. We consider the following function
$f:\mathbb{R}^{m}\times\mathbb{R}^{n}\to\mathbb{R}$ with $m=n=10$:

$f_{2}(x,y)=f_{1}(x,y)-\exp(-\lVert x\rVert^{2}/2)+\exp(-\lVert
y\rVert^{2}/2)$

On such a function, $H_{x,x}=(1+(1-x^{2})\exp(-x^{2}/2))I$,
$-H_{y,y}=(1+(1-y^{2})\exp(-y^{2}/2))I$, $H_{x,y}=H_{y,x}=bI$ and
$H_{x,x}^{*}=-H_{y,y}^{*}=2I$. Then, $G_{x,x}^{*}=G_{y,y}^{*}=((4+b^{2})/2)I$.
We have $\delta>b/2$. Therefore, $b=10$ violates the condition for Theorem
4.1. For $b=10$, ignoring the effect of $\delta$ (i.e., considering local
convergence), we have $\bar{\eta}=4/(4+b^{2})$ and $\eta^{*}=2/(4+b^{2})$ if
$\bar{\epsilon}=0$. We conducted experiments with the same setting as Section
6.1. We measure the progress by the approximate suboptimality error
$\tilde{G}$ around the saddle point.

Figure 2. The number of $f$-calls until $\tilde{G}(x_{t},y_{t})\leqslant
10^{-5}$ is obtained with AES (top) and ASLSQP (bottom) with
$\eta=\bar{\eta}\times 10^{-k/10}$ for $k=0,\dots,15$ and $\eta$ adaptation on
$f_{2}$ with $b=10$ and $n=m=10$.

Figure 2 shows the boxplot of the number of $f$-calls until
$\tilde{G}(x_{t},y_{t})\leqslant 10^{-5}$ is reached. As suggested in Section
4.3, both algorithms were able to locate a nearly optimal saddle point with a
sufficiently small $\eta$. This reveals room for improvement in Theorem 4.1.
We conjecture that the upper bound on $\eta$ to guarantee linear convergence
will be similar to those obtained in (Liang and Stokes, 2019) for the
simultaneous gradient update.

###  6.3. Ex 3: NonConvex–concave Case

A function with nonoptimal critical points

$f_{3}(x,y)=2x^{2}+4xy+y^{2}+(4/3)y^{3}-(1/4)y^{4}$

has three critical points $z_{0}=(0,0)$, $z_{1}=(-2-\sqrt{2},2+\sqrt{2})$, and
$z_{2}=(-2+\sqrt{2},2-\sqrt{2})$, of which only $z_{1}$ is the locally optimal
saddle point, whereas the others are local minima. The simultaneous gradient
method (2) is reportedly attracted by local minima (Adolphs et al., 2019).
ASLSQP is also gradient based, but the update step is different from (2), and
it performs $\tau_{\mathrm{slsqp}}$ steps for each oracle call. A question
arises as to whether the same undesired convergence is observed for ASLSQP. To
investigate this question, we run ASLSQP with $\eta=0.1$ and
$\tau_{\mathrm{slsqp}}=1$ (single step) and $5$ (multi-step) from $51^{2}$
different initial points on $[-5,3]\times[-3,5]$. The progress is measured by
$\tilde{G}(x,y)$ around the saddle point $z_{1}=(x^{*},y^{*})$. The results
are shown in Figure 3. All the trials succeeded in locating an approximate
saddle point. From this result, we conjecture that nonsaddle critical points
are not attractors of ASLSQP, which is, however, not covered by Theorem 4.2.
Future work in this line is required.

(a) $\eta=0.1$, $\tau_{\mathrm{slsqp}}=1$

(b) $\eta=0.1$, $\tau_{\mathrm{slsqp}}=5$

Figure 3. $\\#f$-calls until $\tilde{G}(x_{t},y_{t})\leqslant 10^{-5}$ is
obtained from each initial point. The three circles are $z_{1}$, $z_{2}$ and
$z_{3}$.

###  6.4. Ex 4: Limitation of Existing Coevolutionary Approaches

Typical existing coevolutionary approaches do not fit into our framework (10).
Therefore, the limitations of existing coevolutionary approaches are not
revealed in our analysis. To show the nonconvergent behavior of these
approaches, we run mmde (Qiu et al., 2018), coeva, and reckless (Al-Dujaili et
al., 2019) on $f_{1}$ on $m=n=2,10$ with $b=0,2$. We used the implementation
of these three approaches published by the authors of (Al-Dujaili et al.,
2019). Because their approaches assume a rectangular search space $[0,1]^{m}$
and $[0,1]^{n}$, the function is shifted such that the optimal saddle point is
at $x^{*}=(0.5,\cdots,0.5)$ and $y^{*}=(0.5,\cdots,0.5)$. We also run AES with
$\eta$-adaptation, where the initial solution is generated uniform-randomly in
$[0,1]^{m}\times[0,1]^{n}$, $\sigma_{x}$ and $\sigma_{y}$ are initialized to
$1/4$, and their maximum values are set to $1$. The box constraint is treated
with the mirroring technique, i.e., a solution outside the feasible domain is
mapped by applying the transformation $x\mapsto
1-\lvert\text{mod}(x,2)-1\rvert$ for each coordinate.

The results are shown in Figure 4. AES converged linearly towards the global
min–max saddle point. Except for reckless on $b=0$, the existing approaches,
reckless, mmde, and coeva, failed to converge. The case $b=0$ means that there
is no interaction between $x$ and $y$, and these variables can be optimized
separately. Even in such an easy situation, mmde and coeva failed to reach the
target quality within $5\times 10^{5}$ $f$-calls. For cases with an
interaction term ($b>0$), no existing algorithms reached the target quality
within a given budget, whereas AES reached it with $f$-calls less than
$10^{5}$. This demonstrates the difference in convergence behavior between the
proposed framework and existing coevolutionary approaches.

(a) $b=0$, $n=m=2$

(b) $b=2$, $n=m=2$

(c) $b=0$, $n=m=10$

(d) $b=2$, $n=m=10$

Figure 4. 5 runs of existing coevolutionary approaches on $f_{1}$

##  7\. Summary and Final Remarks

A saddle point optimization (10) based solely on an approximate minimization
oracle is studied. Its convergence properties are analyzed in Theorems 4.1 and
4.2, where a sufficient condition on the learning rate $\eta$ to exhibit
linear (geometric) convergence of an approximate suboptimality error is
derived. One important remark is that simultaneous minimization approaches
((10) with $\eta=1$) do not satisfy the sufficient condition. This may reveal
a shortcoming of existing approaches employing simultaneous or alternating
updates with $\eta=1$ (Pinto et al., 2017; Shioya et al., 2018; Al-Dujaili et
al., 2019). As discussed in Section 4.3, this sufficient condition is also
necessary for a convex–concave quadratic function. Hence, it implies that they
do not converge in such a situation and suggests a modification to these
algorithms. A learning rate adaptation heuristic is then proposed and
evaluated on test problems.

The generality of our algorithmic framework (10) is one of its advantages over
analyses on specific algorithms. The minimization oracle can be zero-order,
first-order, etc. Moreover, it can share internal information during different
oracle calls, as it does in AES. The latter point may be useful for
understanding the effect of introducing target networks in deep actor–critic
algorithms (Lillicrap et al., 2016), which may be regarded as (10), where each
oracle call starts from the final solution of the last oracle call. Further
investigation in this direction may be an interesting research topic.

We used the (1+1)-ES algorithm as the minimization oracle in AES for its
simplicity to demonstrate the main idea and its theoretical aspects. For
practical use, it is generally advised to use the CMA-ES (Hansen et al., 2003;
Hansen and Auger, 2014; Hansen and Ostermeier, 2001; Hansen and Kern, 2004;
Akimoto and Hansen, 2020) as the oracle $M_{\epsilon_{x},r_{x}}^{A_{x}}$ in
AES instead of the (1+1)-ES. We would like to emphasize that the difficulty
addressed in this paper is associated with the interaction term between $x$
and $y$ in $f$, i.e., $H_{x,y}$ and $H_{y,x}$, and the CMA-ES does not help
with this. Using the CMA-ES helps to solve poor conditioning of $H_{x,x}$ and
$H_{y,y}$, which is not considered in this study.

As an important application of our approach, we are interested in a robust
design of a feedback controller for automatic berthing. When we design a
controller, we first model the state equation of a control target, which is a
ship in this case, and optimize the parameters of the controller on
simulation. However, there always exist modelling error and uncertainty in
environment. For example, coefficients of a state equation model are typically
estimated by water tank tests and weather conditions are uncertain when
modelling. To design a controller reliable in a real environment, a task is
formulated as a simulation-based min–max optimization problem (1). We apply
our approach to robust automatic berthing to demonstrate the usefulness of the
proposed approach. The results will be reported in another occasion.

We close the paper with open questions on this topic. (I) A sufficient
condition on $\eta$ derived in Theorem 4.1 for global linear convergence is
not tight for $\delta>0$ (see Section 6). We conjecture that we may be able to
obtain a condition similar to that obtained for a simultaneous gradient update
(Mescheder et al., 2017). (II) The current analysis does not guarantee that
our approach does not converge to nonsaddle critical points, whereas our
experiments in Section 6 demonstrated that it does not. (III) The proposed
learning rate adaptation is yet to be theoretically analyzed for its
convergence. (IV) This study did not investigate methods for improving the
convergence rate $\gamma$ when the interaction terms $H_{x,y}$ and $H_{y,x}$
are ill conditioned. The proposed approach introduced a learning rate $\eta<1$
to guarantee convergence, but the resulting convergence rate depends on
$\eta$. (V) There can exist a local (or global) optimum to the min–max
optimization problem (1) that is not a min–max saddle point. This corresponds
to the case where the optimum is not a critical point of any of $y$. Though it
is out of the scope of the current paper, developing an algorithm that can
locate such a local optimum of (1) is an important direction of the future
work.

###### Acknowledgements.

The author would like to thank the authors of (Al-Dujaili et al., 2019) for
publishing their code. This work is partially supported by JSPS KAKENHI Grant
Number 19H04179.

## References

  * (1)
  * Adolphs et al. (2019) Leonard Adolphs, Hadi Daneshmand, Aurelien Lucchi, and Thomas Hofmann. 2019.  Local Saddle Point Optimization: A Curvature Exploitation Approach. In International Conference on Artificial Intelligence and Statistics. 486–495.
  * Akimoto et al. (2018) Youhei Akimoto, Anne Auger, and Tobias Glasmachers. 2018\.  Drift theory in continuous search spaces: expected hitting time of the (1+1)-ES with 1/5 success rule. In Genetic and Evolutionary Computation Conference. 801–808.
  * Akimoto and Hansen (2020) Youhei Akimoto and Nikolaus Hansen. 2020.  Diagonal acceleration for covariance matrix adaptation evolution strategies.  Evolutionary computation 28, 3 (2020), 405–435.
  * Al-Dujaili et al. (2019) Abdullah Al-Dujaili, Shashank Srikant, Erik Hemberg, and Una-May O’Reilly. 2019.  On the application of Danskin’s theorem to derivative-free minimax problems.  AIP Conference Proceedings 2070, 1 (2019), 20–26.
  * Bogunovic et al. (2018) Ilija Bogunovic, Jonathan Scarlett, Stefanie Jegelka, and Volkan Cevher. 2018.  Adversarially Robust Optimization with Gaussian Processes.  In Advances in Neural Information Processing Systems. 5760–5770.
  * Branke and Rosenbusch (2008) Jürgen Branke and Johanna Rosenbusch. 2008.  New Approaches to Coevolutionary Worst-Case Optimization. In International Conference on Parallel Problem Solving from Nature. 144–153.
  * Cherukuri et al. (2017) Ashish Cherukuri, Bahman Gharesifard, and Jorge Cortés. 2017\.  Saddle-Point Dynamics: Conditions for Asymptotic Stability of Saddle Points.  SIAM Journal on Control and Optimization 55, 1 (2017), 486–511.
  * Conn and Vicente (2012) A. R. Conn and L. N. Vicente. 2012.  Bilevel Derivative-Free Optimization and Its Application to Robust Optimization.  Optimization Methods Software 27, 3 (2012), 561–577.
  * Daskalakis et al. (2018) Constantinos Daskalakis, Andrew Ilyas, Vasilis Syrgkanis, and Haoyang Zeng. 2018.  Training GANs with Optimism. In International Conference on Learning Representations.
  * de Oliveira (2013) Oswaldo de Oliveira. 2013\.  The Implicit and Inverse Function Theorems: Easy Proofs.  Real Anal. Exchange 39, 1 (2013), 207–218.
  * Devroye (1972) L. Devroye. 1972\.  The compound random search. In International Symposium on Systems Engineering and Analysis. 195–110.
  * Gidel et al. (2017) Gauthier Gidel, Tony Jebara, and Simon Lacoste-Julien. 2017\.  Frank-Wolfe Algorithms for Saddle Point Problems. In International Conference on Artificial Intelligence and Statistics. 362–371.
  * Goodfellow et al. (2014) Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio. 2014\.  Generative adversarial nets. In Advances in neural information processing systems. 2672–2680.
  * Hansen and Auger (2014) Nikolaus Hansen and Anne Auger. 2014.  Principled design of continuous stochastic search: From theory to practice.  In Theory and principled methods for the design of metaheuristics. Springer, 145–180.
  * Hansen and Kern (2004) Nikolaus Hansen and Stefan Kern. 2004.  Evaluating the CMA evolution strategy on multimodal test functions. In International Conference on Parallel Problem Solving from Nature. Springer, 282–291.
  * Hansen et al. (2003) Nikolaus Hansen, Sibylle D Müller, and Petros Koumoutsakos. 2003\.  Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES).  Evolutionary computation 11, 1 (2003), 1–18.
  * Hansen and Ostermeier (2001) Nikolaus Hansen and Andreas Ostermeier. 2001.  Completely derandomized self-adaptation in evolution strategies.  Evolutionary computation 9, 2 (2001), 159–195.
  * Jensen (2004) Mikkel T. Jensen. 2004\.  A New Look at Solving Minimax Problems with Coevolutionary Genetic Algorithms.  Kluwer Academic Publishers, 369–384.
  * Kern et al. (2004) S. Kern, S. D. Müller, N. Hansen, D. Büche, J. Ocenasek, and P. Koumoutsakos. 2004\.  Learning probability distributions in continuous evolutionary algorithms–a comparative review.  Natural Computing 3, 1 (2004), 77–112.
  * Kraft (1988) Dieter Kraft. 1988\.  A software package for sequential quadratic programming.  Technical Report.
  * Liang and Stokes (2019) Tengyuan Liang and James Stokes. 2019.  Interaction matters: A note on non-asymptotic local convergence of generative adversarial networks. In International Conference on Artificial Intelligence and Statistics. 907–915.
  * Lillicrap et al. (2016) Timothy P Lillicrap, Jonathan J Hunt, Alexander Pritzel, Nicolas Heess, Tom Erez, Yuval Tassa, David Silver, and Daan Wierstra. 2016.  Continuous control with deep reinforcement learning.  In International Conference on Learning Representations.
  * Liu et al. (2020) Sijia Liu, Songtao Lu, Xiangyi Chen, Yao Feng, Kaidi Xu, Abdullah Al-Dujaili, Mingyi Hong, and Una-May O’Reilly. 2020\.  Min-Max Optimization without Gradients: Convergence and Applications to Black-Box Evasion and Poisoning Attacks. In International Conference on Machine Learning. 2307–2318.
  * Mescheder et al. (2017) Lars Mescheder, Sebastian Nowozin, and Andreas Geiger. 2017\.  The Numerics of GANs. In Advances in Neural Information Processing Systems. 1823–1833.
  * Morinaga and Akimoto (2019) Daiki Morinaga and Youhei Akimoto. 2019.  Generalized drift analysis in continuous domain: linear convergence of (1+1)-ES on strongly convex functions with Lipschitz continuous gradients. In Foundations of Genetic Algorithms. 13–24.
  * Morinaga et al. (2021) Daiki Morinaga, Kazuto Fukuchi, Jun Sakuma, and Youhei Akimoto. 2021.  Convergence Rate of the (1+1)-Evolution Strategy with Success-Based Step-Size Adaptation on Convex Quadratic Function. In Genetic and Evolutionary Computation Conference.  (Accepted as a full paper).
  * Nagarajan and Kolter (2017) Vaishnavh Nagarajan and J. Zico Kolter. 2017.  Gradient Descent GAN Optimization is Locally Stable. In Advances in Neural Information Processing Systems. 5591–5600.
  * Nouiehed et al. (2019) Maher Nouiehed, Maziar Sanjabi, Tianjian Huang, Jason D Lee, and Meisam Razaviyayn. 2019.  Solving a Class of Non-Convex Min-Max Games Using Iterative First Order Methods.  In Advances in Neural Information Processing Systems. 14934–14942.
  * Picheny et al. (2019) Victor Picheny, Mickael Binois, and Abderrahmane Habbal. 2019\.  A Bayesian optimization approach to find Nash equilibria.  Journal of Global Optimization 73, 1 (2019), 171–192.
  * Pinto et al. (2017) Lerrel Pinto, James Davidson, Rahul Sukthankar, and Abhinav Gupta. 2017.  Robust Adversarial Reinforcement Learning. In International Conference on Machine Learning. 2817–2826.
  * Qiu et al. (2018) X. Qiu, J. Xu, Y. Xu, and K. C. Tan. 2018\.  A New Differential Evolution Algorithm for Minimax Optimization in Robust Design.  IEEE Transactions on Cybernetics 48, 5 (2018), 1355–1368.
  * Rechenberg (1973) Ingo Rechenberg. 1973\.  Evolutionsstrategie: Optimierung technisher Systeme nach Prinzipien der biologischen Evolution.  Frommann-Holzboog.
  * Salimans et al. (2016) Tim Salimans, Ian Goodfellow, Wojciech Zaremba, Vicki Cheung, Alec Radford, Xi Chen, and Xi Chen. 2016.  Improved Techniques for Training GANs.  In Advances in Neural Information Processing Systems. 2234–2242.
  * Schumer and Steiglitz (1968) M. Schumer and K. Steiglitz. 1968.  Adaptive step size random search.  Automatic Control, IEEE Transactions on 13 (1968), 270–276.
  * Shioya et al. (2018) Hiroaki Shioya, Yusuke Iwasawa, and Yutaka Matsuo. 2018\.  Extending Robust Adversarial Reinforcement Learning Considering Adaptation and Diversity. In International Conference on Learning Representations, Workshop Track Proceedings.
  * Yadav et al. (2018) Abhay Yadav, Sohil Shah, Zheng Xu, David Jacobs, and Tom Goldstein. 2018.  Stabilizing Adversarial Nets With Prediction Methods. In International Conference on Learning Representations.

##  Appendix A Proofs

The following proposition is derived immediately from the implicit function
theorem, Theorem 5 of (de Oliveira, 2013). The proofs of Theorems 4.1 and 4.2
are based on this result.

###### Proposition A.1 (Implicit Function Theorem).

In the following, let $(x^{*},y^{*})$ be a strict saddle point of $f$ and $f$
be locally strong convex–concave around $(x^{*},y^{*})$. There exist open sets
$\mathcal{E}_{x,x}\subseteq\mathbb{R}^{m}$ including $x^{*}$ and
$\mathcal{E}_{x,y}\subseteq\mathbb{R}^{n}$ including $y^{*}$, such that there
is a unique $\hat{y}:\mathcal{E}_{x,x}\to\mathcal{E}_{x,y}$ such that
$\nabla_{y}f(x,\hat{y}(x))=0$. Moreover, $y^{*}=\hat{y}(x^{*})$ and
$J_{\hat{y}}(x)=-(H_{y,y}(x,\hat{y}(x)))^{-1}H_{y,x}(x,\hat{y}(x))$ for all
$x\in\mathcal{E}_{x,x}$. Analogously, there exist open sets
$\mathcal{E}_{y,y}\subseteq\mathbb{R}^{n}$ including $y^{*}$ and
$\mathcal{E}_{y,x}\subseteq\mathbb{R}^{m}$ including $x^{*}$, such that there
is a unique $\hat{x}:\mathcal{E}_{y,y}\to\mathcal{E}_{y,x}$ such that
$\nabla_{x}f(\hat{x}(y),y)=0$. Moreover, $x^{*}=\hat{x}(y^{*})$ and
$J_{\hat{x}}(y)=-(H_{x,x}(\hat{x}(y),y))^{-1}H_{x,y}(\hat{x}(y),y)$ for all
$y\in\mathcal{E}_{y,y}$. If $f$ is globally strong convex–concave, one can
take $\mathcal{E}_{x,x}=\mathcal{E}_{y,x}=\mathbb{R}^{m}$ and
$\mathcal{E}_{y,y}=\mathcal{E}_{x,y}=\mathbb{R}^{n}$ on the above statements.

Proposition A.1 shows that under the condition stated in the proposition, the
local minimal solutions to
$\operatornamewithlimits{argmin}_{x^{\prime}\in\mathbb{R}^{m}}f(x^{\prime},y)$
and
$\operatornamewithlimits{argmin}_{y^{\prime}\in\mathbb{R}^{n}}-f(x,y^{\prime})$
are uniquely determined by $\hat{x}(y)$ and $\hat{y}(x)$. Moreover, by the
mean value theorem we obtain

(18)
$\hat{x}(y)-x^{*}=-(H_{x,x}(\hat{x}(\tilde{y}),\tilde{y}))^{-1}H_{x,y}(\hat{x}(\tilde{y}),\tilde{y})(y-y^{*})\enspace,$

where $\tilde{y}=(1-b)y+by^{*}$ for some $b\in(0,1)$. The approximate
minimization oracle returns its approximate solutions $x^{\prime}$ and
$y^{\prime}$. Then, we can write

(19) $x_{t+1}-x^{*}=(1-\eta)(x_{t}-x^{*})\\\
+\eta(x_{t}^{\prime}-\hat{x}(y_{t}))+\eta(\hat{x}(y_{t})-x^{*})\enspace.$

Analogously we can obtain the decomposition of $y_{t+1}-y^{*}$.

###  A.1. Proof of Theorem 4.1

Let $x^{\prime}=M^{A_{x}}_{\epsilon_{x},r_{x}}(f(\cdot,y_{t}),x_{t})$ and
$y^{\prime}=M^{A_{y}}_{\epsilon_{y},r_{y}}(-f(x_{t},\cdot),y_{t})$. Let
$\hat{x}(y):\mathbb{R}^{n}\to\mathbb{R}^{m}$ and
$\hat{y}(x):\mathbb{R}^{m}\to\mathbb{R}^{n}$ be as shown in Proposition A.1.
We are going to investigate each term of (19).

First, we study $x^{\prime}-\hat{x}(y_{t})$ and $y^{\prime}-\hat{y}(x_{t})$,
which are the deviations of the outputs of the approximate minimization
oracles from the exact solutions. For
$(x_{t},y_{t})\in\mathbb{R}^{m}\times\mathbb{R}^{n}$, we can write
$\sqrt{\smash[b]{G_{x,x}^{*}}}(x^{\prime}-\hat{x}(y_{t}))=E_{x}\sqrt{\smash[b]{G_{x,x}^{*}}}(x_{t}-\hat{x}(y_{t}))$
and
$\sqrt{\smash[b]{G_{y,y}^{*}}}(y^{\prime}-\hat{y}(x_{t}))=E_{y}\sqrt{\smash[b]{G_{y,y}^{*}}}(y_{t}-\hat{y}(x_{t}))$,
where $E_{x}$ and $E_{y}$ are square matrices defined as

$\begin{split}E_{x}&=\sqrt{\smash[b]{G_{x,x}^{*}}}(x^{\prime}-\hat{x}(y_{t}))(x_{t}-\hat{x}(y_{t}))^{\mathrm{T}}\sqrt{\smash[b]{G_{x,x}^{*}}}/\lVert
x_{t}-\hat{x}(y_{t})\rVert_{G_{x,x}^{*}}^{2}\\\
E_{y}&=\sqrt{\smash[b]{G_{y,y}^{*}}}(y^{\prime}-\hat{y}(x_{t}))(y_{t}-\hat{y}(x_{t}))^{\mathrm{T}}\sqrt{\smash[b]{G_{y,y}^{*}}}/\lVert
y_{t}-\hat{y}(x_{t})\rVert_{G_{y,y}^{*}}^{2}\enspace.\end{split}$

From (9), we have

(20) $\begin{split}&\textstyle\lVert
x^{\prime}-\hat{x}(y_{t})\rVert_{G_{x,x}^{*}}^{2}\leqslant\epsilon_{x}\lVert
x_{t}-\hat{x}(y_{t})\rVert_{G_{x,x}^{*}}^{2}\enspace,\\\ &\textstyle\lVert
y^{\prime}-\hat{y}(x_{t})\rVert_{G_{y,y}^{*}}^{2}\leqslant\epsilon_{y}\lVert
y_{t}-\hat{y}(x_{t})\rVert_{G_{y,y}^{*}}^{2}\enspace.\end{split}$

Then, from (10) and (20), we have

$\begin{split}\sigma(E_{x})&=\lVert
x^{\prime}-\hat{x}(y_{t})\rVert_{G_{x,x}^{*}}^{2}/\lVert
x_{t}-\hat{x}(y_{t})\rVert_{G_{x,x}^{*}}^{2}\leqslant\epsilon_{x}\\\
\sigma(E_{y})&=\lVert y^{\prime}-\hat{y}(x_{t})\rVert_{G_{y,y}^{*}}^{2}/\lVert
y_{t}-\hat{y}(x_{t})\rVert_{G_{y,y}^{*}}^{2}\leqslant\epsilon_{y}\enspace.\end{split}$

Next, we study $\hat{x}(y_{t})-x^{*}$ and $\hat{y}(x_{t})-y^{*}$, which are
the differences between the saddle point $(x^{*},y^{*})$ and the solutions to
the minimization problems that the approximate minimization oracles are
solving. Let

$\displaystyle\Delta_{x,y}=\textstyle(H_{x,x}(\hat{x}(\tilde{y}),\tilde{y}))^{-1}H_{x,y}(\hat{x}(\tilde{y}),\tilde{y})-(H_{x,x}^{*})^{-1}H_{x,y}^{*},$

where $\tilde{y}$ is as defined below (18). From (18), we can write

$\displaystyle\hat{x}(y_{t})-x^{*}=-((H_{x,x}^{*})^{-1}H_{x,y}^{*}+\Delta_{x,y})\cdot(y_{t}-y^{*})\enspace.$

Then we obtain

$\sqrt{\smash[b]{G_{x,x}^{*}}}(\hat{x}(y_{t})-x^{*})=-\Big{(}\sqrt{\smash[b]{G_{x,x}^{*}}}(H_{x,x}^{*})^{-1}H_{x,y}^{*}\sqrt{\smash[b]{G_{y,y}^{*}}}^{-1}\\\
+\sqrt{\smash[b]{G_{x,x}^{*}}}\Delta_{x,y}\sqrt{\smash[b]{G_{y,y}^{*}}}^{-1}\Big{)}\cdot\sqrt{\smash[b]{G_{y,y}^{*}}}(y_{t}-y^{*})\enspace.$

Analogously, in light of Proposition A.1, we have

$\displaystyle\hat{y}(x_{t})-y^{*}=-((H_{y,y}^{*})^{-1}H_{y,x}^{*}+\Delta_{y,x})\cdot(x_{t}-x^{*})\enspace,$

where

$\displaystyle\Delta_{y,x}$
$\displaystyle=\textstyle(H_{y,y}(\hat{y}(\tilde{x}),\tilde{x}))^{-1}H_{y,x}(\hat{y}(\tilde{x}),\tilde{x})-(H_{y,y}^{*})^{-1}H_{y,x}^{*},$

where $\tilde{x}=(1-a)x+ax^{*}$ for some $a\in(0,1)$. We obtain

$\sqrt{\smash[b]{G_{y,y}^{*}}}(\hat{y}(x_{t})-y^{*})=-\Big{(}\sqrt{\smash[b]{G_{y,y}^{*}}}(H_{y,y}^{*})^{-1}H_{y,x}^{*}\sqrt{\smash[b]{G_{x,x}^{*}}}^{-1}\\\
+\sqrt{\smash[b]{G_{y,y}^{*}}}\Delta_{y,x}\sqrt{\smash[b]{G_{x,x}^{*}}}^{-1}\Big{)}\cdot\sqrt{\smash[b]{G_{x,x}^{*}}}(x_{t}-x^{*})\enspace.$

Because of the decomposition (19) and the same decomposition for $y$, i.e.,

(21) $y_{t+1}-y^{*}=(1-\eta)(y_{t}-y^{*})\\\
+\eta(y_{t}^{\prime}-\hat{y}(x_{t}))+\eta(\hat{y}(x_{t})-y^{*})\enspace,$

we have

(22) $\begin{bmatrix}\sqrt{\smash[b]{G_{x,x}^{*}}}(x_{t+1}-x^{*})\\\
\sqrt{\smash[b]{G_{y,y}^{*}}}(y_{t+1}-y^{*})\end{bmatrix}=(F+\eta\tilde{R}_{1}+\eta\tilde{R}_{2})\begin{bmatrix}\sqrt{\smash[b]{G_{x,x}}}(x_{t}-x^{*})\\\
\sqrt{\smash[b]{G_{y,y}}}(y_{t}-y^{*})\end{bmatrix},$

where

$\displaystyle F$
$\displaystyle=\begin{bmatrix}(1-\eta)I\quad-\eta\sqrt{\smash[b]{G_{x,x}^{*}}}(H_{x,x}^{*})^{-1}H_{x,y}^{*}\sqrt{\smash[b]{G_{y,y}^{*}}}^{-1}\\\
-\eta\sqrt{\smash[b]{G_{y,y}^{*}}}(H_{y,y}^{*})^{-1}H_{y,x}^{*}\sqrt{\smash[b]{G_{x,x}^{*}}}^{-1}\quad(1-\eta)I\end{bmatrix},$
$\displaystyle\tilde{R}_{1}$ $\displaystyle=\begin{bmatrix}E_{x}&0\\\
0&E_{y}\end{bmatrix},\
\tilde{R}_{2}=\begin{bmatrix}0\quad-\sqrt{\smash[b]{G_{x,x}^{*}}}\Delta_{x,y}\sqrt{\smash[b]{G_{y,y}^{*}}}^{-1}\\\
-\sqrt{\smash[b]{G_{y,y}^{*}}}\Delta_{y,x}\sqrt{\smash[b]{G_{x,x}^{*}}}^{-1}\quad
0\end{bmatrix}.$

Now, we show that there exists $\gamma>0$ such that for any
$(x_{t},y_{t})\in\mathbb{R}^{m}\times\mathbb{R}^{n}$,

(23) $\lVert x_{t+1}-x^{*}\rVert_{G_{x,x}^{*}}^{2}+\lVert
y_{t+1}-y^{*}\rVert_{G_{y,y}^{*}}^{2}\\\ \leqslant\gamma^{2}\big{(}\lVert
x_{t}-x^{*}\rVert_{G_{x,x}^{*}}^{2}+\lVert
y_{t}-y^{*}\rVert_{G_{y,y}^{*}}^{2}\big{)}\enspace.$

Notably, the left-hand side of (23) is the squared Euclidean norm of the left-
hand side of (22). Therefore, to prove (23), it suffices to show that the
greatest singular value of $F+\eta\tilde{R}_{1}+\eta\tilde{R}_{2}$ has an
upper bound $\gamma$. The greatest singular value of
$F+\eta\tilde{R}_{1}+\eta\tilde{R}_{2}$ has an upper bound of the sum of the
greatest singular value of each term, i.e.,
$\sigma(F+\eta\tilde{R}_{1}+\eta\tilde{R}_{2})\leqslant\sigma(F)+\eta\sigma(\tilde{R}_{1})+\eta\sigma(\tilde{R}_{2})$.
The greatest singular value of $F$ is equal to the square root of the greatest
eigenvalue of $FF^{\mathrm{T}}$ in this case. Noting that
$\sqrt{\smash[b]{G_{x,x}}}(H_{x,x}^{*})^{-1}H_{x,y}^{*}\sqrt{\smash[b]{G_{y,y}}}^{-1}=\sqrt{\smash[b]{G_{x,x}}}^{-1}H_{x,y}^{*}(-H_{y,y}^{*})^{-1}\sqrt{\smash[b]{G_{y,y}}}$,
we have

$FF^{\mathrm{T}}=(1-\eta)^{2}I+\eta^{2}\times\\\
\begin{bmatrix}\sqrt{\smash[b]{G_{x,x}}}(H_{x,x}^{*})^{-1}H_{x,y}^{*}G_{y,y}^{-1}H_{y,x}^{*}(H_{x,x}^{*})^{-1}\sqrt{\smash[b]{G_{x,x}}}\qquad
0\\\
0\qquad\sqrt{\smash[b]{G_{y,y}}}(H_{y,y}^{*})^{-1}H_{y,x}^{*}G_{x,x}^{-1}H_{x,y}^{*}(H_{y,y}^{*})^{-1}\sqrt{\smash[b]{G_{y,y}}}\end{bmatrix}\enspace,$

whose greatest eigenvalue is $(1-\eta)^{2}+\eta^{2}\bar{\sigma}^{2}$. Hence,
$\sigma(F)=((1-\eta)^{2}+\eta^{2}\bar{\sigma}^{2})^{1/2}$. The greatest
singular value of $\tilde{R}_{1}$ has an upper bound of
$\bar{\epsilon}=\max(\epsilon_{x},\epsilon_{y})$, as derived below (20). The
greatest singular value of $\tilde{R}_{2}$ has an upper bound of

$\max\big{[}\sigma\big{(}\sqrt{\smash[b]{G_{x,x}^{*}}}\Delta_{x,y}\sqrt{\smash[b]{G_{y,y}^{*}}}^{-1}\big{)},\sigma\big{(}\sqrt{\smash[b]{G_{y,y}^{*}}}\Delta_{y,x}\sqrt{\smash[b]{G_{x,x}^{*}}}^{-1}\big{)}\big{]}\leqslant\delta\enspace.$

Therefore,
$\gamma\leqslant((1-\eta)^{2}+\eta^{2}\bar{\sigma}^{2})^{1/2}+\eta(\bar{\epsilon}+\delta)<1$.

Because (23) holds for any
$(x_{t},y_{t})\in\mathbb{R}^{m}\times\mathbb{R}^{n}$, a repeated application
of (23) shows (14). This completes the proof.

###  A.2. Proof of Theorem 4.2

First, we show that there exists $\bar{r}_{x}>0$ and $\bar{r}_{y}>0$ such that
for any $r_{x}\in(0,\bar{r}_{x})$ and $r_{y}\in(0,\bar{r}_{y})$, there exists
$U_{x}$ and $U_{y}$ such that
$\hat{y}(x)=\operatornamewithlimits{argmin}_{y^{\prime}\in
U_{A_{y}}(y,r_{y})}-f(x,y^{\prime})$ and
$\hat{x}(y)=\operatornamewithlimits{argmin}_{x^{\prime}\in
U_{A_{x}}(x,r_{x})}f(x^{\prime},y)$ hold for all $(x,y)\in U_{x}\times U_{y}$,
where $\hat{y}:U_{x}\to U_{A_{y}}(y,r_{y})$ and $\hat{x}:U_{y}\to
U_{A_{x}}(x,r_{x})$ are as defined in Proposition A.1. In light of Proposition
A.1, there exists a unique $\hat{x}:\mathcal{E}_{x,x}\to\mathcal{E}_{x,y}$ and
a unique $\hat{y}:\mathcal{E}_{y,y}\to\mathcal{E}_{y,x}$ for some
$\mathcal{E}_{x,x}$, $\mathcal{E}_{x,y}$, $\mathcal{E}_{y,x}$, and
$\mathcal{E}_{y,y}$. Because $A_{x}=G_{x,x}^{*}$, we have
$U_{A_{x}}(x,r_{x})\subseteq\\{x^{\prime}\in\mathbb{R}^{m}:\lVert
x^{\prime}-x\rVert_{G_{x,x}^{*}}\leqslant r_{x}\\}$. Hence, there exists a
$\bar{r}_{x}>0$ such that for all $r_{x}\in(0,\bar{r}_{x})$ we have
$U_{A_{x}}(x^{*},r_{x})\subseteq\mathcal{E}_{y,x}$. Analogously, we have that
there exists a $\bar{r}_{y}>0$ such that for all $r_{y}\in(0,\bar{r}_{y})$ we
have $U_{A_{y}}(y^{*},r_{y})\subseteq\mathcal{E}_{x,y}$. Because $\hat{x}$ and
$\hat{y}$ are continuous and $\hat{x}(y^{*})=x^{*}$ and
$\hat{y}(x^{*})=y^{*}$, for any $r_{x}\in(0,\bar{r}_{x})$ and
$r_{y}\in(0,\bar{r}_{y})$, there exist $U_{x}\subseteq\mathcal{E}_{x,x}$
including $x^{*}$ and $U_{y}\subseteq\mathcal{E}_{y,y}$ including $y^{*}$ such
that $\hat{x}(y)\in U_{A_{x}}(x,r_{x})$ and $\hat{y}(x)\in U_{A_{y}}(y,r_{y})$
for all $x\in U_{x}$ and $y\in U_{y}$.

For $(x,y)\in U_{x}\times U_{y}$, we obtain (23) by following the line of the
proof of Theorem 4.2, where
$\gamma\leqslant((1-\eta)^{2}+\eta^{2}\bar{\sigma}^{2})^{1/2}+\eta(\bar{\epsilon}+\delta)$
with $\delta$ being the maximum of
$\sigma\big{(}\sqrt{\smash[b]{G_{x,x}^{*}}}\Delta_{x,y}\sqrt{\smash[b]{G_{y,y}^{*}}}^{-1}\big{)}$
and
$\sigma\big{(}\sqrt{\smash[b]{G_{y,y}^{*}}}\Delta_{y,x}\sqrt{\smash[b]{G_{x,x}^{*}}}^{-1}\big{)}$.
Because of the continuities of $H_{x,x}$, $H_{x,y}$, $H_{y,x}$, and $H_{y,y}$,
for any $\bar{\delta}>0$ there exists a neighborhood $\tilde{U}\subseteq
U_{x}\times U_{y}$ of $(x^{*},y^{*})$ such that $\bar{\delta}<\bar{\delta}$
for all $(x,y)\in\tilde{U}$. That is, for any $\gamma\in(\bar{\gamma},1)$,
where
$\bar{\gamma}=((1-\eta)^{2}+\eta^{2}\bar{\sigma}^{2})^{1/2}+\eta\bar{\epsilon}$,
we can find a neighborhood $\tilde{U}$ of $(x^{*},y^{*})$ such that (23) holds
for all $(x_{t},y_{t})\in\tilde{U}$.

Finally, we show that there exists $U\subseteq\tilde{U}$ such that the
sequence $(x_{t},y_{t})$ starting from $(x_{0},y_{0})\in U$ never leaves $U$.
Define
$B(r):=\\{(x,y)\in\mathbb{R}^{m}\times\mathbb{R}^{n}:\tilde{G}(x,y)<r^{2}/2\\}$.
Then, one can find an $\bar{r}>0$ such that $B(\bar{r})\subseteq\tilde{U}$.
Let $U=B(\bar{r})\subseteq\tilde{U}$. Because $\gamma<1$, (23) implies that if
$(x_{0},y_{0})\in U$, we have $(x_{t},y_{t})\in U$. Hence, $(x_{t},y_{t})$
never leaves $U$.

Repeated applications of (23) lead to (16). This completes the proof.