Accelerated Algorithms for Smooth Convex-Concave Minimax Problems with $\mathcal{O}(1/k^2)$ Rate on Squared Gradient Norm

Introduction

Minimax optimization problems, or minimax games, of the form

have recently gained significant interest in the optimization and machine learning communities due to their application in adversarial training (Goodfellow et al., 2015; Madry et al., 2018) and generative adversarial networks (GANs) (Goodfellow et al., 2014).

Prior works on minimax optimization often consider compact domains $X,Y$ for $\mathbf{x},\mathbf{y}$ and use the duality gap

to quantify suboptimality of algorithms’ iterates in solving (1). However, while it is a natural analog of minimization error for minimax problems, the duality gap can be difficult to measure directly in practice, and it is unclear how to generalize the notion to non-convex-concave problems.

In contrast, the squared gradient magnitude $\|\nabla\mathbf{L}(\mathbf{x},\mathbf{y})\|^{2}$, when $\mathbf{L}$ is differentiable, is a more directly observable value for quantifying suboptimality. Moreover, the notion is meaningful for differentiable non-convex-concave minimax games. Interestingly, very few prior works have analyzed convergence rates on the gradient norm for minimax problems, and the optimal convergence rate or corresponding algorithms were hitherto unknown.

In this work, we introduce the extra anchored gradient (EAG) algorithms for smooth convex-concave minimax problems and establish an accelerated $\|\nabla\mathbf{L}(\mathbf{z}^{k})\|^{2}\leq\mathcal{O}(R^{2}/k^{2})$ rate, where $R$ is the Lipschitz constant of $\nabla\mathbf{L}$. The rate improves upon the $\mathcal{O}(R^{2}/k)$ rates of prior algorithms and is the first $\mathcal{O}(R^{2}/k^{2})$ rate in this setup. We then provide a matching $\Omega(R^{2}/k^{2})$ complexity lower bound for gradient-based algorithms and thereby establish optimality of EAG.

Beyond establishing the optimal complexity, our results provide the following observations. First, different suboptimality measures lead to materially different acceleration mechanisms, since reducing the duality gap is done optimally by the extragradient algorithm (Nemirovski, 2004; Nemirovsky, 1992). Also, since our optimal accelerated convergence rate is on the non-ergodic last iterate, neither averaging nor keeping track of the best iterate is necessary for optimally reducing the gradient magnitude in the deterministic setup.

1 Preliminaries and notation

2 Prior work

The first main component of our proposed algorithm is the extragradient (EG) algorithm of Korpelevich (1977). EG and its variants, including the algorithm of Popov (1980), have been studied in the context of saddle point and variational inequality problems and have appeared in the mathematical programming literature (Solodov & Svaiter, 1999; Tseng, 2000; Noor, 2003; Censor et al., 2011; Lyashko et al., 2011; Malitsky & Semenov, 2014; Malitsky, 2015, 2020). More recently in the machine learning literature, similar ideas such as optimism (Chiang et al., 2012; Rakhlin & Sridharan, 2013a), prediction (Yadav et al., 2018), and negative momentum (Gidel et al., 2019; Zhang et al., 2020) have been presented and used in the context of multi-player games (Daskalakis et al., 2011; Rakhlin & Sridharan, 2013b; Syrgkanis et al., 2015; Antonakopoulos et al., 2021) and GANs (Gidel et al., 2018; Mertikopoulos et al., 2019; Liang & Stokes, 2019; Peng et al., 2020).

𝓞​(𝑹/𝒌)𝓞𝑹𝒌\boldsymbol{\mathcal{O}(R/k)} rates on duality gap.

Convergence rates on squared gradient norm.

Using standard arguments (e.g. (Solodov & Svaiter, 1999, Lemma 2.3)), one can show $\underset{i=0,\dots,k}{\min}\,\|\mathbf{G}(\mathbf{z}^{i})\|^{2}\leq\mathcal{O}(R^{2}/k)$ convergence rate of EG, provided that $\mathbf{L}$ is $R$-smooth. Ryu et al. (2019) showed that optimistic descent algorithms also attain $\mathcal{O}(R^{2}/k)$ convergence in terms of the best iterate and proposed simultaneous gradient descent with anchoring, which pulls iterates toward the initial point $\mathbf{z}^{0}$, and established $\mathcal{O}(R^{2}/k^{2-2p})$ convergence rates in terms of squared gradient norm of the last iterate (where $p>\frac{1}{2}$ is an algorithm parameter; see Section A). Notably, anchoring resembles the Halpern iteration (Halpern, 1967; Lieder, 2020), which was used in Diakonikolas (2020) to develop a regularization-based algorithm with near-optimal (optimal up to logarithmic factors) complexity with respect to the gradient norm of the last iterate. Anchoring turns out to be the second main component of the acceleration; combining EG steps with anchoring, we obtain the optimal last-iterate convergence rate of $\mathcal{O}(R^{2}/k^{2})$.

Structured minimax problems.

For structured minimax problems of the form

where $f,g$ are convex and $\mathbf{A}$ is a linear operator, primal-dual splitting algorithms (Chambolle & Pock, 2011; Condat, 2013; Vũ, 2013; Yan, 2018; Ryu & Yin, 2021) and Nesterov’s smoothing technique (Nesterov, 2005a, b) have also been extensively studied (Chen et al., 2014; He & Monteiro, 2016). Notably, when $g$ is of “simple” form, Neterov’s smoothing framework achieves an accelerated rate $\mathcal{O}\left(\frac{\|\mathbf{A}\|}{k}+\frac{L_{f}}{k^{2}}\right)$ on duality gap. Additionally, Chambolle & Pock (2016) have shown that splitting algorithms can achieve $\mathcal{O}(1/k^{2})$ or linear convergence rates under appropriate strong convexity and smoothness assumptions on $f$ and $g$, although they rely on proximal operations. Kolossoski & Monteiro (2017); Hamedani & Aybat (2018); Zhao (2019); Alkousa et al. (2020) generalized these accelerated algorithms to the setting where the coupling term $\langle\mathbf{A}\mathbf{x},\mathbf{y}\rangle$ is replaced by non-bilinear convex-concave function $\Phi(\mathbf{x},\mathbf{y})$.

Complexity lower bounds.

Ouyang & Xu (2021) presented a $\Omega\left(\frac{\|\mathbf{A}\|}{k}+\frac{L_{f}}{k^{2}}\right)$ complexity lower bound on duality gap for gradient-based algorithms solving bilinear minimax problems with proximable $g$, establishing optimality of Nesterov’s smoothing. Zhang et al. (2019) presented lower bounds for strongly-convex-strongly-concave problems. Golowich et al. (2020) proved that with the narrower class of $1$-SCLI algorithms, which includes EG but not EAG, the squared gradient norm of the last iterate cannot be reduced beyond $\mathcal{O}(R^{2}/k)$ in $R$-smooth minimax problems. These approaches are aligned with the information-based complexity analysis, introduced in (Nemirovsky & Yudin, 1983) and thoroughly studied in (Nemirovsky, 1991, 1992) for the special case of linear equations.

Other problem setups.

Nesterov (2009) and Nedić & Ozdaglar (2009) proposed subgradient algorithms for non-smooth minimax problems. Stochastic minimax and variational inequality problems were studied in (Nemirovski et al., 2009; Juditsky et al., 2011; Lan, 2012; Ghadimi & Lan, 2012, 2013; Chen et al., 2014, 2017; Hsieh et al., 2019). Strongly monotone variational inequality problems or strongly-convex-strongly-concave minimax problems were studied in (Tseng, 1995; Nesterov & Scrimali, 2011; Gidel et al., 2018; Mokhtari et al., 2020a; Lin et al., 2020b; Wang & Li, 2020; Zhang et al., 2020; Azizian et al., 2020). Recently, minimax problems with objectives that are either strongly convex or nonconvex in one variable were studied in (Rafique et al., 2018; Thekumparampil et al., 2019; Jin et al., 2019; Nouiehed et al., 2019; Ostrovskii et al., 2020; Lin et al., 2020a, b; Lu et al., 2020; Wang & Li, 2020; Yang et al., 2020; Chen et al., 2021). Minimax optimization of composite objectives with smooth and nonsmooth-but-proximable convex-concave functions were studied in (Tseng, 2000; Csetnek et al., 2019; Malitsky & Tam, 2020; Bùi & Combettes, 2021).

Accelerated algorithms: Extra anchored gradient

We now present two accelerated EAG algorithms that are qualitatively very similar but differ in the choice of step-sizes. The two algorithms present a tradeoff between the simplicity of the step-size and the simplicity of the convergence proof; one algorithm has a varying step-size but a simpler convergence proof, while the other algorithm has a simpler constant step-size but has a more complicated proof.

The proposed extra anchored gradient (EAG) algorithms have the following general form:

The simplest choice of $\{\alpha_{k}\}_{k\geq 0}$ is the constant one. Together with the choice $\beta_{k}=\frac{1}{k+2}$ (which we clarify later), we get the following simpler algorithm.

In the setup of Theorem 1, $\alpha\in\left(0,\frac{1}{8R}\right]$ satisfies (4), and the particular choice $\alpha=\frac{1}{8R}$ yields

While EAG-C is simple in its form, its convergence proof (presented in the appendix) is complicated. Furthermore, the constant $260$ in Corollary 1 seems large and raises the question of whether it could be reduced. These issues, to some extent, are addressed by the following alternative version of EAG.

EAG with varying step-size (EAG-V)

where $\alpha_{0}\in\left(0,\frac{1}{R}\right)$ and

As the recurrence relation (2.1) may seem unfamiliar, we provide the following lemma describing the behavior of the resulting sequence.

If $\alpha_{0}\in\left(0,\frac{3}{4R}\right)$, then the sequence $\{\alpha_{k}\}_{k\geq 0}$ of (2.1) monotonically decreases to a positive limit. In particular, when $\alpha_{0}=\frac{0.618}{R}$, we have $\lim_{k\rightarrow\infty}\alpha_{k}\approx\frac{0.437}{R}$.

We now state the convergence results for EAG-V.

EAG-V with $\alpha_{0}=\frac{0.618}{R}$ satisfies

2 Proof outline

We now outline the convergence analysis for EAG-V, whose proof is simpler than that of EAG-C. The key ingredient of the proof is a Lyapunov analysis with a nonincreasing Lyapunov function, the $V_{k}$ of the following lemma.

Let $\{\beta_{k}\}_{k\geq 0}\subseteq(0,1)$ and $\alpha_{0}\in\left(0,\frac{1}{R}\right)$ be given. Define the sequences $\{A_{k}\}_{k\geq 0},\{B_{k}\}_{k\geq 0}$ and $\{\alpha_{k}\}_{\geq 0}$ by the recurrence relations

for $k\geq 0$, where $B_{0}=1$. Suppose that $\alpha_{k}\in(0,\frac{1}{R})$ holds for all $k\geq 0$. Assume $\mathbf{L}$ is $R$-smooth and convex-concave. Then the sequence $\{V_{k}\}_{k\geq 0}$ defined as

for EAG iterations in (3) is nonincreasing.

In Lemma 2, the choice of $\beta_{k}=\frac{1}{k+2}$ leads to $B_{k}=k+1$, $A_{k}=\frac{\alpha_{k}(k+2)(k+1)}{2}$, and (2.1). Why the Lyapunov function of Lemma 2 leads to the convergence guarantee of Theorem 2 may not be immediately obvious. The following proof provides the analysis.

Let $\beta_{k}=\frac{1}{k+2}$ as specified by the definition of EAG-V. By Lemma 2, the quantity $V_{k}$ defined by (9) is nonincreasing in $k$. Therefore,

where (2.2) follows from the monotonicity inequality $\langle\mathbf{G}(\mathbf{z}^{k}),\mathbf{z}^{k}-\mathbf{z}^{\star}\rangle\geq 0$, (2.2) follows from Young’s inequality, (2.2) follows from plugging in $A_{k}=\frac{\alpha_{k}(k+1)(k+2)}{2}$ and $B_{k}=k+1$, and (2.2) follows from Lemma 1 ($\alpha_{k}\downarrow\alpha_{\infty}$). Reorganize to get

and divide both sides by $\frac{\alpha_{\infty}}{4}(k+1)(k+2)$. ∎

3 Discussion of further generalizations

The algorithms and results of Sections 2.1 and 2.2 remain valid when we replace $\mathbf{G}$ with an $R$-Lipschitz continuous monotone operator; neither the definition of the EAG algorithms nor any part of the proofs of Theorems 1 and 2 utilize properties of saddle functions beyond the monotonicity of their subdifferentials.

For EAG-C, the step-size conditions (4) in Theorem 1 can be relaxed to accommodate larger values of $\alpha$. However, we do not pursue such generalizations to keep the already complicated and arduous analysis of EAG-C manageable. Also, larger step-sizes are more naturally allowed in EAG-V and Theorem 2. Finally, although (4) holds for values of $\alpha$ up to $\frac{0.1265}{R}$, we present a slightly smaller range $\left(0,\frac{1}{8R}\right]$ in Corollary 1 for simplicity.

For EAG-V, the choice $\beta_{k}=\frac{1}{k+2}$ was obtained by roughly, but not fully, optimizing the bound on EAG-V originating from Lemma 2. If one chooses $\beta_{k}=\frac{1}{k+\delta}$ with $\delta>1$, then (6) and (7) become

As the proof of Theorem 2 illustrates, linear growth of $B_{k}$ and quadratic growth of $A_{k}$ leads to $\mathcal{O}(1/k^{2})$ convergence of $\|\mathbf{G}(\mathbf{z}^{k})\|^{2}$. The value $\alpha_{0}=\frac{0.618}{R}$ in Lemma 1 and Corollary 2 was obtained by numerically minimizing the constant $\frac{4}{\alpha_{\infty}^{2}}\left(1+\alpha_{0}\alpha_{\infty}R^{2}\right)$ in Theorem 2 in the case of $\delta=2$. The choice $\delta=2$, however, is not optimal. Indeed, the constant $27$ of Corollary 2 can be reduced to $24.44$ with $(\delta^{\star},\alpha_{0}^{\star})\approx(2.697,0.690/R)$, which was obtained by numerically optimizing over $\delta$ and $\alpha_{0}$. Finally, there is a possibility that a choice of $\beta_{k}$ not in the form of $\beta_{k}=\frac{1}{k+\delta}$ leads to an improved constant.

In the end, we choose to present EAG-C and EAG-V with the simple choice $\beta_{k}=\frac{1}{k+2}$. As we establish in Section 3, the EAG algorithms are optimal up to a constant.

Optimality of EAG via a matching complexity lower bound

Upon seeing an accelerated algorithm, it is natural to ask whether the algorithm is optimal. In this section, we present a $\Omega(R^{2}/k^{2})$ complexity lower bound for the class of deterministic gradient-based algorithms for smooth convex-concave minimax problems. This result establishes that EAG is indeed optimal.

For the class of smooth minimax optimization problems, a deterministic algorithm $\mathcal{A}$ produces iterates $(\mathbf{x}^{k},\mathbf{y}^{k})=\mathbf{z}^{k}$ for $k\geq 1$ given a starting point $(\mathbf{x}^{0},\mathbf{y}^{0})=\mathbf{z}^{0}$ and a saddle function $\mathbf{L}$, and we write $\mathbf{z}^{k}=\mathcal{A}(\mathbf{z}^{0},\dots,\mathbf{z}^{k-1};\mathbf{L})$ for $k\geq 1$. Define $\mathfrak{A}_{\textrm{sim}}$ as the class of algorithms satisfying

and $\mathfrak{A}_{\textrm{sep}}$ as the class of algorithms satisfying

To clarify, algorithms in $\mathfrak{A}_{\textrm{sim}}$ access and utilize the $\mathbf{x}$- and $\mathbf{y}$-subgradients simultaneously. So $\mathfrak{A}_{\textrm{sim}}$ contains simultaneous gradient descent, extragradient, Popov, and EAG (if we also count intermediate sequences $\mathbf{z}^{k+1/2}$ as algorithms’ iterates). On the other hand, algorithms in $\mathfrak{A}_{\textrm{sep}}$ can access and utilize the $\mathbf{x}$- and $\mathbf{y}$-subgradients separately. So $\mathfrak{A}_{\textrm{sim}}\subset\mathfrak{A}_{\textrm{sep}}$, and alternating gradient descent-ascent belongs to $\mathfrak{A}_{\textrm{sep}}$ but not to $\mathfrak{A}_{\textrm{sim}}$.

In this section, we present a complexity lower bound that applies to all algorithms in $\mathfrak{A}_{\textrm{sep}}$, not just the algorithms in $\mathfrak{A}_{\textrm{sim}}$. Although EAG-C and EAG-V are in $\mathfrak{A}_{\textrm{sim}}$, we consider the broader class $\mathfrak{A}_{\textrm{sep}}$ to rule out the possibility that separately updating the $\mathbf{x}$- and $\mathbf{y}$-variables provides an improvement beyond a constant factor.

We say $\mathbf{L}(\mathbf{x},\mathbf{y})$ is biaffine if it is an affine function of $\mathbf{x}$ for any fixed $\mathbf{y}$ and an affine function of $\mathbf{y}$ for any fixed $\mathbf{x}$. Biaffine functions are, of course, convex-concave. We first establish a complexity lower bound on minimiax optimization problems with biaffine loss functions.

holds for any algorithm in $\mathfrak{A}_{\textrm{sep}}$, where $\left\lfloor\cdot\right\rfloor$ is the floor function and $\mathbf{z}^{\star}$ is the saddle point of $\mathbf{L}$ closest to $\mathbf{z}^{0}$. Moreover, this lower bound is optimal in the sense that it cannot be improved with biaffine functions.

Since smooth biaffine functions are special cases of smooth convex-concave functions, Theorem 3 implies the optimality of EAG applied to smooth convex-concave mimimax optimization problems.

For $R$-smooth convex-concave minimax problems, an algorithm in $\mathfrak{A}_{\textrm{sep}}$ cannot attain a worst-case convergence rate better than

with respect to $\|\nabla\mathbf{L}(\mathbf{z}^{k})\|^{2}$. Since EAG-C and EAG-V have rates $\mathcal{O}(R^{2}\|\mathbf{z}^{0}-\mathbf{z}^{\star}\|^{2}/k^{2})$, they are optimal, up to a constant factor, in $\mathfrak{A}_{\textrm{sep}}$.

are characterized by $\mathbf{A}\mathbf{x}-\mathbf{b}=0$ and $\mathbf{A}^{\intercal}(\mathbf{y}-\mathbf{c})=0$.

Through translation, we may assume without loss of generality that $\mathbf{x}^{0}=0,\mathbf{y}^{0}=0$. In this case, (11) becomes

for $k\geq 2$. (We detail these arguments in the appendix.) Furthermore let $\mathbf{A}=\mathbf{A}^{\intercal}$ and $\mathbf{b}=\mathbf{A}^{\intercal}\mathbf{c}=\mathbf{A}\mathbf{c}$. Then the characterization of $\mathfrak{A}_{\textrm{sep}}$ further simplifies to

Note that $\mathcal{K}_{k-1}(\mathbf{A};\mathbf{b})$ is the order-$(k-1)$ Krylov subspace.

Consider the following lemma. Its proof, deferred to the appendix, combines arguments from Nemirovsky (1991, 1992).

for any $\mathbf{x}\in\mathcal{K}_{k-1}(\mathbf{A};\mathbf{b})$, where $\mathbf{x}^{\star}$ is the minimum norm solution to the equation $\mathbf{A}\mathbf{x}=\mathbf{b}$.

Take $\mathbf{A}$ and $\mathbf{b}$ as in Lemma 3 and $\mathbf{c}=\mathbf{x}^{\star}$. Then $\mathbf{z}^{\star}=(\mathbf{x}^{\star},\mathbf{x}^{\star})$ is the saddle point of $\mathbf{L}(\mathbf{x},\mathbf{y})=\langle\mathbf{A}\mathbf{x}-\mathbf{b},\mathbf{y}-\mathbf{c}\rangle$ with minimum norm. Finally,

for any $\mathbf{x}^{k},\mathbf{y}^{k}\in\mathcal{K}_{k-1}(\mathbf{A};\mathbf{b})$. This completes the construction of the biaffine $\mathbf{L}$ of Theorem 3.

2 Optimal complexity lower bound

We now formalize the notion of complexity lower bounds. This formulation will allow us to precisely state and prove the second statement of Theorem 3 regarding the optimality of the lower bound.

Let $\mathcal{F}$ be a function class, $\mathcal{P}_{\mathcal{F}}=\{\mathcal{P}_{f}\}_{f\in\mathcal{F}}$ a class of optimization problems (with some common form), and $\mathcal{E}(\cdot;\mathcal{P}_{f})$ a suboptimality measure for the problem $\mathcal{P}_{f}$. Define the worst-case complexity of an algorithm $\mathcal{A}$ for $\mathcal{P}_{\mathcal{F}}$ at the $k$-th iteration given the initial condition $\|\mathbf{z}^{0}-\mathbf{z}^{\star}\|\leq D$, as

where $\mathbf{z}^{j}=\mathcal{A}(\mathbf{z}^{0},\dots,\mathbf{z}^{j-1};f)$ for $j=1,\dots,k$ and $B(\mathbf{z};D)$ denotes the closed ball of radius $D$ centered at $\mathbf{z}$. The optimal complexity lower bound with respect to an algorithm class $\mathfrak{A}$ is

A complexity lower bound is a lower bound on the optimal complexity lower bound.

As an aside, the argument of Corollary 3 can be expressed as: for any $\mathcal{A}\in\mathfrak{A}_{\textrm{sep}}$, we have

The first inequality follows from $\mathcal{A}\in\mathfrak{A}_{\textrm{sep}}$, the second from $\mathcal{L}_{R}^{\textrm{biaff}}\subset\mathcal{L}_{R}$, and the third from Theorem 3.

Using above notations, our goal is to prove that for $n\geq k+2$,

We establish this claim with the chain of inequalities:

Inequality (17) is what we established in Section 3.1. Inequality (18) follows from $\mathfrak{A}_{\textrm{sim}}\subset\mathfrak{A}_{\textrm{sep}}$ and the fact that the infimum over a larger class is smaller. Roughly speaking, the quantities in lines (19) and (20) are the complexity lower bounds for solving linear equations using only matrix-vector products, which were studied thoroughly in (Nemirovsky, 1991, 1992). We will show inequalities (19), (20), and (21) by establishing the connection of Nemirovsky’s work with our setup of biaffine saddle problems. Once this is done, equality holds throughout and (16) is proved.

We first provide the definitions. Let $\mathcal{P}^{2n}_{R,D}$ be the collection of linear equations with $2n\times 2n$ matrices $\mathbf{B}$ satisfying $\|\mathbf{B}\|\leq R$ and $\mathbf{v}=\mathbf{B}\mathbf{z}^{\star}$ for some $\mathbf{z}^{\star}\in B(0;D)$. Let $\mathcal{P}^{2n,\textrm{skew}}_{R,D}\subset\mathcal{P}^{2n}_{R,D}$ be the subclass of equations with skew-symmetric $\mathbf{B}$. Let $\mathfrak{A}_{\textrm{lin}}$ be the class of iterative algorithms solving linear equations $\mathbf{B}\mathbf{z}=\mathbf{v}$ using only matrix multiplication by $\mathbf{B}$ and $\mathbf{B}^{\intercal}$ in the sense that

where $\mathbf{v}^{0}=0$, $\mathbf{v}^{1}=\mathbf{v}$, and for $k\geq 2$,

The optimal complexity lower bound for a class of linear equation instances is defined as

Define $\mathcal{C}\left(\mathfrak{A}_{\textrm{lin}};\mathcal{P}^{2n,\textrm{skew}}_{R,D},k\right)$ analogously.

Now we relate the optimal complexity lower bounds for biaffine minimax problems to those for linear equations. For $\mathbf{L}(\mathbf{x},\mathbf{y})=\mathbf{b}^{\intercal}\mathbf{x}+\mathbf{x}^{\intercal}\mathbf{A}\mathbf{y}-\mathbf{c}^{\intercal}\mathbf{y}$, we have

For both algorithm classes $\mathfrak{A}_{\textrm{sim}}$ and $\mathfrak{A}_{\textrm{lin}}$, we may assume without loss of generality that $\mathbf{z}^{0}=0$ through translation. Then, the span condition (10) for $\mathfrak{A}_{\textrm{sim}}$ becomes

Since the supremum over a larger class of problems is larger, inequality (19) holds. Similarly, inequality (20) follows from $\mathcal{P}^{2n,\textrm{skew}}_{R,D}\subset\mathcal{P}^{2n}_{R,D}$.

Finally, (21) follows from the following lemma, using arguments based on Chebyshev-type matrix polynomials from Nemirovsky (1992). Its proof is deferred to the appendix.

3 Broader algorithm classes via resisting oracles

In (10) and (11), we assumed the subgradient queries are made within the span of the gradients at the previous iterates. This requirement (the linear span assumption) can be removed, i.e., a similar analysis can be done on general deterministic black-box gradient-based algorithms (formally defined in the appendix, Section C.5), using the resisting oracle technique (Nemirovsky & Yudin, 1983) at the cost of slightly enlarging the required problem dimension. We informally state the generalized result below and provide details in the appendix.

Although we do not formally pursue this, the requirement that the algorithm is not randomized can also be removed using the techniques of Woodworth & Srebro (2016), which exploit near-orthogonality of random vectors in high dimensions.

4 Discussion

We established that one cannot improve the lower bound of Theorem 3 using biaffine functions, arguably the simplest family of convex-concave functions. Furthermore, this optimality statement holds for both algorithm classes $\mathfrak{A}_{\textrm{sep}}$ and $\mathfrak{A}_{\textrm{sim}}$ as established through the chain of inequalities in Section 3.2. However, as demonstrated by Drori (2017), who introduced a non-quadratic lower bound for smooth convex minimization that improves upon the classical quadratic lower bounds of Nemirovsky (1992) and Nesterov (2013), a non-biaffine construction may improve the constant. In our setup, there is a factor-near-$100$ difference between the upper and lower bounds. (Note that each EAG iteration requires $2$ evaluations of the saddle subdifferential oracle.) We suspect that both the algorithm and the lower bound can be improved upon, but we leave this to future work.

Golowich et al. (2020) establishes that for the class of 1-SCLI algorithms (S is for stationary), a subclass of $\mathfrak{A}_{\textrm{sim}}$ for biaffine objectives, one cannot achieve a rate faster than $\|\nabla\mathbf{L}(\mathbf{z}^{k})\|^{2}\leq\mathcal{O}(1/k)$. This lower bound applies to EG but not EAG; EAG is not 1-SCLI, as its anchoring coefficients $\frac{1}{k+2}$ vary over iterations, and its convergence rate breaks the 1-SCLI lower bound. On the other hand, we can view EAG as a non-stationary CLI algorithm (Arjevani & Shamir, 2016, Definition 2). We further discuss these connections in the appendix, Section E.

Experiments

We now present experiments illustrating the accelerated rate of EAG. We compare EAG-C and EAG-V against the prior algorithms with convergence guarantees: EG, Popov’s algorithm (or optimistic descent) and simultaneous gradient descent with anchoring (SimGD-A). The precise forms of the algorithms are restated in the appendix.

Figure 1(a) presents experiments on our first example, constructed as follows. For $\epsilon>0$, define

Next, for $0<\epsilon\ll\delta\ll 1$, define

Figure 1(b) presents experiments on our second example

Figure 2(a) illustrates the algorithms applied to (24). For $|x|,|y|\gg\epsilon$,

so the algorithms roughly behave as if the objective is the bilinear function $\delta xy$. When $\delta$ is sufficiently small, trajectories of the algorithms closely resemble the corresponding continuous-time flows with $\mathbf{L}(x,y)=xy$.

Csetnek et al. (2019) demonstrated that Popov’s algorithm can be viewed as discretization of the Moreau–Yosida regularized flow $\dot{\mathbf{z}}(t)=-\frac{\mathbf{G}-(\textrm{Id}+\lambda\mathbf{G})^{-1}}{\lambda}\left(\mathbf{z}(t)\right)$ for some $\lambda>0$, and a similar analysis can be performed with EG. This connection explains why EG’s trajectory in Figure 2(a) and the regularized flow depicted in Figure 2(b) are similar.

On the other hand, EAG and SimGD-A can be viewed as a discretization of the anchored flow ODE

The anchored flow depicted in Figure 2(b) approaches the solution much more quickly due to the anchoring term dampening the cycling behavior. The trajectories of EAG and SimGD-A iterates in Figure 2(a) are very similar to the anchored flow. However, SimGD-A requires diminishing step-sizes $\frac{1-p}{(k+1)^{p}}$ (both theoretically and experimentally) and therefore progresses much slower.

Conclusion

This work presents the extra anchored gradient (EAG) algorithms, which exhibit accelerated $\mathcal{O}(1/k^{2})$ rates on the squared gradient magnitude for smooth convex-concave minimax problems. The acceleration combines the extragradient and anchoring mechanisms, which separately achieve $\mathcal{O}(1/k)$ or slower rates. We complement the $\mathcal{O}(1/k^{2})$ rate with a matching $\Omega(1/k^{2})$ complexity lower bound, thereby establishing optimality of EAG.

At a superficial level, the acceleration mechanism of EAG seems to be distinct from that of Nesterov; anchoring dampens oscillations, but momentum provides the opposite effect of dampening. However, are the two accelerations phenomena entirely unrelated? Finding a common structure, a connection, between the two acceleration phenomena would be an interesting direction of future work.

Acknowledgements

TY and EKR were supported by the National Research Foundation of Korea (NRF) Grant funded by the Korean Government (MSIP) [No. 2020R1F1A1A01072877], the National Research Foundation of Korea (NRF) Grant funded by the Korean Government (MSIP) [No. 2017R1A5A1015626], by the New Faculty Startup Fund from Seoul National University, and by the AI Institute of Seoul National University (AIIS) through its AI Frontier Research Grant (No. 0670-20200015) in 2020. We thank Jaewook Suh and Jongmin Lee for reviewing the manuscript and providing valuable feedback. We thank Jelena Diakonikolas for the discussion on the prior work on parameter-free near-optimal methods for the smooth minimax setup. Finally, we thank the anonymous referees for bringing to our attention the recent complexity lower bound on the class of $1$-SCLI algorithms by Golowich et al. (2020).

References

Appendix A Algorithm specifications

For the sake of clarity, we precisely specify all the algorithms discussed in this work.

Simultaneous gradient descent for smooth minimax optimization is defined as

The notation becomes more concise with the joint variable notation $\mathbf{z}^{k}=(\mathbf{x}^{k},\mathbf{y}^{k})$ and the saddle operator (2), where the sign change in $\mathbf{y}$-gradient is already included:

Alternating gradient descent-ascent is defined as

Note that we update the $\mathbf{x}$ variable first and then use it to update the $\mathbf{y}$-iterate.

The extragradient (EG) algorithm is defined as

Popov’s algorithm, or optimistic descent, is defined as

Simultaneous gradient descent with anchoring (SimGD-A) (Ryu et al., 2019) is defined as

where $p\in(1/2,1)$ and $\gamma>0$. It has been proved in Ryu et al. (2019) that SimGD-A converges at $\mathcal{O}(1/k^{2-2p})$ rate. In this paper, we always used $\gamma=1$ and $p=\frac{1}{2}+10^{-2}$.

Appendix B Omitted proofs of Section 2

The following identities follow directly from the definition of EAG iterates:

Recall that $\mathbf{G}$ is a monotone operator, so that

where (B.1) follows from (26) and (28), and (29) results from cancellation and collection of terms using (7). Next, we have

from $R$-Lipschitzness of $\mathbf{G}$ and (27). Now multiplying the factor $\frac{A_{k}}{\alpha_{k}^{2}R^{2}}$ to (30) and subtracting from (29) gives

Observe that the $\left\langle\mathbf{G}(\mathbf{z}^{k}),\mathbf{G}(\mathbf{z}^{k+1/2})\right\rangle$ term vanishes because of (6), and that

Plugging these identities into (B.1) and simplifying, we get

where the last inequality is an application of Young’s inequality.

B.2 Proof of Lemma 1

We may assume $R=1$ without loss of generality because we can recover the general case by replacing $\alpha_{k}$ with $\alpha_{k}R$. Rewrite (2.1) as

Suppose that we have already established $0<\alpha_{N}<\rho$ for some $N\geq 0$ and $\rho\in(0,1)$, where $\rho$ satisfies

Note that (33) holds true for all $N\geq 0$ if $\rho<\frac{3}{4}$. Now we will show that given (33),

so that $\alpha_{k}\downarrow\alpha$ for some $\alpha\geq(1-\gamma)\alpha_{N}$. It suffices to prove that $(1-\gamma)\alpha_{N}<\alpha_{N+k}<\rho$ for all $k\geq 0$, because it is clear from (32) that $\{\alpha_{k}\}_{k\geq 0}$ is decreasing.

We use induction on $k$ to prove that $\alpha_{N+k}\in((1-\gamma)\alpha_{N},\rho)$. The case $k=0$ is trivial. Now suppose that $(1-\gamma)\alpha_{N}<\alpha_{N+j}<\rho$ holds true for all $j=0,\dots,k$. Then by (32), for each $0\leq j\leq k$ we have

Summing up the inequalities for $j=0,\dots,k$, we obtain

which gives $(1-\gamma)\alpha_{N}<\alpha_{N+k+1}<\alpha_{N}<\rho$, completing the induction.

In particular, when $\alpha_{0}=0.618$, direct calculation gives $0.437>\alpha_{N}>0.4366$ when $N=1000$. With $\rho=0.437$ and $N=1000$, we have $\gamma=\frac{1}{2}\left(\frac{1}{N+1}+\frac{1}{N+2}\right)\frac{\rho^{2}}{1-\rho^{2}}<2.5\times 10^{-4}$, which gives $\alpha\geq(1-\gamma)\alpha_{N}\approx 0.4365$.

B.3 Proof of Theorem 1

As in the proof of Theorem 2, assume without loss of generality that $R=1$. The strategy of the proof is basically the same as in Theorem 2; we construct a nonincreasing Lyapunov function by combining the same set of inequalities, but with different (more intricate) coefficients. For $k\geq 0$, let

As in Lemma 2, we will use $B_{k}=\frac{1}{1-\beta_{k}}=k+1$, and $a_{k}\geq 0$ will be specified later. Because we have the fixed step-size $\alpha$, the identities (26), (27), and (28) become

Now, subtracting the same inequalities from monotonicity and Lipschitzness from $V_{k}-V_{k+1}$ as in Lemma 2, each with coefficients $(k+1)(k+2)$ and $\tau_{k}\geq 0$ (to be specified later), we obtain

where we define $\mathbf{M}_{k}:=\begin{bmatrix}\mathbf{G}(\mathbf{z}^{k})&\mathbf{G}(\mathbf{z}^{k+1/2})&\mathbf{G}(\mathbf{z}^{k+1})\end{bmatrix}$ and

Subdivide the interval $I_{k}$ into two parts:

We divide cases: $A_{k}\in I_{k}^{-}$ and $A_{k}\in I_{k}^{+}$. However, the latter case is in fact not needed unless we wish to extend the proof for $\alpha$ beyond $\frac{0.1265}{R}$. If that is not the case, we recommend the readers to refer to Case 1 only. Nevertheless, we exhibit analysis of both cases because Case 2 might provide useful data for enlarging or even completely determining the range of convergent step-sizes for EAG-C.

Case 1. Suppose that $A_{k}\in I_{k}^{-}$. In this case, we choose

The denominator and numerator of (40) are both positive because $u_{k}>A_{k}>\frac{\alpha(k+1)(k+1-\alpha(k+2))}{2(1-\alpha)}$ (see (38)). Thus, $\tau_{k}>0$. Next, define $A_{k+1}$ as

The expressions seem ridiculously complicated, but there are a number of repeating terms. Let

Because $A_{k}\leq\frac{\alpha(k+1)(k+2)}{2}<u_{k}$ (see (35)), we have $E_{1}>0,E_{2}\geq 0$. (Note that $E_{2}=0$ only in the boundary case $A_{k}=\sup I_{k}^{-}$.) Next, put

which is a factor that appears within the definition of $\tau_{k}$ (40); we have already seen that $E_{3}>0$. Further, let

It is obvious that $E_{5},E_{7}>0$, and $E_{6}>0$ follows directly from (37). To see that $E_{4}>0$, observe that $k+1-\alpha(k+2)=(k+2)\left(\frac{k+1}{k+2}-\alpha\right)\geq(k+2)\left(\frac{1}{2}-\alpha\right)\geq 0$, provided that $\alpha\leq\frac{1}{2}$. This implies

This immediately shows that the diagonal entries $s_{ii}$ are nonnegative for $i=1,2,3$. By brute-force calculation, it is not difficult to verify the identity

Using this, we see that $\mathbf{v}:=\begin{bmatrix}\frac{\alpha(k+2)E_{7}}{2E_{5}}&\frac{E_{4}}{2(1-\alpha)E_{5}}&1\end{bmatrix}^{\intercal}$ satisfies $\mathbf{S}_{k}\mathbf{v}=0$, and this implies $\det\mathbf{S}_{k}=0$. The cofactor-expansion of $\det\mathbf{S}_{k}$ along the first row gives

when $s_{11}>0$, and via continuity argument we can argue that $\begin{vmatrix}s_{22}&s_{23}\\ s_{23}&s_{33}\end{vmatrix}\geq 0$ even in the boundary case $s_{11}=0$. Similarly one can show that $\begin{vmatrix}s_{11}&s_{12}\\ s_{12}&s_{22}\end{vmatrix}\geq 0$. Therefore, we have shown that all diagonal submatrices of $\mathbf{S}_{k}$ (including the trivial case $\begin{vmatrix}s_{11}&0\\ 0&s_{33}\end{vmatrix}=s_{11}s_{33}\geq 0$) have nonnegative determinants, that is, $\mathbf{S}_{k}\succeq\mathbf{O}$.

Finally, (41) shows that $A_{k+1}$ is increasing with respect to $A_{k}$. We see that

and the last expression is nonnegative because of the assumption (4), which we restate here for the case $R=1$ for convenience: $1-3\alpha-\alpha^{2}-\alpha^{3}\geq 0$ and $1-8\alpha+\alpha^{2}-2\alpha^{3}\geq 0$. This proves that $A_{k+1}\in I_{k+1}^{-}\subset I_{k+1}$, as desired.

Case 2. Suppose that $A_{k}\in I_{k}^{+}$. The proof would be similar to Case 1, but choices of $\tau_{k}$ and $A_{k+1}$ are different. We let

and so on. (Note that $2A_{k}-\alpha(k+1)(k+2)\geq 0$ because now we are assuming that $A_{k}\in I_{k}^{+}$.) We omit further details of calculations, but with the above choices of $\tau_{k}$ and $A_{k+1}$ it can be shown that $\det\mathbf{S}_{k}=0$ and $s_{11},s_{33}\geq 0$, using (36) through (39). As in Case 1, this implies $\mathbf{S}_{k}\succeq\mathbf{O}$.

and the last term is positive for any $\alpha\in(0,1)$, i.e., $A_{k+1}<u_{k+1}$. This completes Case 2.

where the first inequality follows from Lipschitzness of $\mathbf{G}$ (recall that we are assuming that $R=1$). Also by (35) and (36),

where (B.3) follows from (50) and the monotonicity inequality $\langle\mathbf{z}^{k}-\mathbf{z}^{\star},\mathbf{G}(\mathbf{z}^{k})\rangle\geq 0$, and (B.3) follows from Young’s inequality. Rearranging terms, we conclude that

where $C=\frac{4(1+\alpha+\alpha^{2})}{\alpha^{2}(1+\alpha)}$.

Proof of Lemma 5. Direct calculation gives

because $k+1-\alpha(k+2)=(k+2)(\frac{k+1}{k+2}-\alpha)\geq(k+2)(\frac{1}{2}-\alpha)\geq 0$, which shows (36). Similarly, we observe that

and each line corresponds to an inequality within (37), (38) and (39).

Appendix C Omitted proofs of Section 3

In this section, we provide a self-contained discussion on the complexity lower bound results for linear operator equations from Nemirovsky (1991, 1992).

The proof of Theorem 3 was essentially completed in the main body of the paper, except the argument regarding translation, (13), and the proof of Lemma 3.

Let $\mathbf{z}^{0}=(\mathbf{x}^{0},\mathbf{y}^{0})$ and $\mathbf{L}(\mathbf{x},\mathbf{y})=\langle\mathbf{A}\mathbf{x}-\mathbf{b},\mathbf{y}-\mathbf{c}\rangle$ be given, and assume that $\|\mathbf{z}_{\mathbf{L}}^{\star}(\mathbf{z}^{0})-\mathbf{z}^{0}\|\leq D$. Let $\mathbf{b}_{0}=\mathbf{b}-\mathbf{A}\mathbf{x}^{0}$ and $\mathbf{c}_{0}=\mathbf{c}-\mathbf{y}^{0}$. Then

As one can see, we have $\mathbf{x}^{k}-\mathbf{x}^{0}\in\mathcal{X}_{k}(\mathbf{A};\mathbf{b}_{0},\mathbf{c}_{0})$ and $\mathbf{y}^{k}-\mathbf{y}^{0}\in\mathcal{Y}_{k}(\mathbf{A};\mathbf{b}_{0},\mathbf{c}_{0})$, where we inductively define

Then it is not difficult to see that for $k\geq 2$,

Now consider $\mathbf{L}_{0}(\mathbf{x},\mathbf{y}):=\langle\mathbf{A}\mathbf{x}-\mathbf{b}_{0},\mathbf{y}-\mathbf{c}_{0}\rangle=\left\langle\mathbf{A}(\mathbf{x}+\mathbf{x}^{0})-\mathbf{b},\mathbf{y}+\mathbf{y}^{0}-\mathbf{c}\right\rangle$. Because $\mathbf{z}^{\star}_{\mathbf{L}_{0}}$ is a saddle point of $\mathbf{L}_{0}$ if and only if $\mathbf{z}^{\star}_{\mathbf{L}_{0}}+\mathbf{z}^{0}$ is a saddle point of $\mathbf{L}$, we have $\mathbf{z}^{\star}_{\mathbf{L}_{0}}(0)=\mathbf{z}^{\star}_{\mathbf{L}}(\mathbf{z}^{0})-\mathbf{z}^{0}$, and thus $\|\mathbf{z}^{\star}_{\mathbf{L}_{0}}(0)\|\leq D$. Therefore, if we let

C.2 Complexity of solving linear operator equations and minimax polynomials

We define the problem class by $\|\mathbf{A}\|\leq R$, which is equivalent to $\lambda_{j}\in[-R,R]$ for all $j=1,\dots,n$. Therefore, we consider a method corresponding to a polynomial $q(t)$ such that $p(t)=1-tq(t)$ minimizes

More precisely, if $p_{k}^{\star}(t)=1-tq_{k}^{\star}(t)$ minimizes the last quantity among all $p(t)$ such that $\deg p\leq k$ and $p(0)=1$, and if we put $\mathbf{x}^{k}=q_{k}^{\star}(\mathbf{A})\mathbf{b}$, then (52) implies

for all $\mathbf{A}$ whose spectrum belongs to $[-R,R]$ and $\mathbf{b}=\mathbf{A}\mathbf{x}^{\star}$ with $\|\mathbf{x}^{\star}\|\leq D$. As $p_{k}^{\star}$ solves (53), it is called a minimax polynomial.

In order to establish Lemma 3, we present a two-fold analysis in the following. First, we compute the quantity (53) by explicitly naming $p_{k}^{\star}$ for each $k\geq 1$. (This was given by Nemirovsky (1992), but without a proof.) Then, following the exposition from (Nemirovsky, 1991), we show that there exists an instance of $(\mathbf{A},\mathbf{b})$ such that

holds for any polynomial $q$ of degree $\leq k-1$.

C.3 Proof of Lemma 3

The solutions to (53) are characterized using the Chebyshev polynomials of first kind, defined by

or equivalently by $T_{N}(t)=\cos(N\arccos t)$. If $N=2d$ for some nonnegative integer $d$, then $T_{N}$ is an even polynomial satisfying $T_{N}(0)=\cos(d\pi)=(-1)^{d}$. On the other hand, if $N=2d+1$, then $T_{N}$ is an odd polynomial of the form

which can be shown via induction using the recurrence relation $T_{N+1}(t)=2tT_{N}(t)-T_{N-1}(t)$, which follows from the trigonometric identity

Based on arguments from (Nemirovsky, 1992; Mason & Handscomb, 2002), we will show that given $k\geq 1$ and $m:=\lfloor\frac{k}{2}\rfloor$,

The Chebyshev polynomials satisfy the equioscillation property which makes them so special: the extrema of $T_{N}$ within $$occur at$t_{j}=\cos\frac{(N-j)\pi}{N}$for$j=0,\dots,N$, and the signs of the extremal values alternate. Indeed, we have$|T_{N}(t)=\cos(N\arccos t)|\leq 1$for all$t\in$, and for each$j=0,\dots,N$,

Also, we have $T_{N}(t_{j})=-T_{N}(t_{j-1})$ for each $j=1,\dots,n$.

Given $k\geq 1$, we denote by $\mathcal{P}_{k}$ the collection of all polynomials $p$ of degree $\leq k$ with $p(0)=1$. Recall that we are to minimize

Observe that $p_{k}^{\star}\in\mathcal{P}_{k}$ due to (54). Next, note that $\lambda p_{k}^{\star}(\lambda)=\frac{(-1)^{m}R}{2m+1}T_{2m+1}(\frac{\lambda}{R})$ has extrema of alternating signs and same magnitude within $[-R,R]$, which occur precisely at $\lambda_{j}:=R\cos\frac{(2m+1-j)\pi}{2m+1}$, where $j=0,\dots,2m+1$. Suppose that $p_{k}^{\star}$ is not a minimizer of $M(p,R)$ over $\mathcal{P}_{k}$, so that there exists $p\in\mathcal{P}_{k}$ such that

As $p$ and $p_{k}^{\star}$ are both polynomials of degree $\leq 2m$ and constant terms 1, we can write

for some polynomial $q$ of degree $\leq 2m-1$. But then $|p(\lambda_{j})|=|p_{k}^{\star}(\lambda_{j})-\lambda_{j}q(\lambda_{j})|<|p_{k}^{\star}(\lambda_{j})|$, which implies that $p_{k}^{\star}(\lambda_{j})$ and $\lambda_{j}q(\lambda_{j})$ have same signs for $j=0,\dots,2m+1$. Now, because $p_{k}^{\star}(\lambda_{j})$ have alternating signs and

we see that the signs of $q(\lambda_{j})$ alternate over $j=0,\dots,m$ and over $j=m+1,\dots,2m+1$, respectively. Therefore, $q$ must have at least one zero in each open interval $(\lambda_{j},\lambda_{j+1})$ for $j=0,\dots,m-1,m+1,\dots,2m$. This implies that $q(t)\equiv 0$ since $\deg q\leq 2m-1$, while $q$ has at least $2m$ zeros. Therefore, we arrive at $p_{k}^{\star}=p$, which is a contradiction.

Furthermore, the above arguments show that the minimization of (55) over $p\in\mathcal{P}_{k}$ is in fact the same as the minimization of

and the final problem from the line (60) is equivalent to

so that $\|\mathbf{x}^{\star}\|=D$. For any given $\mathbf{x}=q(\mathbf{A})\mathbf{b}$ with $\deg q\leq k-1$, we use (52) to rewrite $\|\mathbf{A}\mathbf{x}-\mathbf{b}\|^{2}$ as

where $p(t)=1-tq(t)\in\mathcal{P}_{k}$. But since $(p_{k}^{\star},\boldsymbol{\mu}^{\star})$ is the primal-dual solution pair to the problems (63) and (67), $p_{k}^{\star}$ minimizes $\sum_{j=0}^{2m+1}\mu_{j}^{\star}\lambda_{j}^{2}p(\lambda_{j})^{2}$ within $\mathcal{P}_{k}$. Therefore,

C.4 Proof of Lemma 4

Let $k\geq 0$ be a given (fixed) integer. Consider the polynomial $p_{k}^{\star}$ we defined in the previous section. It is an even polynomial of degree $2\lfloor\frac{k}{2}\rfloor$, and thus $p_{k}^{\star}\left(\sqrt{t}\right)$ is a polynomial in $t$ of degree $\lfloor\frac{k}{2}\rfloor$, whose constant term is $p_{k}^{\star}(0)=1$. Therefore, we can write $p_{k}^{\star}\left(\sqrt{t}\right)=1-tq_{k}(t)$ for some polynomial $q_{k}$. We will show that

We proceed via arguments similar to derivations in C.2. First, observe that

where $|\mathbf{B}|$ is the matrix square root of the positive semidefinite matrix $\mathbf{B}^{\intercal}\mathbf{B}$. Rewriting (68) in terms of $|\mathbf{B}|$, we obtain

Plugging the last equation into (69) gives

Finally, because $|\mathbf{B}|$ is a symmetric matrix whose eigenvalues are within $[0,R]$, we can apply (52) with $|\mathbf{B}|,\mathbf{z}^{\star}$ in places of $\mathbf{A},\mathbf{x}^{\star}$, and use (58) to conclude that

C.5 Proof of Theorem 4

We first describe the general class $\mathfrak{A}$ of algorithms without the linear span assumption. An algorithm $\mathcal{A}$ within $\mathfrak{A}$ is a sequence of deterministic functions $\mathcal{A}_{1},\mathcal{A}_{2},\dots$, each of which having the form

The sequence $\{\mathbf{z}^{i}\}_{i\geq 0}$ are the inquiry points, and $\{\overline{\mathbf{z}}^{i}\}_{i\geq 0}$ are the approximate solutions produced by $\mathcal{A}$. When $k\geq 1$ is the predefined maximum number of iterations, then we assume $\overline{\mathbf{z}}^{k}=\mathbf{z}^{k}$ without loss of generality. Similar definitions for deterministic algorithms have been considered in (Nemirovsky, 1991; Ouyang & Xu, 2021).

Then $\left(\nabla_{\mathbf{x}}\mathbf{L}_{0}(\mathbf{x},\mathbf{y}),\nabla_{\mathbf{y}}\mathbf{L}_{0}(\mathbf{x},\mathbf{y})\right)=\left(\mathbf{A}(\mathbf{y}-\mathbf{y}^{0})-\mathbf{b},\mathbf{A}(\mathbf{x}-\mathbf{x}^{0})-\mathbf{b}\right)$, and $\mathbf{z}^{0}+\left(\mathbf{x}^{\textrm{min}},\mathbf{x}^{\textrm{min}}\right)$ is a saddle point of $\mathbf{L}_{0}$.

We follow the oracle-resisting proof strategy of Nemirovsky (1991), described as follows. For each $i=1,\dots,k$, we inductively define a rotated biaffine function

for $j=0,\dots,i$, where $\mathcal{N}_{i}$ is a subspace of $\ker(\mathbf{A})$ such that $\dim(\mathcal{N}_{i})\leq 2i$. Note that (70) implies that the algorithm iterates $(\mathbf{z}^{j},\overline{\mathbf{z}}^{j})$ for $j=1,\dots,i$ do not change when $\mathbf{L}_{i-1}$ is replaced by $\mathbf{L}_{i}$. Hence, this process sequentially adjusts the objective function $\mathbf{L}$ upon observing an iterate $\mathbf{z}^{i}$ to resist the algorithm from optimizing it efficiently. Indeed, if (71) holds with $i=j=k$, then

for some polynomials $q_{\mathbf{x}},q_{\mathbf{y}}$ of degree $\leq k-1$ and $\mathbf{v}_{\mathbf{x}}^{k},\mathbf{v}_{\mathbf{y}}^{k}\in\mathcal{N}_{i}\subseteq\ker(\mathbf{A})$. Thus

Then the theorem statement follows from the fact that $\mathbf{z}^{\star}=\mathbf{z}^{0}+(\mathbf{U}_{k}\mathbf{x}^{\textrm{min}},\mathbf{U}_{k}\mathbf{x}^{\textrm{min}})$ is a saddle point of $\mathbf{L}_{k}$.

It remains to provide an inductive scheme for choosing $\mathbf{U}_{i}$. We set $\mathbf{U}_{0}=\mathbf{I}$ (so that $\mathbf{A}_{0}=\mathbf{A}$), $\mathcal{N}_{0}=\{0\}$, and define $\mathcal{K}_{-1}(\mathbf{A};\mathbf{b})=\{0\}$ for convenience. Let $1\leq i\leq k$, and suppose that we already have an orthogonal matrix $\mathbf{U}_{i-1}$ and $\mathcal{N}_{i-1}\subseteq\ker(\mathbf{A})$ for which $\mathbf{U}_{i-1}\mathbf{b}=\mathbf{b}$, $\dim(\mathcal{N}_{i-1})\leq 2i-2$, and (71) holds with $i-1$ (which is vacuously true when $i=1$). Let

We want $\mathbf{U}_{i}$ (to be defined) to satisfy $\mathbf{s}_{\mathbf{x}}^{i},\mathbf{s}_{\mathbf{y}}^{i}\in\mathbf{U}_{i}\ker(\mathbf{A})$ while $\mathcal{K}_{i-1}(\mathbf{A}_{i-1};\mathbf{b})=\mathcal{K}_{i-1}(\mathbf{A}_{i};\mathbf{b})$. The latter condition is satisfied if $\mathbf{U}_{i}=\mathbf{Q}_{i}\mathbf{U}_{i-1}$ for some orthogonal matrix $\mathbf{Q}_{i}$ which preserves every element within

because then it follows that $\mathbf{U}_{i}\mathbf{b}=\mathbf{Q}_{i}\mathbf{U}_{i-1}\mathbf{b}=\mathbf{Q}_{i}\mathbf{b}=\mathbf{b}$ and

where $\Pi$ denotes the orthogonal projection, $\mathbf{r}^{i}_{\mathbf{x}},\mathbf{r}^{i}_{\mathbf{y}}\in\mathcal{N}_{i-1}$ and $\mathbf{s}^{i}_{\mathbf{x}},\mathbf{s}^{i}_{\mathbf{y}}\in\mathcal{J}_{i-1}^{\perp}$. Since $\dim\ker(\mathbf{A})=n-2m-2\geq n-k-2$ and $\dim\left(\mathcal{N}_{i-1}\right)^{\perp}\geq n-(2i-2)\geq n-2k+2$, we have

Then clearly $\mathbf{U}_{i}\mathbf{b}=\mathbf{b}$, $\mathcal{N}_{i}\subseteq\ker(\mathbf{A})$, and $\dim\mathcal{N}_{i}\leq 2i$. Next, for each $j=0,\dots,i-1$, we have

and similarly $\mathbf{y}^{i}-\mathbf{y}^{0}\in\mathcal{K}_{i-1}(\mathbf{A}_{i};\mathbf{b})\oplus\mathbf{U}_{i}\mathcal{N}_{i}$. This proves (71).

But $\mathbf{Q}_{i}^{\intercal}(\mathbf{y}^{j}-\mathbf{y}^{0})=\mathbf{y}^{j}-\mathbf{y}^{0}$ because $\mathbf{y}^{j}-\mathbf{y}^{0}\in\mathcal{K}_{j-1}(\mathbf{A}_{i-1};\mathbf{b})\oplus\mathbf{U}_{i-1}\mathcal{N}_{i-1}\subseteq\mathcal{J}_{i-1}$, and

which shows that $\nabla_{\mathbf{x}}\mathbf{L}_{i}(\mathbf{x}^{j},\mathbf{y}^{j})=\mathbf{Q}_{i}\mathbf{A}_{i-1}\mathbf{Q}_{i}^{\intercal}(\mathbf{y}^{j}-\mathbf{y}^{0})-\mathbf{b}=\mathbf{A}_{i-1}(\mathbf{y}^{j}-\mathbf{y}^{0})-\mathbf{b}=\nabla_{\mathbf{x}}\mathbf{L}_{i-1}(\mathbf{x}^{j},\mathbf{y}^{j})$. Arguing analogously for the $\mathbf{y}$-variable gives $\nabla_{\mathbf{y}}\mathbf{L}_{i}(\mathbf{x}^{j},\mathbf{y}^{j})=\nabla_{\mathbf{y}}\mathbf{L}_{i-1}(\mathbf{x}^{j},\mathbf{y}^{j})$, proving (70). This completes the induction step, and hence the proof. ∎

Appendix D Experimental details

and $\mathbf{H}=2\mathbf{A}^{\intercal}\mathbf{A}$. Ouyang & Xu (2021) shows that $\|\mathbf{A}\|\leq\frac{1}{2}$, which implies $\|\mathbf{H}\|\leq\frac{1}{2}$. Therefore (25) is a $1$-smooth saddle function.

D.2 Best-iterate gradient norm bound for EG

In Figure 1, we indicated theoretical upper bounds for EG. To clarify, there is no known last-iterate convergence result for EG with respect to $\|\mathbf{G}(\cdot)\|^{2}$. However, it is straightforward to derive $\mathcal{O}(R^{2}/k)$ best-iterate convergence via standard summability arguments in weak convergence proofs for EG. Although there is no theoretical guarantee that $\|\mathbf{G}(\mathbf{z}^{k})\|^{2}$ will monotonically decrease with EG, in our experiments on both examples, they did monotonically decrease (see Figures 1(a), 1(b)). Therefore, we safely used the best-iterate bounds to visualize the upper bound for EG in Figure 1. For the sake of completeness, we derive the best-iterate bound below.

The last inequality is just monotonicity: $\langle\mathbf{z}-\mathbf{z}^{+},\mathbf{w}-\mathbf{z}^{\star}\rangle=\alpha\langle\mathbf{G}(\mathbf{w}),\mathbf{w}-\mathbf{z}^{\star}\rangle\geq 0$. Now the conclusion follows from

where the last inequality follows from $R$-Lipschitzness of $\mathbf{G}$. ∎

Now fix an integer $k\geq 0$, and consider the EG iterations

for $i=0,\dots,k$. Applying Lemma 6 with $\mathbf{z}=\mathbf{z}^{i}$, $\mathbf{w}=\mathbf{z}^{i+1/2}$ and $\mathbf{z}^{+}=\mathbf{z}^{i+1}$, we have

for $i=0,\dots,k$. Summing up the inequalities (72) for all $i=0,\dots,k$, we obtain

The left hand side is at most $\|\mathbf{z}^{0}-\mathbf{z}^{\star}\|^{2}$, while the right hand side is lower bounded by

where $C=\frac{1}{\alpha^{2}(1-\alpha^{2}R^{2})}$.

D.3 ODE flows for 𝐋​(𝒙,𝒚)=𝒙​𝒚𝐋𝒙𝒚𝒙𝒚\boldsymbol{\mathbf{L}(x,y)=xy}

Interestingly, the continuous-time flows with $\mathbf{L}(x,y)=xy$ have exact closed-form solutions.

Note that $\mathbf{G}(x,y)=\begin{bmatrix}0&1\\ -1&0\end{bmatrix}\begin{bmatrix}x\\ y\end{bmatrix}$. Therefore,

The solution to the Moreau–Yosida regularized flow

can be obtained with the matrix exponent. The results are

The anchored flow ODE for $\mathbf{L}(x,y)=xy$ is given by

From the first equation, we have $\frac{d}{dt}(tx(t))=t\dot{x}(t)+x(t)=-ty(t)+x^{0}$, while similar manipulation of the second equation gives $\frac{d}{dt}(ty(t))=tx(t)+y^{0}$. Therefore,

Using the initial conditions to determine the coefficients $c_{1},c_{2}$, we obtain

Appendix E Connection to CLI lower bounds

In this section, we discuss how EAG relates to the prior work on complexity lower bounds on the class of CLI and SCLI algorithms, introduced and studied in (Arjevani et al., 2016; Arjevani & Shamir, 2016; Azizian et al., 2020; Golowich et al., 2020). Specifically, we show that EAG is not SCLI, so it can break the $\Omega(R^{2}/k)$ lower bound on squared gradient norm for the 1-SCLI class derived by Golowich et al. (2020). On the other hand, we show that EAG is 2-CLI in the sense of Golowich et al. (2020), and that EAG belongs to an extended class of 1-CLI algorithms.

We start with the notion of 1-SCLI algorithms by Golowich et al. (2020). Consider an algorithm $\mathcal{A}$ for finding saddle points of biaffine functions of the form

Following the convention of Azizian et al. (2020) and Golowich et al. (2020), we also require that $\mathbf{C},\mathbf{N}$ are matrix polynomials. The classical extragradient method (EG) is an 1-SCLI algorithm: with $\mathbf{G}(\mathbf{z})=\mathbf{B}\mathbf{z}+\mathbf{v}$, we can express EG as

If a 1-SCLI algorithm $\mathcal{A}$ described by (73) is consistent with respect to $\mathbf{B}$, then

Indeed, the 1-SCLI formulation of EG satisfies (74).

For the class of consistent 1-SCLI algorithms, Golowich et al. (2020) established $\Omega(1/k)$ a complexity lower bound on squared gradient norm.

where $\mathbf{z}^{\star}$ is the unique saddle point of $\mathbf{L}$.

To clarify, $\deg\mathbf{N}$ refers to the degree of the matrix polynomial defining $\mathbf{N}$. 1-SCLI algorithms with $d_{\mathbf{C}}=\deg\mathbf{C}=1$ forms a subclass of $\mathfrak{A}_{\textrm{sim}}$ and $\mathfrak{A}_{\textrm{sep}}$. (Even if $d_{\mathbf{C}}>1$, one can still view 1-SCLI algorithms as instances of $\mathfrak{A}_{\textrm{sim}}$ or $\mathfrak{A}_{\textrm{sep}}$ by introducing $d_{\mathbf{C}}-1$ dummy iterates for each 1-SCLI iteration.) However, EAG is an algorithm that belongs to $\mathfrak{A}_{\textrm{sim}}$ but is not 1-SCLI; if it was, a contradiction would occur, as $\|\nabla\mathbf{L}(\mathbf{z}^{k})\|^{2}\leq\mathcal{O}(1/k^{2})$ for EAG. In fact, it is intuitively clear that EAG is not 1-SCLI; the S in 1-SCLI stands for stationary, but EAG has anchoring coefficients $\frac{1}{k+2}$ that vary over iterations.