∎

11institutetext:  Long Chen 22institutetext: Department of Mathematics,
University of California at Irvine, Irvine, CA 92697, USA
22email<EMAIL_ADDRESS>33institutetext: Jingrong Wei 44institutetext:
Department of Mathematics, University of California at Irvine, Irvine, CA
92697, USA
44email<EMAIL_ADDRESS>
# Accelerated Gradient and Skew-Symmetric Splitting Methods for A Class of
Monotone Operator Equations

Long Chen  Jingrong Wei

(Received: date / Accepted: date)

###### Abstract

A class of monotone operator equations, which can be decomposed into sum of a
gradient of a strongly convex function and a linear and skew-symmetric
operator, is considered in this work. Based on discretization of the
generalized gradient flow, gradient and skew-symmetric splitting (GSS) methods
are proposed and proved to convergent in linear rate. To further accelerate
the convergence, an accelerated gradient flow is proposed and accelerated
gradient and skew-symmetric splitting (AGSS) methods are developed, which
extends the acceleration among the existing works on the convex minimization
to a more general class of monotone operator equations. In particular, when
applied to smooth saddle point systems with bilinear coupling, an accelerated
transformed primal-dual (ATPD) method is proposed and shown to achieve linear
rates with optimal lower iteration complexity.

###### Keywords:

Monotone operator accelerated iterative method dynamical system convergence
analysis Lyapunov function saddle point problem

###### MSC:

37N40 47H05 47J25 65B99 65J15 65L20

##  1 Introduction

Consider the nonlinear equation

$\mathcal{A}(x)=0,\quad x\in\mathbb{R}^{n},$ (1)

where $\mathcal{A}$ is strongly monotone: there exists $\mu>0$, for any
$x,y\in\mathbb{R}^{n},$

$(\mathcal{A}(x)-\mathcal{A}(y),x-y)\geq\mu\|x-y\|^{2}.$ (2)

We further assume $\mathcal{A}(x)$ can be decomposed into

$\mathcal{A}(x)=\nabla F(x)+\mathcal{N}x,$ (3)

with a strongly convex $C^{1}$ function $F$, and a linear and skew-symmetric
operator $\mathcal{N}$, i.e., $\mathcal{N}^{\intercal}=-\mathcal{N}$. Then
$\mathcal{A}$ is Lipschitz continuous with constant $L_{\mathcal{A}}\leq
L_{F}+\|\mathcal{N}\|$, where $L_{F}$ is the Lipschitz constant of $\nabla F$.
Therefore $\mathcal{A}$ is monotone and Lipschitz continuous which is also
known as inverse-strongly monotonicity browder1967construction ;
liu1998regularization  or co-coercitivity zhu1996co . Define the condition
numbers:

$\kappa(\mathcal{A})=L_{\mathcal{A}}/\mu,\quad\kappa(F)=L_{F}/\mu,\quad\kappa(\mathcal{N})=\|\mathcal{N}\|/\mu.$

By definition $\kappa(\mathcal{A})\leq\kappa(F)+\kappa(\mathcal{N})$.

Since $\mathcal{A}$ is strongly monotone and Lipschitz continuous, (1) has a
unique solution $x^{\star}$ rockafellar1976monotone . The goal of this paper
is to design iterative methods for computing $x^{\star}$ based on the
splitting (3), which will be called Gradient and skew-Symmetric Splitting
(GSS) methods and Accelerated Gradient and skew-Symmetric Splitting (AGSS)
method for its acceleration.

A fundamental tool to solve (1) is the proximal point method
rockafellar1976monotone . More specific class of monotone operator equation is
particularly of interest, since it is closely related to convex optimization
and saddle point problems. For instance, when $\mathcal{N}=0$, (1) becomes the
convex minimization problem

$\min_{x}F(x).$ (4)

We refer to the book nesterov2003introductory  on the first-order methods for
solving (4) especially the accelerated gradient method which achieved the
linear rate $1-1/\sqrt{\kappa(F)}$ with optimal lower complexity bound.

Another important application is the strongly-convex-strongly-concave saddle
point problem with bilinear coupling:

$\min_{u\in\mathbb{R}^{m}}\max_{p\in\mathbb{R}^{n}}\mathcal{L}(u,p):=f(u)-g(p)+(Bu,p),$
(5)

where $f,g$ are smooth, strongly convex functions with Liptschitz continuous
gradients and $B\in\mathbb{R}^{n\times m},m\geq n,$ is a full rank matrix. By
a sign change, the saddle point system $\nabla\mathcal{L}(u,p)=0$ becomes (1)
with $x=(u,p)^{\intercal}$, $F(u,p)=f(u)+g(p)$ and
$\mathcal{N}=\begin{pmatrix}0&B^{\intercal}\\\ -B&0\end{pmatrix}$. The primal-
dual algorithms for solving (5), such as the augmented Lagrangian method (ALM)
hestenes1969multiplier ; powell1969method , the primal-dual hybrid gradient
(PDHG) method chambolle2011first ; chambolle2016ergodic ; he2012convergence ;
zhu2008efficient  and the alternating direction method of multipliers (ADMM)
gabay1976dual ; glowinski1975approximation , are all based on the proximal
point method applying to the primal or dual problems. Under the strong
monotonicity condition, the accelerated proximal point method developed in
kim2021accelerated , combining restart techniques, has a linear rate. Consider
the pure first-order iterative methods: the sequence
$\left\\{\left(u_{k},p_{k}\right)\right\\}_{k=0,1,\ldots}$ is generated so
that
$\left(u_{k},p_{k}\right)\in\mathcal{H}^{u}_{k}\times\mathcal{H}^{p}_{k}$,
with $\mathcal{H}^{u}_{0}=\operatorname{Span}\left\\{u_{0}\right\\}$,
$\mathcal{H}^{p}_{0}=\operatorname{Span}\left\\{p_{0}\right\\}$, and

$\begin{array}[]{l}\mathcal{H}^{u}_{k+1}:=\operatorname{Span}\left\\{u_{i},\nabla_{u}\mathcal{L}\left(\tilde{u}_{i},\tilde{p}_{i}\right):\forall\tilde{u}_{i}\in\mathcal{H}^{u}_{i},\tilde{p}_{i}\in\mathcal{H}^{p}_{i},0\leq
i\leq k\right\\},\\\
\mathcal{H}^{p}_{k+1}:=\operatorname{Span}\left\\{p_{i},\nabla_{p}\mathcal{L}\left(\tilde{u}_{i},\tilde{p}_{i}\right):\forall\tilde{u}_{i}\in\mathcal{H}^{p}_{i},\tilde{p}_{i}\in\mathcal{H}^{u}_{i},0\leq
i\leq k\right\\}.\end{array}$

In zhang2022lower , it is shown that if the duality gap
$\Delta(u,p):=\max_{q}\mathcal{L}(u,q)-\min_{v}\mathcal{L}(v,p)$ is required
to be bounded by $\epsilon_{\rm out}$, the number of iterations needed is at
least

$\Omega\left(\sqrt{\kappa(f)+\kappa^{2}(\mathcal{N})+\kappa(g)}\cdot|\ln\epsilon_{\rm
out}|\right)$

where $\kappa(f),\kappa(g)$ are condition numbers of $f,g$ respectively,
$\kappa(\mathcal{N})=\|B\|/\sqrt{\mu_{f}\mu_{g}}$ measuring the coupling, and
$\Omega$ means the iteration complexity bounded below asymptotically, up to a
constant. A few optimal first-order algorithms are developed in recent
literature kovalev2021accelerated ; thekumparampil2022lifted . If further the
proximal oracles for $f,g$ are allowed, the lower complexity bound of the
class of first-order algorithms is
$\displaystyle\Omega\left(\kappa(\mathcal{N})|\ln\epsilon_{\rm out}|\right)$.
We refer to chambolle2011first ; chambolle2016ergodic  for an optimal proximal
method.

One more important class is for linear operator $\mathcal{A}$ which has the
following decomposition:

$\mathcal{A}=\mathcal{A}^{\rm sym}+\mathcal{A}^{\rm skew},$

where $\mathcal{A}^{\rm sym}=(\mathcal{A}+\mathcal{A}^{\intercal})/2$ is the
symmetric (Hermitian for complex matrices) part and $\mathcal{A}^{\rm
skew}=(\mathcal{A}-\mathcal{A}^{\intercal})/2$ is the skew-symmetric part. The
condition $\mathcal{A}$ is monotone is equivalent to $\mathcal{A}^{\rm sym}$
is symmetric and positive definite (SPD) and $\lambda_{\min}(\mathcal{A}^{\rm
sym})\geq\mu$. Bai, Golub, and Ng in bai2003hermitian  proposed the
Hermitian/skew-Hermitian splitting method (HSS) for solving general non-
Hermitian positive definite linear systems $\mathcal{A}x=b$:

$\displaystyle(\alpha I+\mathcal{A}^{\rm sym})x_{k+\frac{1}{2}}$
$\displaystyle=(\alpha I-\mathcal{A}^{\rm skew})x_{k}+b$ (6)
$\displaystyle(\alpha I+\mathcal{A}^{\rm skew})x_{k+1}$ $\displaystyle=(\alpha
I-\mathcal{A}^{\rm sym})x_{k+\frac{1}{2}}+b.$

The iterative method (6) solves the equations for the symmetric (Hermitian)
part and skew-symmetric (skew-Hermitian) part alternatively. For the HSS
method (6), efficient solvers for linear operators $(\alpha I+\mathcal{A}^{\rm
sym})^{-1}$ and $(\alpha I+\mathcal{A}^{\rm skew})^{-1}$ are needed. A linear
convergence rate of $\displaystyle\frac{\sqrt{\kappa(\mathcal{A}^{\rm
sym})}-1}{\sqrt{\kappa(\mathcal{A}^{\rm sym})}+1}$ can be achieved for an
optimal choice of parameter $\alpha$. Several variants of the method are
derived and analyzed in bai2007accelerated ; bai2007successive ;
bai2008inexact . Benzi and Golub benzi2004preconditioner  applied the HSS
method for solving the linear saddle point problems with preconditioning
strategy.

For a general non-singular matrix $\mathcal{A}$, $\mathcal{A}^{\rm sym}$ may
not be positive definite. But we can always find a transformation operator
$\mathcal{M}$ s.t. $\mathcal{A}\mathcal{M}$ admits the decomposition (3) and
the transformed system has a similar or better condition number. Indeed the
existence of the solution to the equation $\mathcal{A}x=b$ can be interpreted
as the following inf-sup condition:

$\inf_{y\in\mathbb{R}^{n}}\sup_{x\in\mathbb{R}^{n}}\frac{(\mathcal{A}x,y)}{\|x\|\|y\|}\geq\mu.$
(7)

The inf-sup condition is equivalent to: for any $y\in\mathbb{R}^{n}$, there
exists $x\in\mathbb{R}^{n}$, such that

$(\mathcal{A}x,y)\geq C_{1}\|y\|^{2},\quad\hbox{ and }\;\|x\|\leq C_{2}\|y\|.$
(8)

If we can find a linear operator $\mathcal{M}$ s.t. $x=\mathcal{M}y$
satisfying (8), then $\mathcal{A}\mathcal{M}$ will satisfy (3) with the
symmetric and skew-symmetric decomposition. A large class of iterative methods
for solving the linear saddle point systems can be derived by different
choices of $\mathcal{M}$; see Benzi.M;Golub.G;Liesen.J2005  and the references
therein. A transformation has also been proposed for nonlinear saddle point
problems recently in ChenWei2022Transformed  and the resulting algorithm is
called the transformed primal-dual (TPD) method.

We shall develop our iterative methods based on discretization of ODE flows.
Namely we treat $x(t)$ as a continuous function of an artificial time variable
$t$ and design ODE systems $x^{\prime}(t)=\mathcal{G}(x(t))$ so that the
$x^{\star}$ is a stable equilibrium point of the corresponding dynamic system,
i.e., $\mathcal{G}(x^{*})=0$ and $\lim_{t\to\infty}x(t)=x^{*}$. Then we apply
ODE solvers to obtain various iterative methods. By doing this way, we can
borrow the stability analysis for dynamic systems and convergence theory of
ODE solvers. In particular, the dynamic systems and accelerated iterative
methods for the convex minimization (4) is analyzed via a unified Lyapunov-
based approach in chen2021unified ; see also attouch2019convergence  for
maximally monotone operators. However, there is few investigation on the ODE
flows for (1), and iterative methods with accelerated linear rates.

The most popular choice is the generalized gradient flow:

$x^{\prime}(t)=-\mathcal{A}(x(t)).$ (9)

Although $\mathcal{A}$ is not the gradient of some scalar function when
$\mathcal{N}\neq 0$, we still use the name convention for the convex case
$\mathcal{N}=0$. One can easily show the exponential stability of $x^{\star}$
using the monotonicity of $\mathcal{A}$. Discretization for the generalized
gradient flow (9) leads to iteration methods for solving (1). The implicit
Euler scheme with step size $\alpha_{k}$ is the unconditionally stable
proximal method $(I+\alpha_{k}\mathcal{A})(x_{k+1})=x_{k}$. The explicit Euler
scheme gives the generalized gradient descent method:

$x_{k+1}=x_{k}-\alpha_{k}\mathcal{A}(x_{k}).$ (10)

For $\displaystyle\alpha_{k}=\mu/L_{\mathcal{A}}^{2}$, we have the linear
convergence

$\|x_{k+1}-x^{\star}\|^{2}\leq\left(1-1/\kappa^{2}(\mathcal{A})\right)\|x_{k+1}-x^{\star}\|^{2}.$

But this linear rate is pretty slow when $\kappa(\mathcal{A})\gg 1$ since the
dependence is $\kappa^{2}(\mathcal{A})$. To achieve the accuracy
$\|x_{k+1}-x^{\star}\|^{2}\leq\epsilon_{\rm out}$,
$\mathcal{O}(\kappa^{2}(\mathcal{A})|\ln\epsilon_{\rm out}|)$ iterations are
needed. In contrast, for convex optimization, i.e., $\mathcal{N}=0$, the rate
for the gradient descent method is $1-1/\kappa(F)$ for step size
$\alpha_{k}=1/L_{F}$.

To speed up the convergence, we explore more property from the decomposition
(3). Let $B^{\intercal}={\rm upper}(\mathcal{N})$ be the upper triangular part
of $\mathcal{N}$ and $B^{\rm sym}=B+B^{\intercal}$ be a symmetrization of $B$.
Based on the splitting $\mathcal{N}=B^{\rm sym}-2B$ or
$\mathcal{N}=2B^{\intercal}-B^{\rm sym}$, we develop Gradient and skew-
Symmetric Splitting (GSS) methods:

$\frac{x_{k+1}-x_{k}}{\alpha}=-\left(\nabla F(x_{k})+B^{\rm
sym}x_{k}-2Bx_{k+1}\right),$ (11)

or

$\frac{x_{k+1}-x_{k}}{\alpha}=-\left(\nabla F(x_{k})-B^{\rm
sym}x_{k}+2B^{\intercal}x_{k+1}\right).$ (12)

Notice that $B$ is lower triangular, both (11) and (12) are explicit schemes
as $x_{k+1}$ can be computed via forward (for (11)) or backward (for (12))
substitution. We prove that they achieve the linear convergence rate

$\|x_{k}-x^{\star}\|^{2}\leq\left(\frac{1}{1+4/\max\\{\kappa(B^{\rm
sym}),\kappa(F)\\}}\right)^{k}6\|x_{0}-x^{\star}\|^{2}.$

where $\kappa(B^{\rm sym})=L_{B^{\rm sym}}/\mu$ with $L_{B^{\rm sym}}=\|B^{\rm
sym}\|$ for $\alpha=\displaystyle\min\left\\{\frac{1}{4L_{B^{\rm
sym}}},\frac{1}{4L_{F}}\right\\}$. In particular for the linear problem
$\nabla F=\mu I$, this can be viewed as accelerated overrelaxation (AOR)
methods Hadjidimos:1978Accelerated ; hadjidimos2000successive  for a linear
shifted skew-symmetric system.

To further accelerate the convergence rate, following
luoDifferentialEquationSolvers2021a , we introduce an accelerated gradient
flow

$\left\\{\begin{aligned} x^{\prime}&=y-x,\\\
y^{\prime}&=x-y-\frac{1}{\mu}(\nabla F(x)+\mathcal{N}y).\end{aligned}\right.$
(13)

Comparing with the accelerated gradient flow in
luoDifferentialEquationSolvers2021a  for convex optimization, the difference
is the gradient and skew-symmetric splitting $\mathcal{A}(x)\rightarrow\nabla
F(x)+\mathcal{N}y$.

We propose several iterative schemes based on various discretization of (13).
Provided that a fast solver for computing the shifted skew-symmetric operator
$(\beta I+\mathcal{N})^{-1}$ is available, we give an implicit-explicit (IMEX)
Euler scheme for (13) as a version of Accelerated Gradient and skew-Symmetric
Splitting (AGSS) methods:

$\displaystyle\frac{\hat{x}_{k+1}-x_{k}}{\alpha_{k}}$
$\displaystyle=y_{k}-\hat{x}_{k+1},$ (14)
$\displaystyle\frac{y_{k+1}-y_{k}}{\alpha_{k}}$
$\displaystyle=\hat{x}_{k+1}-y_{k+1}-\frac{1}{\mu}\left(\nabla
F(\hat{x}_{k+1})+\mathcal{N}y_{k+1}\right),$
$\displaystyle\frac{x_{k+1}-x_{k}}{\alpha_{k}}$
$\displaystyle=y_{k+1}-x_{k+1}.$

The scheme (14) is implicit for $y_{k+1}$ as each iteration needs to solve a
shifted skew-symmetric equation $(\beta
I+\mathcal{N})y_{k+1}=b(\hat{x}_{k+1},y_{k})$ with $\beta=1+\mu/\alpha_{k}$.
For the fixed step size $\alpha_{k}=1/\sqrt{\kappa(F)}$, the convergence rate
is accelerated to

$\|x_{k+1}-x^{\star}\|^{2}+\|y_{k+1}-x^{\star}\|^{2}\leq\left(\frac{1}{1+1/\sqrt{\kappa(F)}}\right)^{k}2\mathcal{E}_{0}/\mu,$

where
$\mathcal{E}_{0}=D_{F}(x_{0},x^{\star})+\frac{\mu}{2}\|y_{0}-x^{\star}\|^{2}$
and $D_{F}$ is the Bregman divergence of $F$. We also allow an inexact solver
for approximating $(\beta I+\mathcal{N})^{-1}$ using perturbation argument to
control the inner solve error. For linear systems, comparing with HSS, we
achieve the same accelerated rate without treating the symmetric part
implicitly, i.e., no need to compute $(\alpha I+\mathcal{A}^{\rm sym})^{-1}$.

To fully avoid computing $(\beta I+\mathcal{N})^{-1}$, we provide an explicit
AGSS scheme combining acceleration of $F$ and AOR technique for $\mathcal{N}$:

$\displaystyle\frac{\hat{x}_{k+1}-x_{k}}{\alpha}$
$\displaystyle=y_{k}-\hat{x}_{k+1},$ (15)
$\displaystyle\frac{y_{k+1}-y_{k}}{\alpha}$
$\displaystyle=\hat{x}_{k+1}-y_{k+1}-\frac{1}{\mu}\left(\nabla
F(\hat{x}_{k+1})+B^{\rm sym}y_{k}-2By_{k+1}\right),$
$\displaystyle\frac{x_{k+1}-x_{k}}{\alpha}$
$\displaystyle=y_{k+1}-\frac{1}{2}(x_{k+1}+\hat{x}_{k+1}).$

For the step size $\displaystyle\alpha=\min\left\\{\frac{\mu}{2L_{B^{\rm
sym}}},\sqrt{\frac{\mu}{2L_{F}}}\right\\}$, we obtain the accelerated linear
convergence for scheme (15)

$\|x_{k+1}-x^{\star}\|^{2}\leq\left(\frac{1}{1+1/\max\\{4\kappa(B^{\rm
sym}),\sqrt{8\kappa(F)}\\}}\right)^{k}\frac{2\mathcal{E}_{0}^{\alpha B}}{\mu}$

where $\mathcal{E}^{\alpha
B}_{0}=D_{F}(x_{0},x^{\star})+\frac{1}{2}\|y-x^{\star}\|^{2}_{\mu I-2\alpha
B^{\rm sym}}$ is nonnegative according to our choice of the step size.

The proposed accelerated methods has wide applications. This extends the
acceleration among the existing works on the convex minimization guler1992new
; nesterov1983method ; polyak1964some  to a more general class of monotone
operator equations. As an example, we apply to smooth strongly-convex-
strongly-concave saddle point systems with bilinear coupling (5) and derive
first-order algorithms achieving optimal lower iteration complexity
$\Omega\left(\sqrt{\kappa(f)+\kappa^{2}(\mathcal{N})+\kappa(g)}\cdot|\ln\epsilon_{\rm
out}|\right)$; see Section 5.2 for variants and details.

For affinely constrained optimization problems:

$\displaystyle\min_{u\in\mathbb{R}^{m}}f(u)$ $\displaystyle\text{subject
to}\quad Bu=b,$

we combine the transformed primal-dual method developed in
ChenWei2022Transformed  and the accelerated gradient flow to propose an
accelerated transformed primal-dual (ATPD) method:

$\displaystyle\frac{\hat{x}_{k+1}-x_{k}}{\alpha}$
$\displaystyle=y_{k}-\hat{x}_{k+1},$ (16)
$\displaystyle\frac{v_{k+1}-v_{k}}{\alpha}$
$\displaystyle=\frac{1}{2}(\hat{u}_{k+1}-v_{k+1})-\frac{1}{\mu_{f}}\mathcal{I}_{\mathcal{V}}^{-1}(\nabla
f(\hat{u}_{k+1})+B^{\intercal}q_{k}),$
$\displaystyle\frac{q_{k+1}-q_{k}}{\alpha}$
$\displaystyle=\hat{p}_{k+1}-q_{k+1}-\mathcal{I}_{\mathcal{Q}}^{-1}(Sp_{k+1}+Bv_{k}-2Bv_{k+1}+B\mathcal{I}_{\mathcal{V}}^{-1}f(u_{k+1})),$
$\displaystyle\frac{x_{k+1}-\hat{x}_{k+1}}{\alpha}$
$\displaystyle=(y_{k+1}-y_{k})-\frac{1}{4}(x_{k+1}-\hat{x}_{k+1}).$

where $x=(u,p)$, $y=(v,q)$,
$\mathcal{I}_{\mathcal{V}},\mathcal{I}_{\mathcal{Q}}$ are SPD operators and
$S=B\mathcal{I}_{\mathcal{V}}^{-1}B^{\intercal}$. The complexity of the
explicit scheme (16) is

$\mathcal{O}(\sqrt{\kappa_{\mathcal{I}_{\mathcal{V}}}(f)}\,\kappa_{\mathcal{I}_{\mathcal{Q}}}(S)|\ln\epsilon_{\rm
out}|)$

which suggests that $\mathcal{I}_{\mathcal{Q}}^{-1}$ is better chosen as a
preconditioner for the Schur complement
$S=B\mathcal{I}_{\mathcal{V}}^{-1}B^{\intercal}$. In particular, when
$\mathcal{I}_{\mathcal{V}}$ is a scaled identity and
$\mathcal{I}_{\mathcal{Q}}=S$, the outer iteration complexity is
$\mathcal{O}(\sqrt{\kappa_{\mathcal{I}_{\mathcal{V}}}(f)}|\ln\epsilon_{\rm
out}|)$. Iterative methods for computing $S^{-1}r$ can be thought of as an
inner iteration which can be a linear inner solver or a non-linear one. The
linear solver will enter the estimate by the factor
$\kappa_{\mathcal{I}_{\mathcal{Q}}}(S)$. The non-linear inner solver for
$S^{-1}$, e.g., the conjugate gradient method, can be estimated by the
perturbation argument on controlling the residual $\varepsilon_{\rm in}$
developed in Section 4.4. The total iteration complexity for the accelerated
transformed primal-dual method is

$\mathcal{O}\left(|\ln\epsilon_{\rm out}||\ln\varepsilon_{\rm
in}|\sqrt{\kappa_{\mathcal{I}_{\mathcal{V}}}(f)\,\kappa_{\mathcal{I}_{\mathcal{Q}}}(S)}\right),$

which achieves the optimal complexity bound for affinely constrained problems
salim2022optimal .

The strong convexity of $f$ can be further relaxed to

$f_{\beta}(u)=f(u)+\frac{\beta}{2}\|Bu-b\|^{2}$

by using augmented Lagrangian method (ALM). Then $\mu_{f}>0$ can be relaxed to
$\mu_{f_{\beta}}>0$; see discussion in Section 5.3.

For the rest of the paper, we first review convex function theory and Lyapunov
analysis theory in Section 2, as basic preliminaries for our proofs. The
generalized gradient flow and convergence of GSS methods for the considered
monotone equation are proposed in Section 3. Then in Section 4, we extend the
accelerated gradient flow in convex minimization and derive AGSS methods with
accelerated linear rates. As applications, we propose accelerated first-order
algorithms for convex-concave saddle point problems with bilinear coupling in
Section 5. In particular, we give optimal algorithms for strongly-convex-
strongly-concave saddle point systems. Concluding remarks are addressed in
Section 6.

##  2 Preliminary

In this section, we follow nesterov2003introductory  to introduce notation and
preliminaries in convex function theory. We also briefly review a unified
convergence analysis of first order convex optimization methods via strong
Lyapunov functions established in chen2021unified .

###  2.1 Convex function

Let $\mathcal{X}$ be a finite-dimensional Hilbert space with inner product
$(\cdot,\cdot)$ and norm $\|\cdot\|$. $\mathcal{X}^{\prime}$ is the linear
space of all linear and continuous mappings
$T:\mathcal{X}\rightarrow\mathbb{R}$, which is called the dual space of
$\mathcal{X}$, and $\langle\cdot,\cdot\rangle$ denotes the duality pair
between $\mathcal{X}$ and $\mathcal{X}^{\prime}$. For any proper closed convex
and $C^{1}$ function $F:\mathcal{X}\rightarrow\mathbb{R}$, the Bregman
divergence of $F$ is defined as

$D_{F}(x,y):=F(x)-F(y)-\langle\nabla F(y),x-y\rangle.$

We say $F\in\mathcal{S}_{\mu}$ or $\mu$-strongly convex with $\mu\geqslant 0$
if $F$ is differentiable and

$D_{F}(x,y)\geqslant\frac{\mu}{2}\left\|x-y\right\|^{2},\quad\forall
x,y\in\mathcal{X}.$

In addition, denote $F\in\mathcal{S}_{\mu,L}$ if $F\in\mathcal{S}_{\mu}$ and
there exists $L>0$ such that

$D_{F}(x,y)\leqslant\frac{L}{2}\left\|x-y\right\|^{2},\quad\forall
x,y\in\mathcal{X}.$

For fixed $y\in\mathcal{X},D_{F}(\cdot,y)$ is convex as $F$ is convex and

$\nabla D_{F}(\cdot,y)=\nabla F(\cdot)-\nabla F(y).$

For $F\in\mathcal{S}_{\mu,L}$, we have

$\frac{\mu}{2}\|x-y\|^{2}\leq D_{F}(x,y)\leq\frac{L}{2}\|x-y\|^{2}.$

Especially for $F(x)=\frac{1}{2}\|x\|^{2}$, Bregman divergence reduces to the
half of the squared distance $D_{F}(x,y)=D_{F}(y,x)=\frac{1}{2}\|x-y\|^{2}$.
In general $D_{F}(x,y)$ is non-symmetric in terms of $x$ and $y$. A
symmetrized Bregman divergence is defined as

$\langle\nabla F(y)-\nabla F(x),y-x\rangle=D_{F}(y,x)+D_{F}(x,y).$

By direct calculation, we have the following three-terms identity.

######  Lemma 1 (Bregman divergence identity chen1993convergence )

If function $F:\mathcal{X}\rightarrow\mathbb{R}$ is differentiable, then for
any $x,y,z\in\mathcal{X}$, it holds that

$\langle\nabla F(x)-\nabla F(y),y-z\rangle=D_{F}(z,x)-D_{F}(z,y)-D_{F}(y,x).$
(17)

When $F(x)=\frac{1}{2}\|x\|^{2}$, identity (17) becomes

$(x-y,y-z)=\frac{1}{2}\|z-x\|^{2}-\frac{1}{2}\|z-y\|^{2}-\frac{1}{2}\|y-x\|^{2}.$

For a non-smooth function $F$, the subgradient is a set-valued function
defined as

$\partial
F(x):=\left\\{\xi\in\mathcal{X}^{\prime}:F(y)-F(x)\geqslant\langle\xi,y-x\rangle,\quad\forall
y\in\mathcal{X}\right\\}.$

For a proper, closed and convex function $F:\mathcal{X}\to\mathbb{R}$,
$\partial F(x)$ is a nonempty bounded set for $x\in\mathcal{X}$
nesterov2003introductory . It is evident that the subgradient for smooth $F$
is a single-valued function reduced to the gradient, that is $\partial
F(x)=\\{\nabla F(x)\\}$. We extend the set $\mathcal{S}_{\mu}$ using the
notion of subgradient: $F\in\mathcal{S}_{\mu}$ if for all $x,y\in\mathcal{X}$,

$F(y)\geqslant
F(x)+\langle\xi,y-x\rangle+\frac{\mu}{2}\|x-y\|^{2},\quad\forall\xi\in\partial
F(x).$

Similar to the smooth case, if $F\in\mathcal{S}_{\mu}$, then

$\langle\xi-\eta,x-y\rangle\geqslant\mu\|x-y\|^{2},\quad\forall
x,y\in\mathcal{X},\xi\in\partial F(x),\eta\in\partial F(y).$

Given a convex function $F$, define the proximal operator of $F$ as

$\operatorname{prox}_{\gamma
F}(y):=\underset{x}{\operatorname{argmin}}~{}F(x)+\frac{1}{2\gamma}\|x-y\|^{2},$
(18)

for some $\gamma>0$. The proximal operator is well-defined since the function
$F(\cdot)+\frac{1}{2\gamma}\|\cdot-y\|^{2}$ is strongly convex.

###  2.2 Inner product and preconditioners

The standard $l^{2}$ dot product of Euclidean space $(\cdot,\cdot)$ is usually
chosen as the inner product and the norm induced is the Euclidean norm. We now
introduce inner product $(\cdot,\cdot)_{\mathcal{I}_{\mathcal{X}}}$ induced by
a given SPD operator $\mathcal{I}_{\mathcal{X}}:\mathcal{X}\to\mathcal{X}$
defined as follows

$(x,y)_{\mathcal{I}_{\mathcal{X}}}:=(\mathcal{I}_{\mathcal{X}}x,y)=(x,\mathcal{I}_{\mathcal{X}}y),\quad\forall
x,y\in\mathcal{X}$

and associated norm $\|\cdot\|_{\mathcal{I}_{\mathcal{X}}}$, given by
$\|x\|_{\mathcal{I}_{\mathcal{X}}}=(x,x)^{1/2}_{\mathcal{I}_{\mathcal{X}}}.$
The dual norm w.r.t the $\mathcal{I}_{\mathcal{X}}$-norm is defined as: for
$\ell\in\mathcal{X}^{\prime}$

$\|\ell\|_{\mathcal{X}^{\prime}}=\sup_{0\neq
x\in\mathcal{X}}\frac{\langle\ell,x\rangle}{\|x\|_{\mathcal{I}_{\mathcal{X}}}}=\left(\ell,\ell\right)_{\mathcal{I}_{\mathcal{X}}^{-1}}^{1/2}=\left(\mathcal{I}_{\mathcal{X}}^{-1}\ell,\ell\right)^{1/2}.$

We shall generalize the convexity and Lipschitz continuity with respect to
$\mathcal{I}_{\mathcal{X}}$-norm: we say
$F\in\mathcal{S}_{\mu_{F,\mathcal{I}_{\mathcal{X}}}}$ with
$\mu_{F,\mathcal{I}_{\mathcal{X}}}\geqslant 0$ if $F$ is differentiable and

$D_{F}(x,y)\geqslant\frac{\mu_{F,\mathcal{I}_{\mathcal{X}}}}{2}\left\|x-y\right\|^{2}_{\mathcal{I}_{\mathcal{X}}},\quad\forall
x,y\in\mathcal{X}.$

In addition, denote
$f\in\mathcal{S}_{\mu_{F,\mathcal{I}_{\mathcal{X}}},L_{F,\mathcal{I}_{\mathcal{X}}}}$
if $f\in\mathcal{S}_{\mu_{F,\mathcal{I}_{\mathcal{X}}}}$ and there exists
$L_{f,\mathcal{I}_{\mathcal{X}}}>0$ such that

$D_{F}(x,y)\leq\frac{L_{f,\mathcal{I}_{\mathcal{X}}}}{2}\left\|x-y\right\|^{2}_{\mathcal{I}_{\mathcal{X}}},\quad\forall
x,y\in\mathcal{X}.$

The gradient method in the inner product $\mathcal{I}_{\mathcal{X}}$ reads as:

$x_{k+1}=x_{k}-\mathcal{I}_{\mathcal{X}}^{-1}\nabla F(x_{k}),$ (19)

where $\mathcal{I}_{\mathcal{X}}^{-1}$ can be understood as a preconditioner.
The convergence rate of the preconditioned gradient descent iteration (19) is
$1-1/\kappa_{\mathcal{I}_{\mathcal{X}}}(F)$ with condition number
$\kappa_{\mathcal{I}_{\mathcal{X}}}(F)=L_{F,\mathcal{I}_{\mathcal{X}}}/\mu_{F,\mathcal{I}_{\mathcal{X}}}$.
To simplify notation, we skip $\mathcal{I}_{\mathcal{X}}$ in the constants
$\mu$ and $L$, e.g., write $\mu_{F,\mathcal{I}_{\mathcal{X}}}$ as $\mu_{F}$,
but keep in the condition number, e.g.
$\kappa_{\mathcal{I}_{\mathcal{X}}}(F)$, to emphasize that the condition
number is measured in the $\mathcal{I}_{\mathcal{X}}$ inner product.

For two symmetric operators $A,B:\mathcal{X}\to\mathcal{X}$, we say
$A\geq(>)B$ if $A-B$ is positive semidefinite (definite). Therefore $A>0$
means $A$ is SPD. One can easily verify that $A\geq B$ is equivalent to
$\lambda_{\min}(B^{-1}A)\geq 1$ or $\lambda_{\max}(A^{-1}B)\leq 1$ for two
non-singular and symmetric operators.

###  2.3 Lyapunov analysis

In order to study the stability of an equilibrium $x^{*}$ of a dynamical
system defined by an autonomous system

$x^{\prime}=\mathcal{G}(x(t)),$ (20)

Lyapunov introduced the so-called Lyapunov function $\mathcal{E}(x)$
Khalil:1173048 ; Haddad2008 , which is nonnegative and the equilibrium point
$x^{*}$ satisfies $\mathcal{E}\left(x^{*}\right)=0$ and the Lyapunov
condition: $-\nabla\mathcal{E}(x)\cdot\mathcal{G}(x)>0$ for $x$ near the
equilibrium point $x^{*}$. That is the flow $\mathcal{G}(x)$ may not be in the
perfect $-\nabla\mathcal{E}(x)$ direction but contains positive component in
that direction. Then the (local) decay property of $\mathcal{E}(x)$ along the
trajectory $x(t)$ of the autonomous system (20) can be derived immediately:

$\frac{\mathrm{d}}{\mathrm{d}t}\mathcal{E}(x(t))=\nabla\mathcal{E}(x)\cdot
x^{\prime}(t)=\nabla\mathcal{E}(x)\cdot\mathcal{G}(x)<0.$

To further establish the convergence rate of $\mathcal{E}(x(t))$, Chen and Luo
chen2021unified  introduced the strong Lyapunov condition: $\mathcal{E}(x)$ is
a Lyapunov function and there exist constant $q\geqslant 1$, strictly positive
function $c(x)$ and function $p(x)$ such that

$-\nabla\mathcal{E}(x)\cdot\mathcal{G}(x)\geq c(x)\mathcal{E}^{q}(x)+p^{2}(x)$

holds true near $x^{*}$. From this, one can derive the exponential decay
$\mathcal{E}(x(t))=O\left(e^{-ct}\right)$ for $q=1$ and the algebraic decay
$\mathcal{E}(x(t))=O\left(t^{-1/(q-1)}\right)$ for $q>1$. Furthermore if
$\|x-x^{*}\|^{2}\leq C\mathcal{E}(x)$, then we can derive the exponential
stability of $x^{*}$ from the exponential decay of Lyapunov function
$\mathcal{E}(x)$.

Note that for an optimization problem, we have freedom to design the vector
field $\mathcal{G}(x)$ and Lyapunov function $\mathcal{E}(x)$.

##  3 Gradient and skew-Symmetric Splitting Methods

In this section, we shall consider the generalized gradient flow

$x^{\prime}(t)=-\mathcal{A}(x(t)),$ (21)

and derive several iterative methods for solving $\mathcal{A}(x)=0$ by
applying ODE solvers to (21). Based on the gradient and skew-symmetric
splitting (3), we apply Accelerated OverRelaxation (AOR) technique to the non-
symmetric part to obtain explicit schemes with linear convergence rate
$(1+c/\kappa(\mathcal{A}))^{-1}$.

###  3.1 Stability

Define the Lyapunov function:

$\mathcal{E}_{q}(x)=\frac{1}{2}\|x-x^{\star}\|^{2}.$ (22)

For $x(t)$ solving the gradient flow (21):

$\displaystyle\frac{\,{\rm d}}{\,{\rm d}t}\mathcal{E}_{q}(x(t))$
$\displaystyle=\langle\nabla\mathcal{E}_{q}(x),x^{\prime}\rangle=-\langle
x-x^{\star},\mathcal{A}(x)-\mathcal{A}(x^{\star})\rangle\leq-2\mu\mathcal{E}_{q}(x),$

which leads to the exponential decay

$\mathcal{E}_{q}(x(t))\leq e^{-2\mu t}\mathcal{E}_{q}(x(0))$

and consequently the solution $x^{*}$ is an exponentially stable equilibrium
point of the dynamic system defined by (21).

###  3.2 Implicit Euler Schemes

The implicit Euler scheme for the generalized gradient flow (21) with step
size $\alpha_{k}>0$ is

$x_{k+1}-x_{k}=-\alpha_{k}\mathcal{A}(x_{k+1}),$ (23)

which is equivalent to the proximal iteration

$(I+\alpha_{k}\mathcal{A})(x_{k+1})=x_{k}.$ (24)

As $\mathcal{A}$ is strongly monotone and Liptschitz continuous, the operator
$(I+\alpha_{k}\mathcal{A})^{-1}$, which is called the resolvent of
$\mathcal{A}$, is nonexpansive browder1967construction  and has a unique fixed
point rockafellar1976monotone .

Using the identity of squares and the strong monotonicity, we have

$\displaystyle\mathcal{E}_{q}(x_{k+1})-\mathcal{E}_{q}(x_{k})={}$
$\displaystyle(x_{k+1}-x_{k},x_{k+1}-x^{\star})-\frac{1}{2}\|x_{k+1}-x_{k}\|^{2}$
(25) $\displaystyle={}$
$\displaystyle-\alpha_{k}\left\langle{\mathcal{A}(x_{k+1})-\mathcal{A}(x^{\star}),x_{k+1}-x^{\star}}\right\rangle-\frac{1}{2}\|x_{k+1}-x_{k}\|^{2}$
$\displaystyle\leq{}$ $\displaystyle-2\alpha_{k}\mu\mathcal{E}_{q}(x_{k+1}),$

from which the linear convergence follows naturally

$\|x_{k+1}-x^{\star}\|^{2}\leq\frac{1}{1+2\alpha_{k}\mu}\|x_{k}-x^{\star}\|^{2}.$

There is no restriction on the step size which is known as the unconditional
stability of the implicit Euler method. Choosing $\alpha_{k}\gg 1$ will
accelerate the convergence with the price of solving a regularized nonlinear
problem (24), which is usually not practical. The implicit Euler method is
present here as a reference to measure the difference of other methods.

###  3.3 Explicit Euler Schemes

The explicit Euler scheme for the generalized gradient flow (21) leads to the
gradient decent method:

$x_{k+1}=x_{k}-\alpha_{k}\mathcal{A}(x_{k}).$ (26)

We write (26) as a correction of the implicit Euler scheme (23), i.e.,

$x_{k+1}-x_{k}=-\alpha_{k}\mathcal{A}(x_{k+1})+\alpha_{k}(\mathcal{A}(x_{k+1})-\mathcal{A}(x_{k})),$

and follow (25) to obtain

$\displaystyle\mathcal{E}_{q}(x_{k+1})-\mathcal{E}_{q}(x_{k})=$
$\displaystyle(x_{k+1}-x_{k},x_{k+1}-x^{\star})-\frac{1}{2}\|x_{k+1}-x_{k}\|^{2}$
(27) $\displaystyle\leq$
$\displaystyle-2\alpha_{k}\mu\mathcal{E}_{q}(x_{k+1})-\frac{1}{2}\|x_{k+1}-x_{k}\|^{2}$
$\displaystyle-\alpha_{k}\left\langle{\mathcal{A}(x_{k})-\mathcal{A}(x_{k+1}),x_{k+1}-x^{\star}}\right\rangle.$

Here we use the implicit Euler scheme as the reference scheme and call the
term $\alpha_{k}(\mathcal{A}(x_{k})-\mathcal{A}(x_{k+1}),x_{k+1}-x^{\star})$ a
mis-matching term. Using the Cauchy-Schwarz inequality and the Lipschitz
continuity of $\mathcal{A}$,

$\displaystyle\alpha_{k}|\left\langle{\mathcal{A}(x_{k})-\mathcal{A}(x_{k+1}),x_{k+1}-x^{\star}}\right\rangle|$
$\displaystyle\leq{}$
$\displaystyle\frac{\alpha_{k}}{2\mu}\|\mathcal{A}(x_{k})-\mathcal{A}(x_{k+1})\|^{2}+\frac{\alpha_{k}\mu}{2}\|x_{k+1}-x^{\star}\|^{2}$
$\displaystyle\leq{}$
$\displaystyle\frac{\alpha_{k}L_{\mathcal{A}}^{2}}{2\mu}\|x_{k+1}-x_{k}\|^{2}+\alpha_{k}\mu\,\mathcal{E}_{q}(x_{k+1})$

Substitute back to (27) and let
$\displaystyle\alpha_{k}\leq\frac{\mu}{L_{\mathcal{A}}^{2}}$, we have the
linear reduction

$\mathcal{E}_{q}(x_{k+1})\leq\frac{1}{1+\alpha_{k}\mu}\mathcal{E}_{q}(x_{k}).$

In particular when $\displaystyle\alpha_{k}=\frac{\mu}{L_{\mathcal{A}}^{2}}$,

$\mathcal{E}_{q}(x_{k+1})\leq\frac{1}{1+1/\kappa^{2}(\mathcal{A})}\mathcal{E}_{q}(x_{k}).$

###### Remark 1

Using another identity on squares

$\mathcal{E}_{q}(x_{k+1})-\mathcal{E}_{q}(x_{k})=(x_{k+1}-x_{k},x_{k}-x^{\star})+\frac{1}{2}\|x_{k+1}-x_{k}\|^{2},$

and the monotonicity at $x_{k}$, the upper bound of the step size can be
relaxed to $\displaystyle\alpha_{k}<2\mu/L_{\mathcal{A}}^{2}$ and for
$\alpha_{k}=\mu/L_{\mathcal{A}}^{2}$, the rate can be improved to
$1-1/\kappa^{2}(\mathcal{A})$. Here we emphasize the approach of using the
implicit Euler method as the reference and the estimate of the mis-match
terms. The obtained linear rate is slightly worse but still in the same order
of $\kappa(\mathcal{A})$ as
$(1+\delta)^{-1}=1-\delta+\mathcal{O}(\delta^{2})$.

###  3.4 Accelerated Overrelaxation Methods for Shifted skew-Symmetric
Problems

In this subsection, we consider the shifted skew-symmetric linear system

$(\mu I+\mathcal{N})x=b,$ (28)

which is a special case of (1) with $F(x)=\frac{\mu}{2}\|x\|^{2}-(b,x)$ and
$\mathcal{A}(x)=(\mu I+\mathcal{N})x-b$. An efficient solver for (28) is an
ingredient of the accelerated gradient and skew-symmetric splitting method for
the nonlinear equation (1) we shall develop later on. The system (28) can be
solved with Krylov subspace methods. For example, minimal residual methods
idema2007minimal ; jiang2007algorithm  based on the Lanczos algorithm require
short recurrences of vector storage. The convergence theorem using Chebyshev
polynomial jiang2007algorithm  shows that the residual converges with a linear
rate of $O((1+\mu/\|\mathcal{N}\|)^{-1})$. Compared with Krylov subspace
methods for solving (28), our method to be presented enjoys a similar linear
convergence rate but with a much simpler form and in-space vector storage.

Recall that, for the matrix equation

$(D-L-U)x=b,$

where $D$ a diagonal matrix, and $L$ and $U$ are strictly lower and upper
triangular matrix, respectively, the Accelerated OverRelaxation (AOR) method
Hadjidimos:1978Accelerated ; hadjidimos2000successive  with the acceleration
parameter $r$ and relaxation factor $\omega>1$ is in the form

$(D-rL)x_{k+1}=[(1-\omega)D+(\omega-r)L+\omega U]x_{k}+\omega b,$

which is a two-parameter generalization of the successive overrelaxation (SOR)
method. We shall apply AOR to (28) with well-designed $r$ and $\omega$.

As $\mathcal{N}=-\mathcal{N}^{\intercal}$, all diagonal entries of
$\mathcal{N}$ are zero. Let $B^{\intercal}={\rm upper}(\mathcal{N})$ be the
upper triangular part of $\mathcal{N}$. Then $\mathcal{N}=B^{\intercal}-B$.
Let $B^{\rm sym}=B+B^{\intercal}$ be a symmetrization of $B$. We have the
splitting

$\displaystyle\mathcal{N}$ $\displaystyle=B^{\rm sym}-2B,\quad\text{and}$ (29)
$\displaystyle\mathcal{N}$ $\displaystyle=2B^{\intercal}-B^{\rm sym}.$ (30)

As a symmetric matrix, all eigenvalues of $B^{\rm sym}$ are real. Let
$\lambda_{1}\leq\cdots\leq\lambda_{n}$ be the eigenvalues of $B^{\rm sym}$.
Since ${\rm trace}(B^{\rm sym})=0$, we know $\lambda_{1}<0$ and
$\lambda_{n}>0$. That is $B^{\rm sym}$ is symmetric but indefinite. Define
$L_{B^{\rm sym}}=\|B^{\rm sym}\|=\max\\{|\lambda_{1}|,\lambda_{n}\\}$ and
$\kappa(B^{\rm sym})=L_{B^{\rm sym}}/\mu$.

We discretize the gradient flow (21) by the matrix splitting (29) and treat
$B^{\rm sym}$ as the explicit term:

$\frac{x_{k+1}-x_{k}}{\alpha}=-\left(\mu x_{k+1}+B^{\rm
sym}x_{k}-2Bx_{k+1}-b\right).$ (31)

The update of $x_{k+1}$ is equivalent to solve

$\left[(1+\alpha\mu)I-2\alpha B\right]x_{k+1}=x_{k}-\alpha(B^{\rm
sym}x_{k}-b),$ (32)

which is AOR for solving

$\displaystyle(\mu I-B+B^{\intercal})x=b$

with $r=\frac{2\alpha\mu}{1+\alpha\mu}$ and
$\omega=\frac{\alpha\mu}{1+\alpha\mu}$. The left hand side of (32) is an
invertible lower-triangular matrix and $x_{k+1}$ can be efficiently computed
by a forward substitution. We can also use the splitting (30) to get another
AOR scheme

$\frac{x_{k+1}-x_{k}}{\alpha}=-\left(\mu x_{k+1}-B^{\rm
sym}x_{k}+2B^{\intercal}x_{k+1}-b\right),$ (33)

for which $x_{k+1}$ can be efficiently computed by a backward substitution.

For a symmetric matrix $M$, we define

$\|x\|_{M}^{2}:=(x,x)_{M}:=x^{\intercal}Mx.$

When $M$ is SPD, $(\cdot,\cdot)_{M}$ defines an inner product with the induced
norm $\|x\|_{M}$. For a general symmetric matrix, $\|\cdot\|_{M}$ may not be a
norm. However the following identity for squares still holds:

$2(a,b)_{M}=\|a\|_{M}^{2}+\|b\|_{M}^{2}-\|a-b\|_{M}^{2}.$ (34)

Without loss of generality, we focus on the convergence analysis of scheme
(31). Notice that when $0<\alpha<1/L_{B^{\rm sym}}$, $I\pm\alpha B^{\rm sym}$
is SPD. For $0<\alpha<1/L_{B^{\rm sym}}$, consider the Lyapunov function

$\mathcal{E}^{\alpha
B}(x)=\frac{1}{2}\|x-x^{\star}\|^{2}-\frac{\alpha}{2}\|x-x^{\star}\|^{2}_{B^{\rm
sym}}=\frac{1}{2}\|x-x^{\star}\|^{2}_{I-\alpha B^{\rm sym}},$ (35)

which is nonnegative and $\mathcal{E}^{\alpha B}(x)=0$ if and only if
$x=x^{\star}$. Noted that $\mathcal{E}^{\alpha B}$ might be negative without
control on the step size as $B^{\rm sym}$ is indefinite.

###### Theorem 3.1

Let $\\{x_{k}\\}$ be the sequence generated by (31) with arbitrary initial
guess $x_{0}$ and step size $0<\alpha<1/L_{B^{\rm sym}}$. Then for the
Lyapunov function (35),

$\displaystyle\mathcal{E}^{\alpha B}(x_{k+1})\leq$
$\displaystyle~{}\frac{1}{1+\alpha\mu}\mathcal{E}^{\alpha B}(x_{k}).$ (36)

In particular, for $\alpha=\displaystyle\frac{1}{2L_{B^{\rm sym}}}$, we have

$\|x_{k}-x^{\star}\|^{2}\leq\left(\frac{1}{1+\mu/(2L_{B^{\rm
sym}})}\right)^{k}3\|x_{0}-x^{\star}\|^{2}.$ (37)

###### Proof

We use the identity for squares (34):

$\frac{1}{2}\|x_{k+1}-x^{\star}\|^{2}-\frac{1}{2}\|x_{k}-x^{\star}\|^{2}=(x_{k+1}-x^{\star},x_{k+1}-x_{k})-\frac{1}{2}\|x_{k+1}-x_{k}\|^{2}.$
(38)

We write GSS scheme (31) as a correction of the implicit Euler scheme

$\displaystyle x_{k+1}-x_{k}$
$\displaystyle=-\alpha(\mathcal{A}(x_{k+1})-\mathcal{A}(x^{\star}))+\alpha
B^{\rm sym}(x_{k+1}-x_{k}).$

For the first term, we have

$-\alpha\langle
x_{k+1}-x^{\star},\mathcal{A}(x_{k+1})-\mathcal{A}(x^{\star})\rangle=-\alpha\mu\|x_{k+1}-x^{\star}\|^{2}.$

We use the identity (34) to expand the cross term as

$\displaystyle(x_{k+1}-x^{\star},x_{k+1}-x_{k})_{\alpha B^{\rm sym}}={}$
$\displaystyle\frac{1}{2}\|x_{k+1}-x^{\star}\|_{\alpha B^{\rm
sym}}^{2}+\frac{1}{2}\|x_{k+1}-x_{k}\|_{\alpha B^{\rm sym}}^{2}$
$\displaystyle-$ $\displaystyle\frac{1}{2}\|x_{k}-x^{\star}\|_{\alpha B^{\rm
sym}}^{2}.$

Substitute back to (38) and rearrange the terms, we obtain the identity

$\displaystyle\mathcal{E}^{\alpha B}(x_{k+1})-\mathcal{E}^{\alpha B}(x_{k})$
$\displaystyle={}$
$\displaystyle-\alpha\mu\|x_{k+1}-x^{\star}\|^{2}-\frac{1}{2}\|x_{k+1}-x_{k}\|_{I-\alpha
B^{\rm sym}}^{2}$ $\displaystyle={}$
$\displaystyle-\alpha\mu\,\mathcal{E}^{\alpha
B}(x_{k+1})-\frac{\alpha\mu}{2}\|x_{k+1}-x^{\star}\|_{I+\alpha B^{\rm
sym}}^{2}-\frac{1}{2}\|x_{k+1}-x_{k}\|_{I-\alpha B^{\rm sym}}^{2}.$

As $0<\alpha<\displaystyle 1/L_{B^{\rm sym}}$, $I\pm\alpha B^{\rm sym}$ is SPD
and the last two terms are non-positive. Dropping them, we obtain the
inequality

$\mathcal{E}^{\alpha B}(x_{k+1})-\mathcal{E}^{\alpha
B}(x_{k})\leq-\alpha\mu\,\mathcal{E}^{\alpha B}(x_{k+1})$

and (36) follows by arrangement.

When $\alpha=1/(2L_{B^{\rm sym}})$, we have the bound

$\frac{1}{4}\|x-x^{\star}\|^{2}\leq\mathcal{E}^{\alpha
B}(x)\leq\frac{3}{4}\|x-x^{\star}\|^{2}$ (39)

which implies (37).

For SOR type iterative methods, usually spectral analysis bai2007successive ;
bai2005generalized ; hadjidimos2000successive  is applied to the error matrix
which is in general non-symmetric and harder to estimate. The Lyapunov
analysis provides a new and relatively simple tool and more importantly
enables us to study nonlinear systems.

###  3.5 Gradient and skew-Symmetric Splitting Methods for Nonlinear Problems

For non-linear equation (1), we treat $\nabla F$ explicitly and propose the
following Gradient and skew-Symmetric Splitting (GSS) scheme

$\frac{x_{k+1}-x_{k}}{\alpha}=-\left(\nabla F(x_{k})+B^{\rm
sym}x_{k}-2Bx_{k+1}\right).$ (40)

Similarly, modification on (33) gives another GSS scheme

$\frac{x_{k+1}-x_{k}}{\alpha}=-\left(\nabla F(x_{k})-B^{\rm
sym}x_{k}+2B^{\intercal}x_{k+1}\right).$ (41)

Both schemes are explicit as $\nabla F(x_{k})$ is known.

We focus on the GSS method (40). The proof for (41) follows in line with a
sign change in the following Lyapunov function

$\mathcal{E}^{\alpha BD}(x_{k}):=\frac{1}{2}\|x_{k}-x^{\star}\|^{2}_{I-\alpha
B^{\rm sym}}-\alpha D_{F}(x^{\star},x_{k}).$ (42)

###### Lemma 2

For $F\in\mathcal{S}_{\mu,L_{F}}$ and $\alpha<1/\max\\{2L_{B^{\rm
sym}},2L_{F}\\}$, the Lyapunov function $\mathcal{E}^{\alpha BD}(x)\geq 0$ and
$\mathcal{E}^{\alpha BD}(x)=0$ if and only if $x=x^{\star}$.

###### Proof

Since $F\in\mathcal{S}_{\mu,L_{F}}$, the Bregman divergence is non-negative
and

$D_{F}(x^{\star},x)\leq\frac{L_{F}}{2}\|x-x^{\star}\|^{2}.$

Then for $\alpha<1/\max\\{2L_{B^{\rm sym}},2L_{F}\\}$,

$\displaystyle\mathcal{E}^{\alpha BD}(x)$
$\displaystyle\geq\frac{1}{2}\|x-x^{\star}\|^{2}-\frac{\alpha}{2}\|x-x^{\star}\|_{B^{\rm
sym}}-\frac{\alpha L_{F}}{2}\|x-x^{\star}\|^{2}$
$\displaystyle\geq\frac{1}{2}\|x-x^{\star}\|^{2}-\frac{1}{4L_{B^{\rm
sym}}}\|x-x^{\star}\|_{B^{\rm sym}}^{2}-\frac{1}{4}\|x-x^{\star}\|^{2}$
$\displaystyle\geq 0.$

If $x\neq x^{\star}$, then the second $\geq$ becomes $>$. So
$\mathcal{E}^{\alpha BD}(x)=0$ if and only if $x=x^{\star}$.

###### Theorem 3.2

Let $\\{x_{k}\\}$ be the sequence generated by the AOR method (40) with
arbitrary initial guess $x_{0}$ and step size $\alpha<1/\max\\{2L_{B^{\rm
sym}},2L_{F}\\}$. Then for the discrete Lyapunov function (42), we have

$\displaystyle\mathcal{E}^{\alpha BD}(x_{k+1})\leq$
$\displaystyle~{}\frac{1}{1+\alpha\mu}\mathcal{E}^{\alpha BD}(x_{k}).$

In particular, for $\alpha=\displaystyle\min\left\\{\frac{1}{4L_{B^{\rm
sym}}},\frac{1}{4L_{F}}\right\\}$, we have

$\|x_{k}-x^{\star}\|^{2}\leq\left(1+1/\max\left\\{4\kappa(B^{\rm
sym}),4\kappa(F)\right\\}\right)^{-k}6\|x_{0}-x^{\star}\|^{2}.$ (43)

###### Proof

We write the scheme (40) as a correction of implicit Euler scheme

$\displaystyle x_{k+1}-x_{k}$
$\displaystyle=-\alpha(\mathcal{A}(x_{k+1})-\mathcal{A}(x^{\star}))+\alpha
B^{\rm sym}(x_{k+1}-x_{k})+\alpha(\nabla F(x_{k+1})-\nabla F(x_{k})).$

The first two terms are treated as before; see the proof of Theorem 3.1. The
last cross term is expanded using the identity (17) for the Bregman
divergence:

$\displaystyle\langle x_{k+1}-x^{\star},\nabla F(x_{x+1})-\nabla
F(x_{k})\rangle=D_{F}(x^{\star},x_{k+1})+D_{F}(x_{k+1},x_{k})-D_{F}(x^{\star},x_{k}).$

Substituting back to (38) and rearranging the terms, we obtain the identity

$\displaystyle\mathcal{E}^{\alpha BD}(x_{k+1})-\mathcal{E}^{\alpha
BD}(x_{k})=$ $\displaystyle~{}-\alpha\mu\mathcal{E}^{\alpha
BD}(x_{k+1})-\frac{\alpha\mu}{2}\|x_{k+1}-x^{\star}\|_{I+\alpha B^{\rm
sym}}^{2}$ $\displaystyle-\alpha^{2}\mu D_{F}(x^{\star},x_{k+1})$
$\displaystyle-\frac{1}{2}\|x_{k+1}-x_{k}\|_{I-\alpha B^{\rm sym}}^{2}+\alpha
D_{F}(x_{k+1},x_{k}).$

When $\alpha<1/\max\\{2L_{B^{\rm sym}},2L_{F}\\}$, $I\pm\alpha B^{\rm sym}$ is
SPD and

$\displaystyle\alpha D_{F}(x_{k+1},x_{k})$ $\displaystyle\leq\frac{\alpha
L_{F}}{2}\|x_{k+1}-x_{k}\|^{2}\leq\frac{1}{4}\|x_{k+1}-x_{k}\|^{2}\leq\frac{1}{2}\|x_{k+1}-x_{k}\|_{I-\alpha
B^{\rm sym}}^{2}.$

Thus the last two terms are non-positive. Dropping all non-positive terms, we
obtain the inequality

$\mathcal{E}^{\alpha BD}(x_{k+1})-\mathcal{E}^{\alpha
BD}(x_{k})\leq-\alpha\mu\mathcal{E}^{\alpha BD}(x_{k+1})$

and the linear reduction follows by arrangement.

When $\alpha=\displaystyle\min\left\\{\frac{1}{4L_{B^{\rm
sym}}},\frac{1}{4L_{F}}\right\\}$, we have the bound

$\frac{1}{8}\|x-x^{\star}\|^{2}\leq\mathcal{E}^{\alpha
BD}(x)\leq\frac{3}{4}\|x-x^{\star}\|^{2}$ (44)

which implies (43).

Notice the rate in Theorem 3.2 is $(1+c/\kappa)$ for
$\kappa=\max\\{\kappa(F),\kappa(B^{\rm sym})\\}$ which matches the rate of the
gradient descent method for convex optimization problems. We expect a
combination of the accelerated gradient flow and AOR for the skew-symmetric
part will give an accelerated explicit scheme.

##  4 Accelerated Gradient and skew-Symmetric Splitting Methods

In this section, we develop the accelerated gradient flow and propose
Accelerated Gradient and skew-Symmetric Splitting (AGSS) methods with
accelerated linear rates. For the implicit-explicit scheme, we can relax to
inexact inner solvers with computable error tolerance.

###  4.1 Accelerated Gradient Flow

Following luoDifferentialEquationSolvers2021a , we introduce an auxiliary
variable $y$ and an accelerated gradient flow

$\left\\{\begin{aligned} x^{\prime}&=y-x,\\\
y^{\prime}&=x-y-\frac{1}{\mu}(\nabla F(x)+\mathcal{N}y).\end{aligned}\right.$
(45)

Denote the vector field on the right hand side of (45) by $\mathcal{G}(x,y)$.
Then $\mathcal{G}(x^{\star},x^{\star})=0$ and thus $(x^{\star},x^{\star})$ is
an equilibrium point of (45). Comparing with the accelerated flow in
luoDifferentialEquationSolvers2021a  for convex optimization, the difference
is to split $\mathcal{A}(x)\rightarrow\nabla F(x)+\mathcal{N}y$. The gradient
component is accelerated by using the accelerated gradient methods for convex
optimization and the skew-symmetric component is accelerated by AOR.

We first show $(x^{\star},x^{\star})$ is exponentially stable. Consider the
Lyapunov function:

$\mathcal{E}(x,y)=D_{F}(x,x^{\star})+\frac{\mu}{2}\|y-x^{\star}\|^{2}.$ (46)

As $F$ is strongly convex, $\mathcal{E}(x,y)\geq 0$ and $\mathcal{E}(x,y)=0$
iff $x=y=x^{\star}$. Furthermore, function
$D_{F}(\cdot,x^{\star})\in\mathcal{S}_{\mu,L_{F}}$.

###### Lemma 3

For function $F\in\mathcal{S}_{\mu}$, we have

$\langle\nabla F(x)-\nabla F(x^{\star}),x-x^{\star}\rangle\geq
D_{F}(x,x^{\star})+\frac{\mu}{2}\|x-x^{\star}\|^{2},$ (47)

###### Proof

By direct computation, $\langle\nabla F(x)-\nabla
F(x^{\star}),x-x^{\star}\rangle=D_{F}(x,x^{\star})+D_{F}(x^{\star},x)$ and
thus (47) follows from the bound

$\min\\{D_{F}(x,x^{\star}),D_{F}(x^{\star},x)\\}\geq\frac{\mu}{2}\|x-x^{\star}\|^{2},$

for a convex function $F\in\mathcal{S}_{\mu}$.

We then verify the strong Lyapunov property.

######  Lemma 4 (Strong Lyapunov Property)

Assume function $F\in\mathcal{S}_{\mu}$. Then for the Lyapunov function (46)
and the accelerated gradient flow vector field $\mathcal{G}$, the following
strong Lyapunov property holds

$-\nabla\mathcal{E}(x,y)\cdot\mathcal{G}(x,y)\geq\mathcal{E}(x,y)+\frac{\mu}{2}\|y-x\|^{2}.$
(48)

###### Proof

First of all, as $\mathcal{G}(x^{\star},x^{\star})=0$,

$-\nabla\mathcal{E}(x,y)\cdot\mathcal{G}(x,y)=-\nabla\mathcal{E}(x,y)\cdot(\mathcal{G}(x,y)-\mathcal{G}(x^{\star},x^{\star})).$

Direct computation gives

$\displaystyle-\nabla\mathcal{E}(x,y)\cdot\mathcal{G}(x,y)=$
$\displaystyle~{}\langle\nabla
D_{F}(x;x^{\star}),x-x^{\star}-(y-x^{\star})\rangle-\mu(y-x^{\star},x-x^{\star})$
$\displaystyle+\mu\|y-x^{\star}\|^{2}+\langle\nabla F(x)-\nabla
F(x^{\star}),y-x^{\star}\rangle+(y-x^{\star},\mathcal{N}(y-x^{\star}))$
$\displaystyle=$ $\displaystyle~{}\langle\nabla F(x)-\nabla
F(x^{\star}),x-x^{\star}\rangle+\mu\|y-x^{\star}\|^{2}-\mu(y-x^{\star},x-x^{\star})$

where we have used $\nabla D_{F}(x;x^{\star})=\nabla F(x)-\nabla F(x^{\star})$
and $(y-x^{\star},\mathcal{N}(y-x^{\star}))=0$ since $\mathcal{N}$ is skew-
symmetric.

Using the bound (47) and the identity for squares (34):

$\displaystyle\frac{\mu}{2}\|y-x^{\star}\|^{2}-\mu(y-x^{\star},x-x^{\star})=\frac{\mu}{2}\|y-x\|^{2}-\frac{\mu}{2}\|x-x^{\star}\|^{2},$

we obtain (48).

The calculation is more clear when $\nabla F(x)=Hx$ is linear with
$H=\nabla^{2}F\geq\mu I$. Then $-\nabla\mathcal{E}(x,y)\cdot\mathcal{G}(x,y)$
is a quadratic form of $(x-x^{\star},y-x^{\star})^{\intercal}$. We calculate
the corresponding matrix as

$\begin{pmatrix}H&0\\\ 0&\mu I\end{pmatrix}\begin{pmatrix}I&-I\\\
-I+\frac{1}{\mu}H&I+\frac{1}{\mu}\mathcal{N}\end{pmatrix}=\begin{pmatrix}H&-H\\\
-\mu I+H&\mu+\mathcal{N}\end{pmatrix}.$

For a quadratic form, $v^{\intercal}Mv=v^{\intercal}\operatorname{sym}(M)v$.
We have

$\displaystyle\operatorname{sym}\begin{pmatrix}H&-H\\\ -\mu I+H&\mu
I+\mathcal{N}\end{pmatrix}$ $\displaystyle=\begin{pmatrix}H&-\mu I/2\\\ -\mu
I/2&\mu I\end{pmatrix}$ $\displaystyle\geq\begin{pmatrix}H/2&0\\\ 0&\mu
I/2\end{pmatrix}+\frac{\mu}{2}\begin{pmatrix}I&-I\\\ -I&I\end{pmatrix},$

where in the last step we use the convexity $H\geq\mu I$. Then (48) follows.

To see how the condition number changes using the accelerated gradient flow,
we consider the following $2\times 2$ matrix

$G=\begin{pmatrix}-1&1\\\ 1-a&\quad-1+b{\rm i}\end{pmatrix}$ (49)

with $a\geq 1$ representing the eigenvalue of $\nabla^{2}F/\mu$ and $b{\rm i}$
for $\mathcal{N}/\mu$ as the eigenvalue of skew-symmetric matrix is pure
imaginary. Then the eigenvalue of $G$ is

$\lambda(G)=-1+\frac{b\pm\sqrt{b^{2}+4(a-1)}}{2}{\rm i}.$

The real part is always $-1$ which implies the decay property of ODE
$x^{\prime}=Gx$. The spectral radius is

$|\lambda|=\mathcal{O}(\sqrt{a})+\mathcal{O}(|b|),\quad\text{as }a\gg 1,|b|\gg
1.$

As a comparison, $|a+b{\rm
i}|=\sqrt{a^{2}+b^{2}}=\mathcal{O}(a)+\mathcal{O}(|b|)$. We accelerate the
dependence from $\mathcal{O}(a)$ to $\mathcal{O}(\sqrt{a})$.

###  4.2 Implicit Euler Schemes

If we apply the implicit Euler method for the accelerated gradient system
(45), the linear convergence is a direct consequence of the strong Lyapunov
property (48). Consider

$\displaystyle\frac{x_{k+1}-x_{k}}{\alpha_{k}}$
$\displaystyle=\mathcal{G}^{x}(x_{k+1},y_{k+1}):=y_{k+1}-x_{k+1},$ (50a)
$\displaystyle\frac{y_{k+1}-y_{k}}{\alpha_{k}}$
$\displaystyle=\mathcal{G}^{y}(x_{k+1},y_{k+1}):=x_{k+1}-y_{k+1}-\frac{1}{\mu}(\nabla
F(x_{k+1})+\mathcal{N}y_{k+1}).$ (50b)

As we have shown before, the implicit Euler method is unconditionally stable
and yield superlinear convergence with the price of solving a nonlinear
equation system where $x_{k+1}$ and $y_{k+1}$ are coupled together, which is
usually not practical. We present the convergence analysis here as a reference
to measure the difference of other methods.

Denote $\mathcal{E}_{k}=\mathcal{E}(x_{k},y_{k})$ and $z_{k}=(x_{k},y_{k})$.
Using the $\mu$-convexity of the Lyapunov function (46) and the strong
Lyapunov property in Lemma 4, we have

$\displaystyle\mathcal{E}_{k+1}-\mathcal{E}_{k}\leq{}$
$\displaystyle(\nabla\mathcal{E}_{k+1},z_{k+1}-z_{k})-\frac{\mu}{2}\|z_{k+1}-z_{k}\|^{2}$
$\displaystyle={}$
$\displaystyle\alpha_{k}\left\langle{\mathcal{\nabla}\mathcal{E}_{k+1},\mathcal{G}(z_{k+1})}\right\rangle-\frac{\mu}{2}\|z_{k+1}-z_{k}\|^{2}$
$\displaystyle\leq{}$ $\displaystyle-\alpha_{k}\mathcal{E}_{k+1},$

from which the linear convergence follows naturally

$\mathcal{E}_{k+1}\leq\frac{1}{1+\alpha_{k}}\mathcal{E}_{k}$

for arbitrary step size $\alpha_{k}>0$.

###  4.3 Implicit and Explicit Schemes

If we treat the skew symmetric part implicit, we can achieve the acceleration
like the convex optimization problem. Consider the following IMplicit and
EXplicit (IMEX) scheme of the accelerated gradient flow:

$\displaystyle\frac{\hat{x}_{k+1}-x_{k}}{\alpha_{k}}$
$\displaystyle=y_{k}-\hat{x}_{k+1},$ (51a)
$\displaystyle\frac{y_{k+1}-y_{k}}{\alpha_{k}}$
$\displaystyle=\hat{x}_{k+1}-y_{k+1}-\frac{1}{\mu}\left(\nabla
F(\hat{x}_{k+1})+\mathcal{N}y_{k+1}\right),$ (51b)
$\displaystyle\frac{x_{k+1}-x_{k}}{\alpha_{k}}$
$\displaystyle=y_{k+1}-x_{k+1}.$ (51c)

We first treat $y$ is known as $y_{k}$ and solve the first equation to get
$\hat{x}_{k+1}$ and then with known $\hat{x}_{k+1}$ to solve the following
shifted skew-symmetric equation

$\left[(1+\alpha_{k})I+\frac{\alpha_{k}}{\mu}\mathcal{N}\right]y_{k+1}=b(\hat{x}_{k+1},y_{k}),$
(52)

with known right hand side
$b(\hat{x}_{k+1},y_{k})=\alpha_{k}\hat{x}_{k+1}-\frac{\alpha_{k}}{\mu}\nabla
F(\hat{x}_{k+1})+y_{k}$. Then with computed $y_{k+1}$, we can solve for
$x_{k+1}$ again using an implicit discretization of $x^{\prime}=y-x$. In terms
of the ODE solvers, (51) is known as the predictor-corrector method. The
intermediate approximation $\hat{x}_{k+1}$ is a predictor and $x_{k+1}$ is a
corrector.

As the skew-symmetric part is treated implicitly, the scheme is expected to
achieve the accelerated linear rate as the accelerated gradient method for
minimizing the convex function $F$. For simplicity, we denote

$\mathcal{E}_{k}=\mathcal{E}(x_{k},y_{k}),\quad\hat{\mathcal{E}}_{k+1}=\mathcal{E}(\hat{x}_{k+1},y_{k+1})$
(53)

for the Lyapunov function $\mathcal{E}$ defined in (46).

###### Lemma 5

Assume function $F\in\mathcal{S}_{\mu}$. Let $(\hat{x}_{k},y_{k})$ be the
sequence generated by the accelerated gradient method (51). Then for $k\geq
0$,

$\displaystyle\hat{\mathcal{E}}_{k+1}-\mathcal{E}_{k}\leq$
$\displaystyle~{}-\alpha_{k}\hat{\mathcal{E}}_{k+1}-\alpha_{k}\left\langle\nabla
D_{F}(\hat{x}_{k+1},x^{\star}),y_{k+1}-y_{k}\right\rangle-\frac{\mu}{2}\left\|y_{k+1}-y_{k}\right\|^{2}$
(54)
$\displaystyle-\frac{\alpha_{k}\mu}{2}\|y_{k+1}-\hat{x}_{k+1}\|^{2}-\frac{\mu}{2}\|\hat{x}_{k+1}-x_{k}\|^{2}.$

###### Proof

The Lyapunov function $\mathcal{E}(x,y)$ is $\mu$-convex. Thus we have

$\displaystyle\hat{\mathcal{E}}_{k+1}-\mathcal{E}_{k}\leq{}$
$\displaystyle\langle\partial_{x}\mathcal{E}(\hat{x}_{k+1},y_{k+1}),\hat{x}_{k+1}-x_{k}\rangle+\langle\partial_{y}\mathcal{E}(\hat{x}_{k+1},y_{k+1}),y_{k+1}-y_{k}\rangle$
$\displaystyle-\frac{\mu}{2}\|\hat{x}_{k+1}-x_{k}\|^{2}-\frac{\mu}{2}\left\|y_{k+1}-y_{k}\right\|^{2}.$

Then substitute

$\displaystyle\hat{x}_{k+1}-x_{k}$
$\displaystyle=\alpha_{k}\mathcal{G}^{x}(\hat{x}_{k+1},y_{k+1})+\alpha_{k}(y_{k}-y_{k+1})$
$\displaystyle y_{k+1}-y_{k}$
$\displaystyle=\alpha_{k}\mathcal{G}^{y}(\hat{x}_{k+1},y_{k+1}),$

and use the strong Lyapunov property (48) at $(\hat{x}_{k+1},y_{k+1})$ to get
the desired result.

The term $-\alpha_{k}\left\langle\nabla
D_{F}(\hat{x}_{k+1};x^{\star}),y_{k+1}-y_{k}\right\rangle$ on the right hand
side of (54) accounts for
$-\alpha_{k}\langle\partial_{x}\mathcal{E}(\hat{x}_{k+1},y_{k+1}),y_{k+1}-y_{k}\rangle$
for using explicit $y_{k}$ in (51a), which is again a mis-match term compare
with $\mathcal{G}(\hat{x}_{k+1},y_{k+1})$ used in the implicit Euler scheme
(50a). The correction step (51c) will be used to bound the mis-match term. We
restate and generalize the result in luoDifferentialEquationSolvers2021a  in
the following lemma.

###### Lemma 6

Assume $\nabla F$ is $L_{F}$-Lipschitz continuous and
$(\hat{x}_{k+1},y_{k+1})$ is generated from $(x_{k},y_{k})$ such that: there
exists $c_{1},c_{2},c_{3}>0$ satisfying

$\hat{\mathcal{E}}_{k+1}-\mathcal{E}_{k}\leq~{}-\alpha_{k}c_{1}\hat{\mathcal{E}}_{k+1}-\alpha_{k}c_{2}\left\langle\nabla
D_{F}(\hat{x}_{k+1},x^{\star}),y_{k+1}-y_{k}\right\rangle-\frac{c_{3}}{2}\left\|y_{k+1}-y_{k}\right\|^{2}.$
(55)

Then set $x_{k+1}$ by the relation

$(1+\alpha_{k}c_{1})(x_{k+1}-\hat{x}_{k+1})=\alpha_{k}c_{2}(y_{k+1}-y_{k}).$
(56)

and choose $\alpha_{k}>0$ satisfying

$\alpha_{k}^{2}L_{F}c_{2}^{2}\leq(1+\alpha_{k}c_{1})c_{3},$ (57)

we have the linear convergence

$\mathcal{E}_{k+1}\leq\frac{1}{1+\alpha_{k}c_{1}}\mathcal{E}_{k}.$

###### Proof

Using the relation (56), we rewrite (55) into

$\displaystyle(1+\alpha_{k}c_{1})\hat{\mathcal{E}}_{k+1}-\mathcal{E}_{k}$ (58)
$\displaystyle\leq$ $\displaystyle~{}-(1+c_{1}\alpha_{k})\left\langle\nabla
D_{F}(\hat{x}_{k+1},x^{\star}),x_{k+1}-\hat{x}_{k+1}\right\rangle-\frac{(1+\alpha_{k}c_{1})^{2}c_{3}}{2\alpha_{k}^{2}c_{2}^{2}}\left\|x_{k+1}-\hat{x}_{k+1}\right\|^{2}.$

As $D_{F}(\cdot,x^{\star})\in\mathcal{S}_{\mu,L_{F}}$, we have

$\displaystyle\mathcal{E}_{k+1}-\hat{\mathcal{E}}_{k+1}$
$\displaystyle=D_{F}(x_{k+1},x^{\star})-D_{F}(\hat{x}_{k+1},x^{\star})$ (59)
$\displaystyle\leq\left\langle\nabla
D_{F}(\hat{x}_{k+1},x^{\star}),x_{k+1}-\hat{x}_{k+1}\right\rangle+\frac{L_{F}}{2}\|x_{k+1}-\hat{x}_{k+1}\|^{2}.$

Summing $(1+c_{1}\alpha_{k})\times$(59) and (58) we get

$\displaystyle(1+c_{1}\alpha_{k})\mathcal{E}_{k+1}-\mathcal{E}_{k}\leq-\left(\frac{(1+\alpha_{k}c_{1})^{2}c_{3}}{2\alpha_{k}^{2}c_{2}^{2}}-\frac{(1+\alpha_{k}c_{1})L_{F}}{2}\right)\left\|x_{k+1}-\hat{x}_{k+1}\right\|^{2}\leq
0$

according to our choice of $\alpha_{k}$. Rearrange the terms to get the
desired inequality.

###### Remark 2

The largest step size $\alpha_{k}$ satisfying (57) is given by

$\alpha_{k}=\frac{c_{1}c_{3}+\sqrt{c_{1}^{2}c_{3}^{2}+4L_{F}c_{2}^{2}c_{3}}}{2L_{F}c_{2}^{2}}\geq\frac{1}{c_{2}}\sqrt{\frac{c_{3}}{L_{F}}},$

where the last one is a simplified formula for the step size satisfying (57).

We showed the decay of Lyapunov function in Lemma 5 satisfying the assumption
of Lemma 6 with $c_{1}=c_{2}=1,c_{3}=\mu$ and the correction step (51c)
matches the relation (56). As a result, we state the following linear
convergence rate of the IMEX scheme (51).

###### Theorem 4.1

Assume function $F\in\mathcal{S}_{\mu,L_{F}}$. Let $(x_{k},y_{k})$ be the
sequence generated by the accelerated gradient method (51) with arbitrary
initial value and step size satisfying

$\alpha_{k}^{2}L_{F}\leq(1+\alpha_{k})\mu,$

then for the Lyapunov function (46),

$\mathcal{E}_{k+1}\leq\frac{1}{1+\alpha_{k}}\mathcal{E}_{k}.$

For $\alpha_{k}=\sqrt{\mu/L_{F}}=1/\sqrt{\kappa(F)}$, we achieve the
accelerated rate

$\mathcal{E}_{k}\leq\left(\frac{1}{1+1/\sqrt{\kappa(F)}}\right)^{k}\mathcal{E}_{0},$

which implies

$\|x_{k}-x^{\star}\|^{2}+\|y_{k}-x^{\star}\|^{2}\leq\left(\frac{1}{1+1/\sqrt{\kappa(F)}}\right)^{k}2\mathcal{E}_{0}/\mu.$

For linear problems, comparing with HSS scheme (6), we achieve the same
accelerated rate for linear and nonlinear systems without computing the
inverse of the symmetric part $(\beta I+\nabla^{2}F)^{-1}$.

###  4.4 Inexact Solvers for the Shifted skew-Symmetric System

In practice, we can relax the inner solver to be an inexact approximation.
That is we solve the equation (52) up to a residual $\varepsilon_{\rm
in}=b-\left[(1+\alpha_{k})I+\frac{\alpha_{k}}{\mu}\mathcal{N}\right]y_{k+1}$.
The scheme can be modified to

$\displaystyle\frac{\hat{x}_{k+1}-x_{k}}{\alpha_{k}}$
$\displaystyle=y_{k}-\hat{x}_{k+1},$ (60a)
$\displaystyle\frac{y_{k+1}-y_{k}}{\alpha_{k}}$
$\displaystyle=\hat{x}_{k+1}-y_{k+1}-\frac{1}{\mu}\left(\nabla
F(\hat{x}_{k+1})+\mathcal{N}y_{k+1}\right)-\frac{\varepsilon_{\rm
in}}{\alpha_{k}},$ (60b) $\displaystyle\frac{x_{k+1}-x_{k}}{\alpha_{k}}$
$\displaystyle=y_{k+1}-\frac{1}{2}(x_{k+1}+\hat{x}_{k+1}).$ (60c)

Notice that in the third step (60c), we use
$\frac{1}{2}(x_{k+1}+\hat{x}_{k+1})$ instead of $x_{k+1}$ for discretization
variable $x$ at $k+1$. The perturbation $\varepsilon_{\rm in}$ is the residual
of the linear equation (52).

###### Corollary 1

If we compute $y_{k+1}$ such that the residual of (60b) satisfies

$\|\varepsilon_{\rm
in}\|^{2}\leq\frac{\alpha_{k}}{2}\left(\|\hat{x}_{k+1}-x_{k}\|^{2}+\alpha_{k}\|y_{k+1}-\hat{x}_{k+1}\|^{2}\right),$
(61)

then for $(x_{k},y_{k})$ be the sequence generated by the inexact accelerated
gradient method (60) with arbitrary initial value and step size satisfying

$\alpha_{k}^{2}L_{F}\leq(1+\alpha_{k}/2)\mu,$

we have the linear convergence

$\mathcal{E}_{k+1}\leq\frac{1}{1+\alpha_{k}/2}\mathcal{E}_{k}.$

For $\alpha_{k}=\sqrt{\mu/L_{F}}=1/\sqrt{\kappa(F)}$, we achieve the
accelerated rate

$\mathcal{E}_{k}\leq\left(\frac{1}{1+1/(2\sqrt{\kappa(F)})}\right)^{k}\mathcal{E}_{0}.$

###### Proof

We write (60b) as

$y_{k+1}-y_{k}=\alpha_{k}\mathcal{G}^{y}(\hat{x}_{k+1},y_{k+1})+\varepsilon_{\rm
in}.$

Compared with Lemma 5, the inexactness introduces a mis-match term
$\mu(y_{k+1}-x^{\star},\epsilon_{\rm out})$ in
$\langle\partial_{y}\mathcal{E}(\hat{x}_{k+1},y_{k+1}),y_{k+1}-y_{k}\rangle$
which can be bounded by

$\displaystyle\mu|(y_{k+1}-x^{\star},\epsilon_{\rm out})|$
$\displaystyle\leq\frac{\alpha_{k}\mu}{4}\|y_{k+1}-x^{\star}\|^{2}+\frac{\mu}{\alpha_{k}}\|\varepsilon_{\rm
in}\|^{2}$
$\displaystyle\leq\frac{\alpha_{k}\mu}{4}\|y_{k+1}-x^{\star}\|^{2}+\frac{\mu}{2}\left(\|\hat{x}_{k+1}-x_{k}\|^{2}+\alpha_{k}\|y_{k+1}-\hat{x}_{k+1}\|^{2}\right).$

Use $-\alpha_{k}\hat{\mathcal{E}}_{k+1}/2$ to cancel the first term and the
additional quadratic term in (54) to cancel the second. Then we have

$\displaystyle\hat{\mathcal{E}}_{k+1}-\mathcal{E}_{k}\leq-\frac{\alpha_{k}}{2}\hat{\mathcal{E}}_{k+1}-\alpha_{k}\left\langle\nabla
D_{F}(\hat{x}_{k+1},x^{\star}),y_{k+1}-y_{k}\right\rangle-\frac{\mu}{2}\left\|y_{k+1}-y_{k}\right\|^{2},$
(62)

which with the correction step (60c) satisfying assumptions in Lemma 6 with
$c_{1}=1/2,c_{2}=1$ and $c_{3}=\mu$.

Solvers for the shifted skew-symmetric equation (52) with the form $(\beta
I+\mathcal{N})y=b(\hat{x}_{k+1},y_{k})$ is discussed in Section 3. In
particular, $\beta=(1+\sqrt{\mu/L_{F}})\sqrt{\mu L_{F}}$ if we choose
$\alpha_{k}=\sqrt{\mu/L_{F}}$. The inner iteration steps are roughly $C_{\rm
in}\sqrt{\kappa(F)}L_{B^{\rm sym}}/L_{F}.$ The outer iteration is $C_{\rm
out}\sqrt{\kappa(F)}$. Constants $C_{\rm in}=\mathcal{O}(|\ln\varepsilon_{\rm
in}|)$ and $C_{\rm out}=\mathcal{O}(|\ln\epsilon_{\rm out}|)$ depend on the
tolerance $\varepsilon_{\rm in}$ for the inner iteration and $\epsilon_{\rm
out}$ for the outer iteration. Therefore

  * •

$\nabla F(x_{k})$ gradient evaluation: $C_{\rm out}\sqrt{\kappa(F)}$;

  * •

$(\beta I+B)^{-1}b$ matrix-vector multiplication: $C_{\rm in}C_{\rm
out}\kappa(B^{\rm sym})$,

where we use the relation $\kappa(F)L_{B^{\rm sym}}/L_{F}=\kappa(B^{\rm
sym})$.

###  4.5 Accelerated Gradient and skew-Symmetric Splitting Methods

Combining the MEX scheme in Section 4.3 and accelerated overrelexation
technique in Section 3.5, we propose the following explicit discretization of
the accelerated gradient flow:

$\displaystyle\frac{\hat{x}_{k+1}-x_{k}}{\alpha}$
$\displaystyle=y_{k}-\hat{x}_{k+1},$ (63a)
$\displaystyle\frac{y_{k+1}-y_{k}}{\alpha}$
$\displaystyle=\hat{x}_{k+1}-y_{k+1}-\frac{1}{\mu}\left(\nabla
F(\hat{x}_{k+1})+B^{\rm sym}y_{k}-2By_{k+1}\right),$ (63b)
$\displaystyle\frac{x_{k+1}-x_{k}}{\alpha}$
$\displaystyle=y_{k+1}-\frac{1}{2}(x_{k+1}+\hat{x}_{k+1}).$ (63c)

The update of $y_{k+1}$ is equivalent to solve a lower triangular linear
algebraic system

$\left((1+\alpha)I+\frac{2\alpha}{\mu}B\right)y_{k+1}=b(\hat{x}_{k+1},y_{k})$

with
$b(\hat{x}_{k+1},y_{k})=y_{k}+\alpha\hat{x}_{k+1}-\frac{\alpha}{\mu}\nabla
F(\hat{x}_{k+1})-\frac{\alpha}{\mu}(B^{\intercal}+B)y_{k}$ which can be
computed efficiently by a forward substitution and thus no inexact inner
solver is needed. Subtracting (63c) from (63a) implies the relation (56) with
$c_{1}=1/2,c_{2}=1$.

Consider the modified Lyapunov function:

$\displaystyle\mathcal{E}^{\alpha B}(x,y)$
$\displaystyle=D_{F}(x,x^{\star})+\frac{1}{2}\|y-x^{\star}\|^{2}_{\mu I-\alpha
B^{\rm sym}}.$ (64)

For $\displaystyle 0<\alpha<\frac{\mu}{L_{B^{\rm sym}}}$, we have $\mu
I-\alpha B^{\rm sym}$ is positive definite. Then $\mathcal{E}\geq 0$ and
$\mathcal{E}=0$ if and only if $x=y=x^{\star}$. We denote

$\mathcal{E}^{\alpha B}_{k}=\mathcal{E}^{\alpha
B}(x_{k},y_{k}),\quad\hat{\mathcal{E}}^{\alpha B}_{k}=\mathcal{E}^{\alpha
B}(\hat{x_{k}},y_{k}).$

###### Lemma 7

Assume function $F\in\mathcal{S}_{\mu,L_{F}}$. Let $(\hat{x}_{k},y_{k})$ be
the sequence generated by the accelerated gradient method (63) with
$0<\alpha\leq\displaystyle\frac{\mu}{2L_{B^{\rm sym}}}$. Then, for $k\geq 0$,
the modified Lyapunov function (64) satisfies

$\displaystyle\hat{\mathcal{E}}^{\alpha B}_{k+1}-\mathcal{E}^{\alpha
B}_{k}\leq$ $\displaystyle~{}-\frac{\alpha}{2}\hat{\mathcal{E}}^{\alpha
B}_{k+1}-\alpha\left\langle\nabla
D_{F}(\hat{x}_{k+1},x^{\star}),y_{k+1}-y_{k}\right\rangle-\frac{\mu}{4}\left\|y_{k+1}-y_{k}\right\|^{2}.$
(65)

###### Proof

We write the difference of the modified Lyapunov function into

$\displaystyle\hat{\mathcal{E}}^{\alpha B}_{k+1}-\mathcal{E}^{\alpha B}_{k}=$
$\displaystyle~{}\hat{\mathcal{E}}_{k+1}-\mathcal{E}_{k}-\frac{1}{2}\|y_{k+1}-x^{\star}\|^{2}_{\alpha
B^{\rm sym}}+\frac{1}{2}\|y_{k}-x^{\star}\|^{2}_{\alpha B^{\rm sym}}.$ (66)

where $\hat{\mathcal{E}}_{k+1},\mathcal{E}_{k}$ are defined in (53) and
$\mathcal{E}$ refers to the Lyapunov function (46). Since $\mathcal{E}$ is
$\mu$-convex,

$\displaystyle\hat{\mathcal{E}}_{k+1}-\mathcal{E}_{k}\leq{}$
$\displaystyle\langle\partial_{x}\mathcal{E}(\hat{x}_{k+1},y_{k+1}),\hat{x}_{k+1}-x_{k}\rangle-\frac{\mu}{2}\|\hat{x}_{k+1}-x_{k}\|^{2}$
$\displaystyle~{}\langle\partial_{y}\mathcal{E}(\hat{x}_{k+1},y_{k+1}),y_{k+1}-y_{k}\rangle-\frac{\mu}{2}\|y_{k+1}-y_{k}\|^{2}.$

Then write the scheme as a correction of the implicit Euler scheme

$\displaystyle\hat{x}_{k+1}-x_{k}$
$\displaystyle=\alpha\mathcal{G}^{x}(\hat{x}_{k+1},y_{k+1})+\alpha(y_{k}-y_{k+1}),$
$\displaystyle y_{k+1}-y_{k}$
$\displaystyle=\alpha\mathcal{G}^{y}(\hat{x}_{k+1},y_{k+1})+\frac{\alpha}{\mu}B^{\rm
sym}(y_{k+1}-y_{k}).$

According to the proof of Lemma 5, we get

$\displaystyle\hat{\mathcal{E}}_{k+1}-\mathcal{E}_{k}\leq$
$\displaystyle~{}-\alpha\hat{\mathcal{E}}_{k+1}-\alpha\left\langle\nabla
D_{F}(\hat{x}_{k+1},x^{\star}),y_{k+1}-y_{k}\right\rangle-\frac{\mu}{2}\left\|y_{k+1}-y_{k}\right\|^{2}$
(67) $\displaystyle+(y_{k+1}-x^{\star},y_{k+1}-y_{k})_{\alpha B^{\rm sym}}.$

We use the identity (34) to expand the last cross term in (67) as

$\displaystyle(y_{k+1}-x^{\star},y_{k+1}-y_{k})_{\alpha B^{\rm sym}}$
$\displaystyle=\frac{1}{2}\|y_{k+1}-x^{\star}\|_{\alpha B^{\rm
sym}}^{2}+\frac{1}{2}\|y_{k+1}-y_{k}\|_{\alpha B^{\rm
sym}}^{2}-\frac{1}{2}\|y_{k}-x^{\star}\|_{\alpha B^{\rm sym}}^{2}.$

Substitute back to (67) and rearrange terms using (66),

$\displaystyle\hat{\mathcal{E}}^{\alpha B}_{k+1}-\mathcal{E}^{\alpha
B}_{k}\leq$
$\displaystyle~{}-\alpha\hat{\mathcal{E}}_{k+1}-\alpha\left\langle\nabla
D_{F}(\hat{x}_{k+1},x^{\star}),y_{k+1}-y_{k}\right\rangle-\frac{1}{2}\left\|y_{k+1}-y_{k}\right\|^{2}_{\mu
I-\alpha B^{\rm sym}}$ $\displaystyle=$
$\displaystyle~{}-\frac{\alpha}{2}\hat{\mathcal{E}}^{\alpha
B}_{k+1}-\frac{\alpha}{2}D_{F}(\hat{x}_{k+1},x^{\star})-\frac{\alpha}{4}\|y_{k+1}-x^{\star}\|_{\mu
I+\alpha B^{\rm sym}}$ $\displaystyle-\alpha\left\langle\nabla
D_{F}(\hat{x}_{k+1},x^{\star}),y_{k+1}-y_{k}\right\rangle-\frac{1}{2}\left\|y_{k+1}-y_{k}\right\|^{2}_{\mu
I-\alpha B^{\rm sym}}$ $\displaystyle=$
$\displaystyle~{}-\frac{\alpha}{2}\hat{\mathcal{E}}^{\alpha
B}_{k+1}-\alpha\left\langle\nabla
D_{F}(\hat{x}_{k+1},x^{\star}),y_{k+1}-y_{k}\right\rangle-\frac{\mu}{4}\left\|y_{k+1}-y_{k}\right\|^{2}$
$\displaystyle-\frac{\alpha}{2}D_{F}(\hat{x}_{k+1},x^{\star})-\frac{\alpha}{4}\|y_{k+1}-x^{\star}\|_{\mu
I+\alpha B^{\rm sym}}-\frac{1}{4}\left\|y_{k+1}-y_{k}\right\|^{2}_{\mu
I-2\alpha B^{\rm sym}}.$

According to our choice $\alpha\leq\displaystyle\frac{\mu}{2L_{B^{\rm sym}}}$,
both $\mu I-2\alpha B^{\rm sym}$ and $\mu I+\alpha B^{\rm sym}$ are SPD.
Dropping the last three non-positive terms we have the desired result.

The decay of the modified Lyapunov function (65) and the appropriate relation
(56) satisfies assumptions of Lemma 6 with $c_{1}=1/2,c_{2}=1$ and
$c_{3}=\mu/2$. Although the Lyapunov function is slightly different, (59)
still holds for $\mathcal{E}^{\alpha B}$. We conclude with the following
linear convergence rate.

###### Theorem 4.2

Assume function $F\in\mathcal{S}_{\mu,L_{F}}$. Let $(x_{k},y_{k})$ be the
sequence generated by the accelerated gradient method (63) with arbitrary
initial value and step size satisfying

$0<\alpha\leq\min\left\\{\frac{\mu}{2L_{B^{\rm
sym}}},\sqrt{\frac{\mu}{2L_{F}}}\right\\},$

then for the modified Lyapunov function (46),

$\mathcal{E}^{\alpha B}_{k+1}\leq\frac{1}{1+\alpha/2}\mathcal{E}^{\alpha
B}_{k}.$

In particular, $\displaystyle\alpha=\min\left\\{\frac{\mu}{2L_{B^{\rm
sym}}},\sqrt{\frac{\mu}{2L_{F}}}\right\\}$, we achieve the accelerated rate

$\|x_{k+1}-x^{\star}\|^{2}\leq\frac{2}{\mu}\mathcal{E}^{\alpha
B}_{k}\leq\left(1+1/\max\left\\{4\kappa(B^{\rm
sym}),\sqrt{8\kappa(F)}\right\\}\right)^{-k}\frac{2\mathcal{E}_{0}^{\alpha
B}}{\mu}.$

As expected, we have developed an explicit scheme (63) which achieves the
accelerated rate. Therefore the cost of gradient $\nabla F(x_{k})$ evaluation
and matrix-vector $(\beta I+B)^{-1}b$ multiplication are both $C_{\rm
out}\max\\{\kappa(B^{\rm sym}),\sqrt{\kappa(F)}\\}$. Compare with the inexact
IMEX scheme (51), we may need more gradient evaluation but less matrix-vector
multiplication. If $\kappa(B^{\rm sym})\gg\sqrt{\kappa(F)}$ and $\nabla
F(x_{k})$ is computationally expensive to evaluate, then the inexact IMEX
scheme (51) or its inexact version (60) is favorable; otherwise if the error
tolerance for the inexact inner solve is small, i.e., $C_{\rm in}$ is large,
and the cost of $(\beta I+B)^{-1}b$ is comparable to the evaluation of $\nabla
F(x_{k})$, then the explicit scheme (63) takes advantage.

##  5 Nonlinear Saddle Point Systems

In this section, we investigate an important class of optimization problems as
an example of the considered monotone operator equation (1). We derive optimal
algorithms for strongly-convex-strongly-concave saddle point systems. Our
algorithms can be relaxed to achieve accelerated linear convergence for
convex-concave saddle point systems and inexact solvers.

###  5.1 Problem Setting

Consider the nonlinear smooth saddle point system with bilinear coupling:

$\min_{u\in\mathcal{V}}\max_{p\in\mathcal{Q}}\mathcal{L}(u,p)=f(u)-g(p)+(Bu,p)$
(68)

where $\mathcal{V}=\mathbb{R}^{m},\mathcal{Q}=\mathbb{R}^{n},m\geq n$ are
finite-dimensional Hilbert spaces with inner product induced by SPD operators
$\mathcal{I}_{\mathcal{V}},\mathcal{I}_{\mathcal{Q}}$, respectively. Functions
$f(u)\in\mathcal{S}_{\mu_{f},L_{f}},$ and
$g(p)\in\mathcal{S}_{\mu_{g},L_{g}}$. The operator $B$ is an $n\times m$
matrix of full rank.

Convex optimization problems with affine equality constraints can be rewritten
into a saddle point system (68):

$\displaystyle\min_{u\in\mathbb{R}^{m}}f(u)$ (69) $\displaystyle\text{subject
to}\quad Bu=b.$

Then $p$ is the Lagrange multiplier to impose the constraint $Bu=b$ and
$\mathcal{L}(u,p)=f(u)-(b,p)+(Bu,p)$. Notice that $g(p)=(b,p)$ is linear and
not strongly convex, i.e., $\mu_{g}=0$.

The point $(u^{\star},p^{\star})$ solves the min-max problem (68) is said to
be a saddle point of $\mathcal{L}(u,p)$, that is

$\mathcal{L}(u^{\star},p)\leq\mathcal{L}(u^{\star},p^{\star})\leq\mathcal{L}(u,p^{\star})\quad\forall\
(u,p)\in\mathbb{R}^{m}\times\mathbb{R}^{n}.$

The saddle point $(u^{*},p^{*})$ satisfies the first order necessary condition
for being the critical point of $\mathcal{L}(u,p)$:

$\displaystyle\nabla f(u^{\star})+B^{\intercal}p^{\star}=0,$ (70a)
$\displaystyle-Bu^{*}+\nabla g(p^{\star})=0.$ (70b)

The first equation (70a) is $\partial_{u}\mathcal{L}(u^{\star},p^{\star})=0$
but the second one (70b) is $-\partial_{p}\mathcal{L}(u^{\star},p^{\star})=0$.
The negative sign is introduced so that (70) is in the form

$\mathcal{A}(x^{\star})=0,$

where $x=(u,p)\in\mathcal{V}\times\mathcal{Q}$, $\nabla
F=\begin{pmatrix}\nabla f&0\\\ 0&\nabla g\end{pmatrix}$, and
$\mathcal{N}=\begin{pmatrix}0&B^{\intercal}\\\ -B&0\end{pmatrix}$. To avoid
confusion, we use the notation $\mathcal{B}^{\rm
sym}=\begin{pmatrix}0&B^{\intercal}\\\ B&0\end{pmatrix}$ in this section. The
splitting (29) becomes

$\mathcal{N}=\begin{pmatrix}0&B^{\intercal}\\\
B&0\end{pmatrix}-\begin{pmatrix}0&0\\\ 2B&0\end{pmatrix},\quad\text{ and
}\quad\mathcal{N}=\begin{pmatrix}I&0\\\ 0&-I\end{pmatrix}\mathcal{B}^{\rm
sym}.$

Therefore $\|\mathcal{N}\|=\|\mathcal{B}^{\rm sym}\|$ for any operator norm.

###  5.2 Strongly-Convex-Strongly-Concave Saddle Point Problem

Given two SPD operators $\mathcal{I}_{\mathcal{V}}$ and
$\mathcal{I}_{\mathcal{Q}}$, we denote by
$\mathcal{I}_{\mu}=\begin{pmatrix}\mu_{f}\mathcal{I}_{\mathcal{V}}&0\\\
0&\mu_{g}\mathcal{I}_{\mathcal{Q}}\end{pmatrix}$ and
$\mathcal{I}_{\mathcal{X}}=\begin{pmatrix}\mathcal{I}_{\mathcal{V}}&0\\\
0&\mathcal{I}_{\mathcal{Q}}\end{pmatrix}$. Then for any
$x=(u,p),y=(v,q)\in\mathcal{V}\times\mathcal{Q},$

$(\mathcal{A}(x)-\mathcal{A}(y),x-y)\geq\|x-y\|_{\mathcal{I}_{\mu}}^{2}\geq\mu\|x-y\|_{\mathcal{I}_{\mathcal{X}}}^{2}$

where $\mu=\min\\{\mu_{f},\mu_{g}\\}$. The accelerated gradient flow and the
discrete schemes follows from discussion in Section 4. The results are
slightly sharp by treating $\mu$ as a block diagonal matrix
$\mathcal{I}_{\mu}$ and including preconditioners
$\mathcal{I}_{\mathcal{V}}^{-1}$ and $\mathcal{I}_{\mathcal{Q}}^{-1}$.

####  5.2.1 The accelerated gradient flow

The component form of the accelerated gradient flow

$\left\\{\begin{aligned} x^{\prime}&=y-x,\\\
y^{\prime}&=x-y-\mathcal{I}_{\mu}^{-1}(\nabla
F(x)+\mathcal{N}y)\end{aligned}\right.$

becomes the preconditioned accelerated gradient flow for the saddle point
system:

$\displaystyle u^{\prime}$ $\displaystyle=v-u,$ (71) $\displaystyle
p^{\prime}$ $\displaystyle=q-p,$ $\displaystyle v^{\prime}$
$\displaystyle=u-v-\frac{1}{\mu_{f}}\mathcal{I}_{\mathcal{V}}^{-1}(\nabla
f(u)+B^{\intercal}q),$ $\displaystyle q^{\prime}$
$\displaystyle=p-q-\frac{1}{\mu_{g}}\mathcal{I}_{\mathcal{Q}}^{-1}(\nabla
g(p)-Bv).$

Consider the Lyapunov function:

$\mathcal{E}(x,y):=D_{F}(x,x^{\star})+\frac{1}{2}\|y-x^{\star}\|_{\mathcal{I}_{\mu}}^{2}=\mathcal{E}^{u}(u,v)+\mathcal{E}^{p}(p,q),$
(72)

with

$\displaystyle\mathcal{E}^{u}(u,v)=D_{f}(u,u^{\star})+\frac{\mu_{f}}{2}\|v-u^{\star}\|_{\mathcal{I}_{\mathcal{V}}}^{2},$
$\displaystyle\mathcal{E}^{p}(p,q)=D_{g}(p,p^{\star})+\frac{\mu_{g}}{2}\|q-p^{\star}\|_{\mathcal{I}_{\mathcal{Q}}}^{2}.$

As $f,g$ are strongly convex, $\mathcal{E}(x,y)\geq 0$ and
$\mathcal{E}(x,y)=0$ iff $x=y=x^{\star}$. Denote the vector field on the right
hand side of (71) by $\mathcal{G}(x,y)$. We have the strong Lyapunov property

$\displaystyle-\nabla\mathcal{E}(x,y)\cdot\mathcal{G}(x,y)\geq\mathcal{E}(x,y)+\frac{1}{2}\|y-x\|^{2}_{\mathcal{I}_{\mu}}.$
(73)

####  5.2.2 Accelerated gradient and skew-symmetric splitting methods

Recall that $x=(u,p),y=(v,q)$. The accelerated gradient and skew-symmetric
splitting method is:

$\displaystyle\frac{\hat{x}_{k+1}-x_{k}}{\alpha}$
$\displaystyle=y_{k}-\hat{x}_{k+1},$ (74a)
$\displaystyle\frac{v_{k+1}-v_{k}}{\alpha}$
$\displaystyle=\hat{u}_{k+1}-v_{k+1}-\frac{1}{\mu_{f}}\mathcal{I}_{\mathcal{V}}^{-1}\left(\nabla
f(\hat{u}_{k+1})+B^{\intercal}q_{k}\right),$ (74b)
$\displaystyle\frac{q_{k+1}-q_{k}}{\alpha}$
$\displaystyle=\hat{p}_{k+1}-q_{k+1}-\frac{1}{\mu_{g}}\mathcal{I}_{\mathcal{Q}}^{-1}\left(\nabla
g(\hat{p}_{k+1})-2Bv_{k+1}+Bv_{k}\right),$ (74c)
$\displaystyle\frac{x_{k+1}-x_{k}}{\alpha}$
$\displaystyle=y_{k+1}-x_{k+1}+\frac{1}{2}(x_{k+1}-\hat{x}_{k+1}).$ (74d)

Each iteration requires $2$ matrix-vector products if we store $Bv_{k}$, $2$
gradient evaluations and the computation of $\mathcal{I}_{\mathcal{V}}^{-1}$
and $\mathcal{I}_{\mathcal{Q}}^{-1}$. Notice that the scheme is explicit as
$v_{k+1}$ can be first updated by (74b) and then used to generate $q_{k+1}$ in
(74c).

Consider the tailored discrete Lyapunov function:

$\displaystyle\mathcal{E}^{\alpha B}(x,y):={}$ $\displaystyle
D_{F}(x,x^{\star})+\frac{1}{2}\|y-x^{\star}\|_{\mathcal{I}_{\mu}-\alpha\mathcal{B}^{\rm
sym}}^{2}.$ (75)

Using the first order necessary conditions (70), the Bregman divergence part
can be related to the duality gap

$\displaystyle D_{F}(x,x^{\star})=$
$\displaystyle~{}D_{f}(u;u^{\star})+D_{g}(p;p^{\star})$ $\displaystyle=$
$\displaystyle~{}f(u)-f(u^{\star})-\langle\nabla
f(u^{\star}),u-u^{\star}\rangle+g(p)-g(p^{\star})-\langle\nabla
g(p^{\star}),p-p^{\star}\rangle$ $\displaystyle=$
$\displaystyle~{}f(u)-f(u^{\star})+\langle
B^{\intercal}p^{\star},u-u^{\star}\rangle+g(p)-g(p^{\star})-\langle
Bu^{\star},p-p^{\star}\rangle$ $\displaystyle=$
$\displaystyle~{}\mathcal{L}(u,p^{\star})-\mathcal{L}(u^{\star},p)\leq\Delta(u,p).$

The convergence analysis will depend on the condition number of the Schur
complement $S=B\mathcal{I}_{\mathcal{V}}^{-1}B^{\intercal}$ in the
$\mathcal{I}_{\mathcal{Q}}$ inner product. We first show $\mathcal{E}^{\alpha
B}$ is non-negative if $\alpha$ is sufficiently small.

###### Lemma 8

Assume $f\in\mathcal{S}_{\mu_{f}}$ and $g\in\mathcal{S}_{\mu_{g}}$. When the
step size

$\alpha\leq\sqrt{\frac{\mu_{f}\mu_{g}}{4L_{S}}}$

where
$L_{S}=\lambda_{\max}(\mathcal{I}_{\mathcal{Q}}^{-1}B\mathcal{I}_{\mathcal{V}}^{-1}B^{\intercal})$,
the matrix

$\mathcal{I}_{\mu}-2\alpha\mathcal{B}^{\rm
sym}=\begin{pmatrix}\mu_{f}\mathcal{I}_{\mathcal{V}}&-2\alpha B^{\intercal}\\\
-2\alpha B&\mu_{g}\mathcal{I}_{\mathcal{Q}}\end{pmatrix}$

is symmetric and positive semidefinite.

###### Proof

We have the block matrix factorization

$\displaystyle\begin{pmatrix}\mu_{f}\mathcal{I}_{\mathcal{V}}&-2\alpha
B^{\intercal}\\\ -2\alpha B&\mu_{g}\mathcal{I}_{\mathcal{Q}}\end{pmatrix}$
(76) $\displaystyle=\begin{pmatrix}I_{m}&0\\\
-\frac{2\alpha}{\mu_{f}}B\mathcal{I}_{\mathcal{V}}^{-1}&I_{n}\end{pmatrix}\begin{pmatrix}\mu_{f}\mathcal{I}_{\mathcal{V}}&0\\\
0&\mu_{g}\mathcal{I}_{\mathcal{Q}}-\frac{4\alpha^{2}}{\mu_{f}}B\mathcal{I}_{\mathcal{V}}^{-1}B^{\intercal}\end{pmatrix}\begin{pmatrix}I_{m}&-\frac{2\alpha}{\mu_{f}}\mathcal{I}_{\mathcal{V}}^{-1}B^{\intercal}\\\
0&I_{n}\end{pmatrix}.$

Then if $\alpha\leq\sqrt{\mu_{f}\mu_{g}/(4L_{S})}$, we have
$\mu_{g}\mathcal{I}_{\mathcal{Q}}-\frac{4\alpha^{2}}{\mu_{f}}B\mathcal{I}_{\mathcal{V}}^{-1}B^{\intercal}\geq
0$ and the results follows.

As a result, $\mathcal{E}^{\alpha B}\geq 0$ if
$\alpha\leq\sqrt{\mu_{f}\mu_{g}/(4L_{S})}$ and $\mathcal{E}^{\alpha B}(x,y)=0$
only if $x=x^{\star}$.

###### Lemma 9

Assume $f\in\mathcal{S}_{\mu_{f}}$ and $g\in\mathcal{S}_{\mu_{g}}$. For
$\mathcal{B}^{\rm sym}=\begin{pmatrix}0&B^{\intercal}\\\ B&0\end{pmatrix}$, we
have

$L_{\mathcal{B}^{\rm sym}}:=\|\mathcal{B}^{\rm
sym}\|_{\mathcal{I}_{\mu}}=\sqrt{\frac{L_{S}}{\mu_{f}\mu_{g}}},\quad\text{and
}\,\kappa_{\mathcal{I}_{\mu}}(\mathcal{N})=\sqrt{\frac{L_{S}}{\mu_{f}\mu_{g}}}.$

###### Proof

Since $\mathcal{B}^{\rm sym}\in\mathbb{R}^{(m+n)^{2}}$ is symmetric,
$\mathcal{B}^{\rm sym}$ has $(m+n)$ real eigenvalues. Let $\lambda$ be the
largest absolute value eigenvalue with respect to $\mathcal{I}_{\mu}$ and
$x=(u,p)$ be the corresponding eigenvector, that is $\|\mathcal{B}^{\rm
sym}\|_{\mathcal{I}_{\mu}}=\lambda$ and

$\mathcal{B}^{\rm sym}x=\lambda\mathcal{I}_{\mu}x.$

The component form follows as

$\displaystyle Bu=\lambda\mu_{g}\mathcal{I}_{\mathcal{Q}}p,\quad
B^{\intercal}p=\lambda\mu_{f}\mathcal{I}_{\mathcal{V}}u.$

Substitute $u$ in the first equation side using the other equation, we get

$\displaystyle
B\mathcal{I}_{\mathcal{V}}^{-1}B^{\intercal}p=\lambda^{2}\mu_{f}\mu_{g}\mathcal{I}_{\mathcal{Q}}p.$

Since $B\mathcal{I}_{\mathcal{V}}^{-1}B^{\intercal}$ is positive definite, we
have $|\lambda|=\sqrt{\frac{L_{S}}{\mu_{f}\mu_{g}}}$ with
$L_{S}=\lambda_{\max}(\mathcal{I}_{\mathcal{Q}}^{-1}B\mathcal{I}_{\mathcal{V}}^{-1}B^{\intercal})$.
Equality on $\kappa_{\mathcal{I}_{\mu}}(\mathcal{N})$ follows since
$\mu=\mu_{f}=\mu_{g}=1$ with respect to $\mathcal{I}_{\mu}$ norm.

Define $\mathcal{E}_{k}^{\alpha B}:=\mathcal{E}^{\alpha B}(x_{k},y_{k})$. With
$\kappa_{\mathcal{I}_{\mu}}(\mathcal{N})$ shown above and $\kappa(F)$ refined
to $\kappa_{\mathcal{I}_{\mathcal{V}}}(f)$ and
$\kappa_{\mathcal{I}_{\mathcal{Q}}}(g)$, we have the following linear
convergence as a direct corollary of Theorem 4.2.

###### Theorem 5.1

Assume $f\in\mathcal{S}_{\mu_{f},L_{f}}$ and $g\in\mathcal{S}_{\mu_{g},L_{g}}$
with $\mu=\min\\{\mu_{f},\mu_{g}\\}>0$. Let $(x_{k},y_{k})$ be the sequence
generated by the accelerated gradient and skew-symmetric splitting method (74)
with arbitrary initial value and step size satisfying

$0\leq\alpha\leq\displaystyle\min\left\\{\sqrt{\frac{\mu_{f}\mu_{g}}{4L_{S}}},\sqrt{\frac{\mu_{f}}{2L_{f}}},\sqrt{\frac{\mu_{g}}{2L_{g}}}\right\\}$

where
$L_{S}=\lambda_{\max}(\mathcal{I}_{\mathcal{Q}}^{-1}B\mathcal{I}_{\mathcal{V}}^{-1}B^{\intercal})$.
Then for the discrete Lyapunov function (75),

$\mathcal{E}_{k+1}^{\alpha B}\leq\frac{1}{1+\alpha/2}\mathcal{E}_{k}^{\alpha
B}.$

In particular, for
$\alpha=1/\max\left\\{2\kappa_{\mathcal{I}_{\mu}}(\mathcal{N}),\sqrt{2\kappa_{\mathcal{I}_{\mathcal{V}}}(f)},\sqrt{2\kappa_{\mathcal{I}_{\mathcal{Q}}}(g)}\right\\}$,
we achieve the accelerated rate

$\displaystyle\|u_{k}-u^{\star}\|_{\mathcal{I}_{\mathcal{V}}}^{2}+\|p_{k}-p^{\star}\|_{\mathcal{I}_{\mathcal{Q}}}^{2}$
$\displaystyle\leq\left(1+1/\max\left\\{4\kappa_{\mathcal{I}_{\mu}}(\mathcal{N}),2\sqrt{2\kappa_{\mathcal{I}_{\mathcal{V}}}(f)},2\sqrt{2\kappa_{\mathcal{I}_{\mathcal{Q}}}(g)}\right\\}\right)^{-k}\frac{2\mathcal{E}_{0}^{\alpha
B}}{\mu}.$

For strongly-convex-strongly-concave saddle point systems, set
$\mathcal{I}_{\mathcal{V}}=I_{m}$ and $\mathcal{I}_{\mathcal{Q}}=I_{n}$,
scheme (74) is an explicit scheme achieving the lower complexity bound:
$\Omega\left(\sqrt{\kappa(f)+\kappa^{2}(\mathcal{N})+\kappa(g)}\cdot|\ln\epsilon_{\rm
out}|\right)$ established in zhang2022lower . Notice that in zhang2022lower ,
only the theoretical lower bound is proved and no algorithms are developed to
match the lower bound. Ours are among the first few explicit schemes achieving
this lower bound.

The condition numbers $\kappa(f),\kappa(g)$, and $\kappa(\mathcal{N})$ are
scaling invariant. They can be improved using appropriate SPD preconditioners
$\mathcal{I}_{\mathcal{V}}$ and $\mathcal{I}_{\mathcal{Q}}$ with the price of
computing $\mathcal{I}_{\mathcal{V}}^{-1}$ and
$\mathcal{I}_{\mathcal{Q}}^{-1}$. The term
$\kappa_{\mathcal{I}_{\mu}}(\mathcal{N})=\sqrt{\frac{L_{S}}{\mu_{f}\mu_{g}}}$
might be the leading term compare with
$\sqrt{\kappa_{\mathcal{I}_{\mathcal{V}}}(f)}$ and
$\sqrt{\kappa_{\mathcal{I}_{\mathcal{Q}}}(g)}$. The observation is that
preconditioner $\mathcal{I}_{\mathcal{Q}}^{-1}$ can be chosen so that $L_{S}$
is small. Namely $\mathcal{I}_{\mathcal{Q}}^{-1}$ is a preconditioner for the
Schur complement $B\mathcal{I}_{\mathcal{V}}^{-1}B^{\intercal}$. For example,
if $\mathcal{I}_{\mathcal{V}}=I_{m}$, then the ideal choice is
$\mathcal{I}_{\mathcal{Q}}=(BB^{\intercal})^{-1}$ so that $L_{S}=1$. Such
choice of $\mathcal{I}_{\mathcal{Q}}$ may increase the condition number
$\kappa_{\mathcal{I}_{\mathcal{Q}}}(g)$ as $(BB^{\intercal})^{-1}$ may not be
a good preconditioner of $g$. Later on we shall show even $\mu_{g}>0$, it is
better to consider the transformed primal dual (TPD) flow proposed in our
recent work ChenWei2022Transformed .

####  5.2.3 Implicit in the skew-symmetric part

The leading term
$\kappa_{\mathcal{I}_{\mu}}(\mathcal{N})=\sqrt{\frac{L_{S}}{\mu_{f}\mu_{g}}}$
is due to the bilinear coupling or equivalently from the explicit treatment
for the skew-symmetric component $\mathcal{N}$. We can treat $\mathcal{N}$
implicitly and obtain the following IMEX scheme:

$\displaystyle\frac{\hat{x}_{k+1}-x_{k}}{\alpha_{k}}$
$\displaystyle=y_{k}-\hat{x}_{k+1},$ (77a)
$\displaystyle\frac{v_{k+1}-v_{k}}{\alpha_{k}}$
$\displaystyle=\hat{u}_{k+1}-v_{k+1}-\frac{1}{\mu_{f}}\mathcal{I}_{\mathcal{V}}^{-1}\left(\nabla
f(\hat{u}_{k+1})+B^{\intercal}q_{k+1}\right),$ (77b)
$\displaystyle\frac{q_{k+1}-q_{k}}{\alpha_{k}}$
$\displaystyle=\hat{p}_{k+1}-q_{k+1}-\frac{1}{\mu_{g}}\mathcal{I}_{\mathcal{Q}}^{-1}\left(\nabla
g(\hat{p}_{k+1})-Bv_{k+1}\right),$ (77c)
$\displaystyle\frac{x_{k+1}-x_{k}}{\alpha_{k}}$
$\displaystyle=y_{k+1}-x_{k+1}.$ (77d)

As the skew-symmetric part is treat implicitly, the restriction
$\alpha\leq\sqrt{\frac{\mu_{f}\mu_{g}}{4L_{S}}}$ can be removed. We state the
convergence theorem directly as the proofs are illustrated with the strong
Lyapunov property in the continuous level and the discussion in Section 4.3
for IMEX schemes.

###### Theorem 5.2

Assume $f\in\mathcal{S}_{\mu_{f},L_{f}}$ and
$g\in\mathcal{S}_{\mu_{g},L_{g}}$. Let $(x_{k},y_{k})$ be the sequence
generated by the preconditioned accelerated gradient method (77) with
arbitrary initial value and step size $\alpha_{k}$ satisfying

$\alpha_{k}^{2}L_{f}\leq(1+\alpha_{k})\mu_{f},\quad\alpha_{k}^{2}L_{g}\leq(1+\alpha_{k})\mu_{g},$

then for the Lyapunov function (72),

$\mathcal{E}_{k+1}\leq\frac{1}{1+\alpha_{k}}\mathcal{E}_{k}.$ (78)

In particualr for
$\alpha_{k}=1/\max\\{\sqrt{\kappa_{\mathcal{I}_{\mathcal{V}}}(f)},\sqrt{\kappa_{\mathcal{I}_{\mathcal{Q}}}(g)}\\}$,
we achieve the accelerated rate

$\mathcal{E}_{k}\leq\left(1+1/\max\\{\sqrt{\kappa_{\mathcal{I}_{\mathcal{V}}}(f)},\sqrt{\kappa_{\mathcal{I}_{\mathcal{Q}}}(g)}\\}\right)^{-k}\mathcal{E}_{0}.$

We discuss the inner solve on (77b)-(77c) which is equivalent to the following
linear algebraic equation

$\begin{pmatrix}(1+\alpha_{k})\mu_{f}\mathcal{I}_{\mathcal{V}}&\alpha_{k}B^{\intercal}\\\
-\alpha_{k}B&(1+\alpha_{k})\mu_{g}\mathcal{I}_{\mathcal{Q}}\end{pmatrix}\begin{pmatrix}v_{k+1}\\\
q_{k+1}\end{pmatrix}=b(\hat{x}_{k+1},y_{k}),$ (79)

with
$b(\hat{x}_{k+1},y_{k})=\mathcal{I}_{\mu}(y_{k}+\alpha_{k}\hat{x}_{k+1})-\alpha_{k}\nabla
F(\hat{x}_{k+1})$. When $\mathcal{I}_{\mathcal{V}},\mathcal{I}_{\mathcal{Q}}$
are identity, solving the equation (79) basically costs the effort of solving

$((1+\alpha_{k})^{2}\mu_{f}\mu_{g}I+\alpha_{k}^{2}BB^{\intercal})x=b,$

which can be solved by conjugate gradient methods with
$\mathcal{O}\left(\frac{\alpha_{k}\|B\|}{(1+\alpha_{k})\sqrt{\mu_{f}\mu_{g}}}\right)$
matrix-vector product. One can also use preconditioned AOR iteration developed
in Section 3.4 for solving (79) but in the $\mathcal{I}_{\mu}$ inner product:

$\displaystyle\frac{v^{\ell+1}-v^{\ell}}{\alpha^{\ell}}$
$\displaystyle=-\frac{1}{\mu_{f}}\mathcal{I}_{\mathcal{V}}^{-1}\left[(1+\alpha_{k})\mu_{f}\mathcal{I}_{\mathcal{V}}v^{\ell+1}+\alpha_{k}B^{\intercal}q^{\ell}-b^{v}(\hat{x}_{k+1},y_{k})\right],$
$\displaystyle\frac{q^{\ell+1}-q^{\ell}}{\alpha^{\ell}}$
$\displaystyle=-\frac{1}{\mu_{g}}\mathcal{I}_{\mathcal{Q}}^{-1}\left[(1+\alpha_{k})\mu_{g}\mathcal{I}_{\mathcal{Q}}q^{\ell+1}+\alpha_{k}B(v^{\ell}-2v^{\ell+1})-b^{q}(\hat{x}_{k+1},y_{k})\right].$

The initial value $(v^{0},q^{0})$ can be set as $(v_{k},q_{k})$. The inner
solver for (79) can be inexact which is a direct application of Section 4.4.
We skip the details here.

Noted that $\hat{x}_{k+1}=(\hat{u}_{k+1},\hat{p}_{k+1}),y_{k}=(v_{k},q_{k})$
are given and not updated in the inner iteration. One AOR iteration is
essentially $2$ matrix-vector products when $\mathcal{I}_{\mathcal{V}}$ and
$\mathcal{I}_{\mathcal{Q}}$ are scaled identities. The inner iteration steps
is again proportional to
$\frac{\alpha_{k}}{1+\alpha_{k}}\sqrt{\frac{L_{S}}{\mu_{f}\mu_{g}}}.$ The
outer iteration complexity is $O(\frac{1}{\alpha_{k}})$. Therefore we conclude
that IMEX scheme (77) requires

  * •

gradient evaluation: $C_{\rm
out}\max\\{\sqrt{\kappa_{\mathcal{I}_{\mathcal{V}}}(f)},\sqrt{\kappa_{\mathcal{I}_{\mathcal{Q}}}(g)}\\}$,

  * •

matrix-vector multiplication: $C_{\rm out}C_{\rm
in}\kappa_{\mathcal{I}_{\mu}}(\mathcal{N})$,

  * •

preconditioners $\mathcal{I}_{\mathcal{V}}^{-1}$ and
$\mathcal{I}_{\mathcal{Q}}^{-1}$: $C_{\rm out}C_{\rm
in}\displaystyle\kappa_{\mathcal{I}_{\mu}}(\mathcal{N})$,

which matches the optimal lower bound in zhang2022lower . The preconditioners
$\mathcal{I}_{\mathcal{V}}^{-1}$ or $\mathcal{I}_{\mathcal{Q}}^{-1}$ are
designed to balance the gradient evaluation cost and matrix-vector
multiplication and preconditioners cost.

####  5.2.4 Implicit in the gradient part

We can use the generalized gradient flow (21) and treat $\nabla F$ implicitly
and $\mathcal{N}$ explicitly with AOR if given the following generalized
proximal operations:

$\displaystyle\operatorname{prox}_{\gamma f}(v)$
$\displaystyle:=\underset{u}{\operatorname{argmin}}~{}f(u)+\frac{1}{2\gamma}\|u-v\|^{2}_{\mathcal{I}_{\mathcal{V}}},$
$\displaystyle\operatorname{prox}_{\sigma g}(q)$
$\displaystyle:=\underset{p}{\operatorname{argmin}}~{}g(p)+\frac{1}{2\sigma}\|p-q\|^{2}_{\mathcal{I}_{\mathcal{Q}}}.$

When $\mathcal{I}_{\mathcal{V}},\mathcal{I}_{\mathcal{Q}}$ are scaled
identity, that matched the typical proximal operators defined in (18). Notice
that the proximal operation works for non-smooth functions. The critical point
of the saddle point system becomes solving for $(u^{\star},p^{\star})$ such
that

$\displaystyle 0$ $\displaystyle\in\partial
f(u^{\star})+B^{\intercal}p^{\star},$ (80) $\displaystyle 0$
$\displaystyle\in-Bu^{*}+\partial g(p^{\star}).$

Here we present the algorithm for non-smooth $f,g$:

$\displaystyle u_{k+1}$
$\displaystyle=\operatorname{prox}_{\alpha/\mu_{f}f}(u_{k}-\mathcal{I}_{\mathcal{V}}^{-1}B^{\intercal}p_{k}),$
(81) $\displaystyle p_{k+1}$
$\displaystyle=\operatorname{prox}_{\alpha/\mu_{g}g}(p_{k}-\mathcal{I}_{\mathcal{Q}}^{-1}B(u_{k}-2u_{k+1})),$

which is equivalent to the iteration

$\displaystyle\frac{u_{k+1}-u_{k}}{\alpha}$
$\displaystyle=-\frac{1}{\mu_{f}}\mathcal{I}_{\mathcal{V}}^{-1}\left(\xi_{k+1}+B^{\intercal}p_{k}\right),\quad\xi_{k+1}\in\partial
f(u_{k+1}),$ (82) $\displaystyle\frac{p_{k+1}-p_{k}}{\alpha}$
$\displaystyle=-\frac{1}{\mu_{g}}\mathcal{I}_{\mathcal{Q}}^{-1}\left(\eta_{k+1}+Bu_{k}-2Bu_{k+1}\right),\quad\eta_{k+1}\in\partial
g(p_{k+1}).$

Consider the Lyapunov function for AOR methods:

$\displaystyle\mathcal{E}^{\alpha B}(u,p)$
$\displaystyle=\frac{\mu_{f}}{2}\|u-u^{\star}\|_{\mathcal{I}_{\mathcal{V}}}^{2}+\frac{\mu_{g}}{2}\|p-p^{\star}\|_{\mathcal{I}_{\mathcal{Q}}}^{2}-2\alpha(B(u-u^{\star}),p-p^{\star})$
(83)
$\displaystyle=\frac{1}{2}\|x-x^{\star}\|^{2}_{\mathcal{I}_{\mu}-2\alpha\mathcal{B}^{\rm
sym}}.$

By Lemma 8, $\mathcal{E}^{\alpha B}(x)\geq 0$ for
$0\leq\alpha<\sqrt{\mu_{f}\mu_{g}/(4L_{S})}$ and $\mathcal{E}^{\alpha B}(x)=0$
if and only if $x=x^{\star}$. As $\nabla F$ is treated implicitly, the
following theorem can be proved similar to Theorem 3.1.

###### Theorem 5.3

Assume $f\in\mathcal{S}_{\mu_{f}}$ and $g\in\mathcal{S}_{\mu_{g}}$. Let
$\\{x_{k}\\}$ be the sequence generated by (81) with arbitrary initial guess
$x_{0}$ and step size
$0\leq\alpha<1/(2\kappa_{\mathcal{I}_{\mu}}(\mathcal{N}))$ with
$\displaystyle\kappa_{\mathcal{I}_{\mu}}(\mathcal{N})=\sqrt{\frac{L_{S}}{\mu_{f}\mu_{g}}}$
and
$L_{S}=\lambda_{\max}(\mathcal{I}_{\mathcal{Q}}^{-1}B\mathcal{I}_{\mathcal{V}}^{-1}B^{\intercal})$.
Then for the Lyapunov function (83),

$\displaystyle\mathcal{E}^{\alpha B}(x_{k+1})\leq$
$\displaystyle~{}\frac{1}{1+\alpha}\mathcal{E}^{\alpha B}(x_{k}).$ (84)

In particular, for $\alpha=\frac{1}{4}\sqrt{\frac{\mu_{f}\mu_{g}}{L_{S}}}$, we
have

$\|x_{k}-x^{\star}\|^{2}\leq\left(1+1/(4\kappa_{\mathcal{I}_{\mu}}(\mathcal{N}))\right)^{-k}3\|x_{0}-x^{\star}\|^{2}.$

This is indeed the optimal bound shown in zhang2022lower . Proximal methods
combining overrelaxation parameters are considered and yielding optimal rates;
see chambolle2011first ; chambolle2016ergodic  for variants of schemes.

###  5.3 Convex-Concave Saddle Point Problems

The acceleration cannot be recovered for the constrained optimization problem
(69) since the corresponding saddle point system (Lagrangian) is not strongly
concave with respect to $p$ (recall that $g(p)=(b,p)$ is linear). However,
(69) admits a unique solution as long as $f$ is strongly convex. In general,
when $f$ is strongly convex, its convex conjugate exists, i.e.,
$f^{*}(\xi)=\max_{u\in\mathbb{R}^{m}}\ \langle\xi,u\rangle-f(u)$ is well
defined and convex. Then (68) is equivalent to the composite optimization
problem without constraints:

$\displaystyle\min_{p\in\mathbb{R}^{n}}f^{*}(-B^{\intercal}p)+g(p),$ (85)

which is a strongly convex optimization problem since $\nabla f$ is Lipschitz
continuous. In this subsection, we adapt the acceleration technique to the
transformed primal-dual flow to develop accelerated gradient methods for
convex-concave saddle point problems.

####  5.3.1 Transformed primal-dual flow

A transformed primal-dual (TPD) flow :

$\displaystyle u^{\prime}$
$\displaystyle=-\mathcal{I}_{\mathcal{V}}^{-1}(\nabla f(u)+B^{\intercal}p),$
(86) $\displaystyle p^{\prime}$
$\displaystyle=-\mathcal{I}_{\mathcal{Q}}^{-1}\left(\nabla
g(p)+B\mathcal{I}_{\mathcal{V}}^{-1}B^{\intercal}p-B(u-\mathcal{I}_{\mathcal{V}}^{-1}\nabla
f(u))\right),$

is proposed in ChenWei2022Transformed , and corresponding numerical
discretization achieve linear convergence for convex-concave saddle point
systems of this kind. The idea is to show with a transformation, it is
equivalent to consider the operator

$\displaystyle\mathcal{A}(x)$ $\displaystyle:=\begin{pmatrix}\nabla
f&B^{\intercal}\\\ -B+B\mathcal{I}_{\mathcal{V}}^{-1}\nabla f&\nabla
g+B\mathcal{I}_{\mathcal{V}}^{-1}B^{\intercal}\end{pmatrix}\begin{pmatrix}u\\\
p\end{pmatrix}$ $\displaystyle=\begin{pmatrix}I&0\\\
B\mathcal{I}_{\mathcal{V}}^{-1}&I\end{pmatrix}\begin{pmatrix}\nabla
f&B^{\intercal}\\\ -B&\nabla g\end{pmatrix}\begin{pmatrix}u\\\ p\end{pmatrix}$

where we extend the matrix-vector product notation to nonlinear version
$(\nabla f,u):=\nabla f(u)$ and $(\nabla g,p):=\nabla g(p)$. Then
$\mathcal{A}$ is a monotone operator under some mild assumption on
$\mathcal{I}_{\mathcal{V}}$. The key is the following estimate on the cross
term.

######  Lemma 10 (Lemma 3.1 in ChenWei2022Transformed )

Suppose $f\in\mathcal{S}_{\mu_{f},L_{f}}$. For any $u_{1},u_{2}\in\mathcal{V}$
and $p_{1},p_{2}\in\mathcal{Q}$, we have

$\displaystyle\langle\nabla f(u_{1})-\nabla
f(u_{2}),\mathcal{I}_{\mathcal{V}}^{-1}B^{\intercal}(p_{1}-p_{2})\rangle$
$\displaystyle\geq{}$
$\displaystyle\frac{\mu_{f}}{2}\|v_{1}-v_{2}\|^{2}_{\mathcal{I}_{\mathcal{V}}}-\frac{L_{f}}{2}\|B^{\intercal}(p_{1}-p_{2})\|^{2}_{\mathcal{I}_{\mathcal{V}}^{-1}}-\frac{1}{2}\langle\nabla
f(u_{1})-\nabla f(u_{2}),u_{1}-u_{2}\rangle,$

where $v=u+\mathcal{I}_{\mathcal{V}}^{-1}B^{\intercal}p$ is a transformed
variable.

For clear illustration, we consider the linear case $\nabla f(u)=H_{f}u$ with
an SPD matrix $H_{f}\geq\mu_{f}\mathcal{I}_{\mathcal{V}}$. Lemma 10 implies
the matrix

$\begin{pmatrix}H_{f}&H_{f}\mathcal{I}_{\mathcal{V}}^{-1}B^{\intercal}\\\
B\mathcal{I}_{\mathcal{V}}^{-1}H_{f}&L_{f}S\end{pmatrix}\geq 0,\quad
S=B\mathcal{I}_{\mathcal{V}}^{-1}B^{\intercal}.$

Then assuming $\nabla g(p)=H_{g}p,H_{g}\geq 0,$ and $L_{f}<2$, we have

$\displaystyle{\rm sym}\mathcal{A}={}$ $\displaystyle{\rm
sym}\begin{pmatrix}H_{f}&B^{\intercal}\\\
-B+B\mathcal{I}_{\mathcal{V}}^{-1}H_{f}&H_{g}+S\end{pmatrix}=\begin{pmatrix}H_{f}&\frac{1}{2}H_{f}\mathcal{I}_{\mathcal{V}}^{-1}B^{\intercal}\\\
\frac{1}{2}B\mathcal{I}_{\mathcal{V}}^{-1}H_{f}&H_{g}+S\end{pmatrix}$
$\displaystyle=\frac{1}{2}\begin{pmatrix}H_{f}&H_{f}\mathcal{I}_{\mathcal{V}}^{-1}B^{\intercal}\\\
B\mathcal{I}_{\mathcal{V}}^{-1}H_{f}&L_{f}S\end{pmatrix}+\frac{1}{2}\begin{pmatrix}H_{f}&0\\\
0&2H_{g}+(2-L_{f})S\end{pmatrix}$
$\displaystyle\geq\frac{\mu}{2}\begin{pmatrix}\mathcal{I}_{\mathcal{V}}&0\\\
0&\mathcal{I}_{\mathcal{Q}}\end{pmatrix}$

with

$\mu=\min\\{\mu_{f},\mu^{+}_{g}\\},$

where the enhanced convexity constant $\mu^{+}_{g}$ is the largest positive
constant s.t.

$2H_{g}+(2-L_{f})S\geq\mu_{g}^{+}\mathcal{I}_{\mathcal{Q}}.$

Switch to the nonlinear case, that is
$2g(\cdot)+(2-L_{f})\frac{1}{2}\|\cdot\|_{S}^{2}\in\mathcal{S}_{\mu_{g}^{+}}$
in the $\mathcal{I}_{\mathcal{Q}}$ inner product. We have the following lower
bounds on $\mu_{g}^{+}$

$\mu_{g}^{+}\geq 2\mu_{g}+(2-L_{f})\mu_{S}\geq(2-L_{f})\mu_{g_{S}},$

where $\mu_{S}=\lambda_{\min}(\mathcal{I}_{\mathcal{Q}}^{-1}S)$ and
$g_{S}(\cdot):=g(\cdot)+\frac{1}{2}\|\cdot\|_{S}^{2}$. Even $\mu_{g}=0$, the
enhanced convexity constant $\mu^{+}_{g}\geq(2-L_{f})\mu_{S}>0$ if given
$L_{f}<2$. The condition $L_{f}<2$ can be always satisfied by rescale $f$ or
$\mathcal{I}_{\mathcal{V}}$.

The TPD flow can be simply written as

$x^{\prime}=-\mathcal{I}_{\mathcal{X}}^{-1}\mathcal{A}(x),$ (87)

where recall that $\mathcal{I}_{\mathcal{X}}={\rm
diag}(\mathcal{I}_{\mathcal{V}},\mathcal{I}_{\mathcal{Q}})$. Consider the
Lyapunov function:

$\mathcal{E}(x)=\mathcal{E}(u,p):=\frac{1}{2}\|u-u^{\star}\|_{\mathcal{I}_{\mathcal{V}}}^{2}+\frac{1}{2}\|p-p^{\star}\|_{\mathcal{I}_{\mathcal{Q}}}^{2}=\frac{1}{2}\|x-x^{\star}\|_{\mathcal{I}_{\mathcal{X}}}^{2}.$
(88)

The strong Lyapunov property has been proved in our recent work
ChenWei2022Transformed . The proof is essentially the way we verify ${\rm
sym}\mathcal{A}\geq\mu\mathcal{I}_{\mathcal{X}}/2$. When adapt to the
nonlinear case, use the convention $(\nabla f,u_{1}-u_{2})=(\nabla
f,u_{1})-(\nabla f,u_{2})=\nabla f(u_{1})-\nabla f(u_{2})$.

######  Lemma 11 (Theorem 3.2 in ChenWei2022Transformed )

Assume function $f\in\mathcal{S}_{\mu_{f},L_{f}}$ with $L_{f}<2$. Then for the
Lyapunov function (88) and the transformed primal-dual gradient flow (87), the
following strong Lyapunov property holds

$\nabla\mathcal{E}(x)\cdot\mathcal{I}_{\mathcal{X}}^{-1}\mathcal{A}(x)\geq\mu\mathcal{E}(x).$
(89)

where $\mu=\min\\{\mu_{f},\mu_{g}^{+}\\}>0$ .

Notice $\mathcal{A}$ cannot be written into the form $\nabla F+\mathcal{N}$
due to the nonlinear transformation term
$B\mathcal{I}_{\mathcal{V}}^{-1}\nabla f$ and thus it is not straightforward
to apply the current framework.

####  5.3.2 Gradient and skew-symmetric splitting methods

We have decomposition $\mathcal{A}=\nabla F+\mathcal{N}+\delta f$ where
$F=f+g_{S}$ and $\delta f=\displaystyle\begin{pmatrix}0&0\\\
B\mathcal{I}_{\mathcal{V}}^{-1}\nabla f&0\end{pmatrix}$. As $\delta f$ is
lower triangular, we can still use the AOR for $\mathcal{N}$ to obtain an
explicit scheme named GSS-TPD:

$\left\\{\begin{aligned}
\frac{u_{k+1}-u_{k}}{\alpha}&=-\mathcal{I}_{\mathcal{V}}^{-1}(\nabla
f(u_{k})+B^{\intercal}p_{k})\\\
\frac{p_{k+1}-p_{k}}{\alpha}&=-\mathcal{I}_{\mathcal{Q}}^{-1}\left[B\mathcal{I}_{\mathcal{V}}^{-1}\nabla
f(u_{k+1})+\nabla g_{S}(p_{k})-B(2u_{k+1}-u_{k})\right].\end{aligned}\right.$
(90)

Consider the Lyapunov function

$\mathcal{E}(x)=\frac{1}{2}\|x-x^{\star}\|^{2}_{\mathcal{I}_{\mathcal{X}}-\alpha\mathcal{B}^{\rm
sym}}-\alpha D_{F}(x^{*},x).$ (91)

######  Theorem 5.4 (Theorem 4.6 in ChenWei2022Transformed )

Suppose $f\in\mathcal{S}_{\mu_{f},L_{f}}$ with $0<\mu_{f}\leq L_{f}<2$ and
$g_{S}(\cdot):=g(\cdot)+\frac{1}{2}\|\cdot\|_{S}^{2}\in\mathcal{S}_{\mu_{g_{S}}}$.
Let $x_{k}=(u_{k},p_{k})$ be generated by AOR iteration (90) with arbitrary
initial value $x_{0}=(u_{0},p_{0})$ and
$\alpha<1/\max\\{2\sqrt{L_{S}},2L_{f},2L_{g_{S}}\\}$. Then for the discrete
Lyapunov function (91), we have

$\displaystyle\mathcal{E}(x_{k+1})\leq$
$\displaystyle~{}\frac{1}{1+\mu\,\alpha/2}\mathcal{E}(x_{k}).$ (92)

where $\mu=\min\\{\mu_{f},\mu_{g}^{+}\\}>0$.

Even for the case $\mu_{g}=0$, we still achieve the linear convergence rate as
$\mu_{g}^{+}>0$. For constrained optimization problem (69), the strong
convexity of $f$ can be further relaxed to the strong convexity of

$f_{\beta}(u)=f(u)+\frac{\beta}{2}\|Bu-b\|^{2}$

by applying TPD to the augmented Lagrangian. Then even $\mu_{f}=0$, if
$\mu_{f_{\beta}}>0$, we can apply GSS-TPD scheme (90) to $f_{\beta}$. However
by adding the augmented term, $\kappa(f)$ is also enlarged to
$\kappa(f_{\beta})$. If we choose relative large $\beta$, $L_{f_{\beta}}$ may
be large and the step size $\alpha$ will be small. See
(ChenWei2022Transformed, , Section 6) for more discussion.

###### Remark 3

We could introduce factor $\mu_{f}^{-1}$ and $(\mu_{g}^{+})^{-1}$ in scheme
(90) and refine the rate to

$\left(1+\min\left\\{\frac{\mu_{f}}{8L_{f}},\frac{\mu_{g}^{+}}{8L_{g_{S}}},\sqrt{\frac{\mu_{f}\mu_{g}^{+}}{16L_{S}}}\right\\}\right)^{-1}$
(93)

for the step size
$\alpha=\min\left\\{\frac{\mu_{f}}{4L_{f}},\frac{\mu_{g}^{+}}{4L_{g_{S}}},\sqrt{\frac{\mu_{f}\mu_{g}^{+}}{4L_{S}}}\right\\}$.
The advantage of using (90) is that only Lipschitz constants are needed to be
estimated while (74) requires estimate on $\mu_{f}$ and $\mu_{g_{S}}$, which
is usually harder to obtain.

####  5.3.3 Accelerated transformed primal-dual gradient flow

We propose the accelerated transformed primal-dual gradient flow:

$\displaystyle u^{\prime}$ $\displaystyle=v-u,$ (94) $\displaystyle
v^{\prime}$
$\displaystyle=\frac{1}{2}(u-v)-\frac{1}{\mu_{f}}\mathcal{I}_{\mathcal{V}}^{-1}(\nabla
f(u)+B^{\intercal}q),$ $\displaystyle p^{\prime}$ $\displaystyle=q-p,$
$\displaystyle q^{\prime}$
$\displaystyle=p-q-\mathcal{I}_{\mathcal{Q}}^{-1}(\nabla
g_{S}(p)-Bv+B\mathcal{I}_{\mathcal{V}}^{-1}\nabla f(u)),$

where recall that $g_{S}(p)=g(p)+\frac{1}{2}\|p\|_{S}^{2}$. Let $x=(u,p)$ and
$y=(v,q)$. Consider the Lyapunov function:

$\mathcal{E}(x,y):=D_{f}(u,u^{\star})+\frac{\mu_{f}}{2}\|v-u^{\star}\|_{\mathcal{I}_{\mathcal{V}}}^{2}+D_{g_{S}}(p,p^{\star})+\frac{1}{2}\|q-p^{\star}\|_{\mathcal{I}_{\mathcal{Q}}}^{2}.$
(95)

As $f,g_{S}$ are strongly convex, $\mathcal{E}(x,y)\geq 0$ and
$\mathcal{E}(x,y)=0$ iff $x=y=x^{\star}$.

Denote the vector field on the right hand side of (94) by $\mathcal{G}(x,y)$.
In order to verify the strong Lyapunov property, we first show the matrix
calculation when $\nabla f(u)=H_{f}u,\nabla g(p)=H_{g}p$ are linear with
$H_{f}=\nabla^{2}f,H_{g}=\nabla^{2}g$ being constant SPD matrices. Then
$H_{g_{S}}=H_{g}+B\mathcal{I}_{\mathcal{V}}^{-1}B^{\intercal}$. As
$-\nabla\mathcal{E}(x,y)\cdot\mathcal{G}(x,y)$ is a quadratic form of
$(x-x^{\star},y-x^{\star})^{\intercal}$, we calculate the corresponding matrix
in the order of $(u,v,p,q)$ as

$\displaystyle\begin{pmatrix}H_{f}&0&0&0\\\
0&\mu_{f}\mathcal{I}_{\mathcal{V}}&0&0\\\ 0&0&H_{g_{S}}&0\\\
0&0&0&\mathcal{I}_{\mathcal{Q}}\end{pmatrix}\begin{pmatrix}I&-I&0&0\\\
-I/2+\mathcal{I}_{\mathcal{V}}^{-1}H_{f}/\mu_{f}&I/2&0&\mathcal{I}_{\mathcal{V}}^{-1}B^{\intercal}/\mu_{f}\\\
0&0&I&-I\\\
\mathcal{I}_{\mathcal{Q}}^{-1}B\mathcal{I}_{\mathcal{V}}^{-1}H_{f}&-\mathcal{I}_{\mathcal{Q}}^{-1}B&-I+\mathcal{I}_{\mathcal{Q}}^{-1}H_{g_{S}}&I\end{pmatrix}$
$\displaystyle=$ $\displaystyle~{}\begin{pmatrix}H_{f}&-H_{f}&0&0\\\
-\mu_{f}\mathcal{I}_{\mathcal{V}}/2+H_{f}&\mu_{f}\mathcal{I}_{\mathcal{V}}/2&0&B^{\intercal}\\\
0&0&H_{g_{S}}&-H_{g_{S}}\\\
B\mathcal{I}_{\mathcal{V}}^{-1}H_{f}&-B&-\mathcal{I}_{\mathcal{Q}}+H_{g_{S}}&\mathcal{I}_{\mathcal{Q}}\end{pmatrix}.$

For a quadratic form, $x^{\intercal}Mx=x^{\intercal}\operatorname{sym}(M)x$.
So we calculate its symmetric part

$\displaystyle\operatorname{sym}\begin{pmatrix}H_{f}&-H_{f}&0&0\\\
-\mu_{f}\mathcal{I}_{\mathcal{V}}/2+H_{f}&\mu_{f}\mathcal{I}_{\mathcal{V}}/2&0&B^{\intercal}\\\
0&0&H_{g_{S}}&-H_{g_{S}}\\\
B\mathcal{I}_{\mathcal{V}}^{-1}H_{f}&-B&-\mathcal{I}_{\mathcal{Q}}+H_{g_{S}}&\mathcal{I}_{\mathcal{Q}}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}H_{f}&-\mu_{f}\mathcal{I}_{\mathcal{V}}/4&0&B\mathcal{I}_{\mathcal{V}}^{-1}H_{f}/2\\\
-\mu_{f}\mathcal{I}_{\mathcal{V}}/4&\mu_{f}\mathcal{I}_{\mathcal{V}}/2&0&0\\\
0&0&H_{g_{S}}&-\mathcal{I}_{\mathcal{Q}}/2\\\
B\mathcal{I}_{\mathcal{V}}^{-1}H_{f}/2&0&-\mathcal{I}_{\mathcal{Q}}/2&\mathcal{I}_{\mathcal{Q}}\end{pmatrix}$
$\displaystyle=$
$\displaystyle~{}\begin{pmatrix}H_{f}/2&-\mu_{f}\mathcal{I}_{\mathcal{V}}/4&0&0\\\
-\mu_{f}\mathcal{I}_{\mathcal{V}}/4&\mu_{f}\mathcal{I}_{\mathcal{V}}/2&0&0\\\
0&0&H_{g_{S}}&-\mathcal{I}_{\mathcal{Q}}/2\\\
0&0&-\mathcal{I}_{\mathcal{Q}}/2&3\mathcal{I}_{\mathcal{Q}}/4\end{pmatrix}+\begin{pmatrix}H_{f}/2&0&0&B\mathcal{I}_{\mathcal{V}}^{-1}H_{f}/2\\\
0&0&0&0\\\ 0&0&0&0\\\
B\mathcal{I}_{\mathcal{V}}^{-1}H_{f}/2&0&0&\mathcal{I}_{\mathcal{Q}}/4\end{pmatrix}$
$\displaystyle\geq$
$\displaystyle~{}\frac{1}{4}\operatorname{diag}\\{H_{f},\mu_{f}\mathcal{I}_{\mathcal{V}},H_{g_{S}},\mathcal{I}_{\mathcal{Q}}\\}.$

For the last inequality, for $(u,v)$-block, as
$H_{f}\geq\mu_{f}\mathcal{I}_{\mathcal{V}}$:

$\displaystyle\begin{pmatrix}H_{f}/2&-\mu_{f}\mathcal{I}_{\mathcal{V}}/4\\\
-\mu_{f}\mathcal{I}_{\mathcal{V}}/4&\mu_{f}\mathcal{I}_{\mathcal{V}}/2\end{pmatrix}$
$\displaystyle=\frac{1}{4}\begin{pmatrix}H_{f}&-\mu_{f}\mathcal{I}_{\mathcal{V}}\\\
-\mu_{f}\mathcal{I}_{\mathcal{V}}&\mu_{f}\mathcal{I}_{\mathcal{V}}\end{pmatrix}+\frac{1}{4}\begin{pmatrix}H_{f}&0\\\
0&\mu_{f}\mathcal{I}_{\mathcal{V}}\end{pmatrix}$ (96)
$\displaystyle\geq\frac{1}{4}\begin{pmatrix}H_{f}&0\\\
0&\mu_{f}\mathcal{I}_{\mathcal{V}}\end{pmatrix}.$

For $(p,q)$-block, if $3H_{g_{S}}/2\geq\mathcal{I}_{\mathcal{Q}}$:

$\displaystyle\begin{pmatrix}H_{g_{S}}&-\mathcal{I}_{\mathcal{Q}}/2\\\
-\mathcal{I}_{\mathcal{Q}}/2&3\mathcal{I}_{\mathcal{Q}}/4\end{pmatrix}$
$\displaystyle=\frac{1}{2}\begin{pmatrix}3H_{g_{S}}/2&-\mathcal{I}_{\mathcal{Q}}\\\
-\mathcal{I}_{\mathcal{Q}}&\mathcal{I}_{\mathcal{Q}}\end{pmatrix}+\frac{1}{4}\begin{pmatrix}H_{g_{S}}&0\\\
0&\mathcal{I}_{\mathcal{Q}}\end{pmatrix}$ (97)
$\displaystyle\geq\frac{1}{4}\begin{pmatrix}H_{g_{S}}&0\\\
0&\mathcal{I}_{\mathcal{Q}}\end{pmatrix}.$

And the rest term for $(u,q)$-block, if $\mathcal{I}_{\mathcal{Q}}\geq
2L_{f}S$, we apply Lemma 10:

$\begin{pmatrix}H_{f}/2&B\mathcal{I}_{\mathcal{V}}^{-1}H_{f}/2\\\
B\mathcal{I}_{\mathcal{V}}^{-1}H_{f}/2&\mathcal{I}_{\mathcal{Q}}/4\end{pmatrix}=\frac{1}{2}\begin{pmatrix}H_{f}&B\mathcal{I}_{\mathcal{V}}^{-1}H_{f}\\\
B\mathcal{I}_{\mathcal{V}}^{-1}H_{f}&L_{f}S\end{pmatrix}+\frac{1}{2}\begin{pmatrix}0&0\\\
0&\mathcal{I}_{\mathcal{Q}}/2-L_{f}S\end{pmatrix}\geq 0.$

We give parameter choice that satisfying the required conditions and verify
the strong Lyapunov property in the following lemma.

###### Lemma 12

Assume $f\in\mathcal{S}_{\mu_{f},L_{f}}$ with $L_{f}\leq 3/4$ and
$g_{S}\in\mathcal{S}_{\mu_{g_{S}}}$. Suppose we choose
$\mathcal{I}_{\mathcal{Q}}$ such that

$2/3-\mu_{g}\leq\lambda_{\min}(\mathcal{I}_{\mathcal{Q}}^{-1}S)\leq\lambda_{\max}(\mathcal{I}_{\mathcal{Q}}^{-1}S)\leq
1/(2L_{f}).$

Then for the Lyapunov function (95) and the accelerated transformed primal-
dual flow vector field $\mathcal{G}$ defined as the right side of (94), the
following strong Lyapunov property holds

$\displaystyle-\nabla\mathcal{E}(x,y)\cdot\mathcal{G}(x,y)$
$\displaystyle\geq\frac{1}{2}\mathcal{E}(x,y).$ (98)

###### Proof

We expand the product in four terms using
$\mathcal{G}=(\mathcal{G}^{u},\mathcal{G}^{v},\mathcal{G}^{p},\mathcal{G}^{q})$:

$\displaystyle-\nabla\mathcal{E}\cdot\mathcal{G}$
$\displaystyle=-\partial_{u}\mathcal{E}\cdot\mathcal{G}^{u}-\partial_{v}\mathcal{E}\cdot\mathcal{G}^{v}-\partial_{p}\mathcal{E}\cdot\mathcal{G}^{p}-\partial_{q}\mathcal{E}\cdot\mathcal{G}^{q}.$

As $\mathcal{G}(x^{\star},x^{\star})=0$, we can insert
$\mathcal{G}(x^{\star},x^{\star})$ into these terms. For the $(u,v)$ terms,
direct computation gives

$\displaystyle-\partial_{u}\mathcal{E}\cdot\mathcal{G}^{u}-\partial_{v}\mathcal{E}\cdot\mathcal{G}^{v}={}$
$\displaystyle\langle\nabla f(u)-\nabla
f(u^{\star}),u-u^{\star}\rangle+\mu_{f}/2\langle
v-u^{\star},u-v\rangle_{\mathcal{I}_{\mathcal{V}}}$ $\displaystyle+\langle
v-u^{\star},B^{\intercal}(q-p^{\star})\rangle$

and

$\displaystyle-\partial_{p}\mathcal{E}\cdot\mathcal{G}^{p}-\partial_{q}\mathcal{E}\cdot\mathcal{G}^{q}={}$
$\displaystyle\langle\nabla g_{S}(p)-\nabla
g_{S}(p^{\star}),p-p^{\star}\rangle+\langle
q-p^{\star},p-q\rangle_{\mathcal{I}_{\mathcal{Q}}}$ $\displaystyle-\langle
q-p^{\star},B(v-u^{\star})\rangle+\langle
q-p^{\star},B\mathcal{I}_{\mathcal{V}}^{-1}(\nabla f(u)-\nabla
f(u^{\star}))\rangle.$

By Lemma 10,

$\displaystyle\langle q-p^{\star},B\mathcal{I}_{\mathcal{V}}^{-1}(\nabla
f(u)-\nabla f(u^{\star}))\rangle$ $\displaystyle\geq$
$\displaystyle~{}\frac{\mu_{f}}{2}\|v-v^{\star}\|_{\mathcal{I}_{\mathcal{V}}}^{2}-\frac{L_{f}}{2}\|q-p^{\star}\|^{2}_{S}-\frac{1}{2}\langle\nabla
f(u)-\nabla f(u^{\star}),u-u^{\star}\rangle.$

Adding them together, we get

$\displaystyle-\nabla\mathcal{E}\cdot\mathcal{G}\geq{}$
$\displaystyle\frac{1}{2}\langle\nabla f(u)-\nabla
f(u^{\star}),u-u^{\star}\rangle-\frac{\mu_{f}}{2}\langle
v-u^{\star},u-v\rangle_{\mathcal{I}_{\mathcal{V}}}$
$\displaystyle+\langle\nabla g_{S}(p)-\nabla
g_{S}(p^{\star}),p-p^{\star}\rangle-\frac{L_{f}}{2}\|q-q^{\star}\|^{2}_{S}$
$\displaystyle-\langle q-p^{\star},p-q\rangle_{\mathcal{I}_{\mathcal{Q}}}$

We use the identity for squares (34) to expand

$(q-p^{\star},p-q)_{\mathcal{I}_{\mathcal{Q}}}=\frac{1}{2}(\|p-p^{\star}\|_{\mathcal{I}_{\mathcal{Q}}}^{2}-\|q-p^{\star}\|_{\mathcal{I}_{\mathcal{Q}}}^{2}-\|q-p\|_{\mathcal{I}_{\mathcal{Q}}}^{2}),$

and split

$\langle\nabla g_{S}(p)-\nabla
g_{S}(p^{\star}),p-p^{\star}\rangle=D_{g_{S}}(p;p^{\star})+D_{g_{S}}(p^{\star};p).$

Observe that $\lambda_{\min}(\mathcal{I}_{\mathcal{Q}}^{-1}S)\geq 2/3-\mu_{g}$
implies $\mu_{g}\mathcal{I}_{\mathcal{Q}}+S\geq 2\mathcal{I}_{\mathcal{Q}}/3,$
and consequently

$D_{g_{S}}(p;p^{\star})+\frac{1}{2}D_{g_{S}}(p^{\star};p)\geq\frac{3}{4}\|p-p^{\star}\|_{\mu_{g}\mathcal{I}_{\mathcal{Q}}+S}^{2}\geq\frac{1}{2}\|p-p^{\star}\|_{\mathcal{I}_{\mathcal{Q}}}^{2}.$

And $\lambda_{\max}(\mathcal{I}_{\mathcal{Q}}^{-1}S)\leq 1/(2L_{f})$ implies
$2L_{f}S\leq\mathcal{I}_{\mathcal{Q}}$. We have

$\displaystyle\frac{1}{4}\|q-p^{\star}\|_{\mathcal{I}_{\mathcal{Q}}}^{2}$
$\displaystyle\geq\frac{L_{f}}{2}\|q-p^{\star}\|_{S}^{2}.$

The terms involving $(u,v)$ can be bounded by (47) and the identity for
squares (34). Therefore,

$\displaystyle-\nabla\mathcal{E}\cdot\mathcal{G}\geq{}$
$\displaystyle\frac{1}{2}D_{f}(u;u^{\star})+\frac{\mu_{f}}{4}\|v-u^{\star}\|_{\mathcal{I}_{\mathcal{V}}}^{2}+\frac{1}{2}D_{g_{S}}(p;p^{\star})+\frac{1}{4}\|q-p^{\star}\|_{\mathcal{I}_{\mathcal{Q}}}^{2}$
$\displaystyle+\frac{\mu_{f}}{4}\|v-u\|_{\mathcal{I}_{\mathcal{V}}}^{2}+\frac{1}{2}\|p-q\|_{\mathcal{I}_{\mathcal{Q}}}^{2}.$

Dropping the last two (nonnegative) quadratic terms we get (98).

Again the condition on $L_{f}$ can be easily fulfilled by rescaling $f$ or
$\mathcal{I}_{\mathcal{V}}$. In the next subsection, we will see that the
condition number of $\mathcal{I}_{\mathcal{Q}}^{-1}S$ will enter the
convergence rates of numerical schemes.

####  5.3.4 Accelerated transformed primal-dual method

We propose an accelerated transformed primal-dual (ATPD) method:

$\displaystyle\frac{\hat{x}_{k+1}-x_{k}}{\alpha}$
$\displaystyle=y_{k}-\hat{x}_{k+1},$ (99a)
$\displaystyle\frac{v_{k+1}-v_{k}}{\alpha}$
$\displaystyle=\frac{1}{2}(\hat{u}_{k+1}-v_{k+1})-\frac{1}{\mu_{f}}\mathcal{I}_{\mathcal{V}}^{-1}(\nabla
f(\hat{u}_{k+1})+B^{\intercal}q_{k}),$ (99b)
$\displaystyle\frac{q_{k+1}-q_{k}}{\alpha}$
$\displaystyle=\hat{p}_{k+1}-q_{k+1}-\mathcal{I}_{\mathcal{Q}}^{-1}\left[\nabla
g_{S}(p_{k+1})+Bv_{k}-2Bv_{k+1}+B\mathcal{I}_{\mathcal{V}}^{-1}f(u_{k+1})\right],$
(99c) $\displaystyle\frac{x_{k+1}-\hat{x}_{k+1}}{\alpha}$
$\displaystyle=(y_{k+1}-y_{k})-\frac{1}{4}(x_{k+1}-\hat{x}_{k+1}).$ (99d)

Recall that $x=(u,p),y=(v,q)$. Approximation $(\hat{u}_{k+1},\hat{p}_{k+1})$
is first updated by (99a) and then used to update $(v_{k+1},q_{k+1})$ by a TPD
iteration. The last step is an extrapolation to produce $(u_{k+1},p_{k+1})$.

Denote $\mathcal{I}_{\mu}=\begin{pmatrix}\mu_{f}\mathcal{I}_{\mathcal{V}}&0\\\
0&\mathcal{I}_{\mathcal{Q}}\end{pmatrix}$. Consider the tailored discrete
Lyapunov function:

$\displaystyle\mathcal{E}^{\alpha B}_{k}=\mathcal{E}^{\alpha
B}(x_{k},y_{k}):={}$ $\displaystyle
D_{f}(u_{k},u^{\star})+D_{g_{S}}(p_{k},p^{\star})+\frac{1}{2}\|y_{k}-x^{\star}\|_{\mathcal{I}_{\mu}-\alpha\mathcal{B}^{\operatorname{sym}}}^{2}.$
(100)

According to Lemma 8 with $\mu_{g}=1$, $\mathcal{E}^{\alpha B}\geq 0$ for
$0\leq\alpha\leq\sqrt{\mu_{f}/(4L_{S})}$ and $\mathcal{E}^{\alpha B}(x)=0$
only if $x=x^{\star}$.

###### Theorem 5.5

Assume $f\in\mathcal{S}_{\mu_{f},L_{f}}$ with $L_{f}\leq 3/4$ and
$g_{S}\in\mathcal{S}_{\mu_{g_{S}},L_{g_{S}}}$. Suppose we can choose
$\mathcal{I}_{\mathcal{Q}}$ such that

$2/3-\mu_{g}\leq\lambda_{\min}(\mathcal{I}_{\mathcal{Q}}^{-1}S)\leq\lambda_{\max}(\mathcal{I}_{\mathcal{Q}}^{-1}S)\leq
1/(2L_{f}).$ (101)

Let $(x_{k},y_{k})$ be the sequence generated by the accelerated transformed
primal-dual gradient method (99) with arbitrary initial guess and

$0<\alpha\leq\displaystyle\min\left\\{\sqrt{\frac{\mu_{f}}{4L_{S}}},\sqrt{\frac{\mu_{f}}{2L_{f}}},\sqrt{\frac{1}{2L_{g_{S}}}}\right\\}.$

Then for the modified Lyapunov function (100) and $k\geq 0$,

$\mathcal{E}^{\alpha B}_{k+1}\leq\frac{1}{1+\alpha/4}\mathcal{E}^{\alpha
B}_{k}.$

Let us discuss the condition (101). Given two SPD operators
$\mathcal{I}_{\mathcal{V}}$ and $\mathcal{I}_{\mathcal{Q}}$, we shall show
that we can always rescale to $c_{v}\mathcal{I}_{\mathcal{V}}$ and
$c_{q}\mathcal{I}_{\mathcal{Q}}$ such that condition (101) holds. For clarity,
we denote the Lipschitz constant of $\nabla f$ w.r.t.
$c_{v}\mathcal{I}_{\mathcal{V}}$ by $L_{f,c_{v}}$ and $L_{f}=L_{f,1}$ for the
non-scaled one. Similar notation apply to $L_{g,c_{q}}$. Then we have the
relation $L_{f,c_{v}}=c_{v}^{-1}L_{f}$ and $L_{g,c_{q}}=c_{q}^{-1}L_{g}$.
Denoted by
$\lambda_{\max}(c_{v},c_{q})=\lambda_{\max}((c_{q}\mathcal{I}_{\mathcal{Q}})^{-1}B(c_{v}\mathcal{I}_{\mathcal{V}})^{-1}B^{\intercal})$
and abbreviate $\lambda_{\max}(1,1)=\lambda_{\max}$. Then we have the scaling
relation $\lambda_{\max}(c_{v},c_{q})=(c_{v}\,c_{q})^{-1}\lambda_{\max}$.
Similar notation is applied to $\lambda_{\min}$. In the sequel
$\kappa_{\mathcal{I}_{\mathcal{Q}}}(S)=\kappa(\mathcal{I}_{\mathcal{Q}}^{-1}B\mathcal{I}_{\mathcal{V}}^{-1}B^{\intercal})$.
Those non-scaled quantities are considered as known and fixed. We investigate
the scaling effect.

First of all, as $L_{f,c_{v}}$ is inversely proportional to $c_{v}$, we can
choose

$c_{v}=\frac{4}{3}L_{f}\kappa_{\mathcal{I}_{\mathcal{Q}}}(S),\text{ s. t.
}\quad
L_{f,c_{v}}=\frac{3}{4}\frac{1}{\kappa_{\mathcal{I}_{\mathcal{Q}}}(S)}\leq\frac{3}{4}.$

With $c_{v}$ determined, we choose $c_{q}$ s.t.

$(c_{v}\,c_{q})^{-1}\lambda_{\min}=\frac{2}{3}.$

By construction, we have the desired bound in the
$c_{v}\mathcal{I}_{\mathcal{V}}$ and $c_{q}\mathcal{I}_{\mathcal{Q}}$ inner
products

$2/3-\mu_{g,c_{q}}\leq\lambda_{\min}(c_{v},c_{q})\leq\lambda_{\max}(c_{v},c_{q})\leq
1/(2L_{f,c_{v}}).$

Then we estimate the bound of iteration steps which is proportional to
$1/\alpha$

$|\ln\epsilon_{\rm
out}|\max\left\\{\sqrt{\frac{4L_{S}(c_{v},c_{q})}{\mu_{f,c_{v}}}},\sqrt{\frac{2L_{f,c_{v}}}{\mu_{f,c_{v}}}},\sqrt{2L_{g_{S},c_{q}}}\right\\},$

where $L_{S}(c_{v},c_{q})=\lambda_{\max}(c_{v},c_{q})$. Write
$L_{S}(c_{v},c_{q})/\mu_{f,c_{v}}=L_{S}(c_{v},c_{q})/L_{f,c_{v}}\kappa_{\mathcal{I}_{\mathcal{V}}}(f)$
and use
$L_{S}(c_{v},c_{q})=\lambda_{\max}(c_{v},c_{q})=(c_{v}c_{q})^{-1}\lambda_{\max}=\frac{2}{3}\kappa_{\mathcal{I}_{\mathcal{Q}}}(S)$
and $1/L_{f,c_{v}}=\frac{4}{3}\kappa_{\mathcal{I}_{\mathcal{Q}}}(S)$ to bound
the term

$\sqrt{\frac{4L_{S}(c_{v},c_{q})}{\mu_{f,c_{v}}}}=\frac{4\sqrt{2}}{3}\sqrt{\kappa_{\mathcal{I}_{\mathcal{V}}}(f)}\,\kappa_{\mathcal{I}_{\mathcal{Q}}}(S).$

Condition (101) implies $\mu_{g_{S},c_{q}}\geq 2/3$ and consequently
$\sqrt{L_{g_{S},c_{q}}}=\frac{2}{\sqrt{3}}\sqrt{\kappa_{\mathcal{I}_{\mathcal{Q}}}(g_{S})}$.
So we get the scaling invariant upper bound of iteration complexity

$|\ln\epsilon_{\rm
out}|\max\left\\{\frac{4\sqrt{2}}{3}\sqrt{\kappa_{\mathcal{I}_{\mathcal{V}}}(f)}\,\kappa_{\mathcal{I}_{\mathcal{Q}}}(S),\frac{2\sqrt{3}}{3}\sqrt{\kappa_{\mathcal{I}_{\mathcal{Q}}}(g_{S})}\right\\}.$
(102)

In particular, for affinely constrained optimization problems, $g_{S}=S$, we
have an explicit scheme with complexity

$\mathcal{O}\left(|\ln\epsilon_{\rm
out}|\sqrt{\kappa_{\mathcal{I}_{\mathcal{V}}}(f)}\,\kappa_{\mathcal{I}_{\mathcal{Q}}}(S)\right).$
(103)

The preconditioner $\mathcal{I}_{\mathcal{V}}$ is designed s.t.
$\kappa_{\mathcal{I}_{\mathcal{V}}}(f)$ is small. The preconditioner
$\mathcal{I}_{\mathcal{Q}}$ is chosen s.t.
$\kappa_{\mathcal{I}_{\mathcal{Q}}}(S)$ is small. Again the ideal choice is
$\mathcal{I}_{\mathcal{Q}}=S$ and $\kappa_{\mathcal{I}_{\mathcal{Q}}}(S)=1$.
For that purpose, sometimes we may choose simple $\mathcal{I}_{\mathcal{V}}$,
so that $S^{-1}=(B\mathcal{I}_{\mathcal{V}}^{-1}B^{\intercal})^{-1}$ is easy
to compute or approximate. For example, $\mathcal{I}_{\mathcal{V}}=I$ and
$\mathcal{I}_{\mathcal{Q}}$ can be obtained by few multigrid cycles or
incomplete LU factorization of the SPD matrix $BB^{\intercal}$.

Iterative methods for computing $S^{-1}r$ can be thought of as an inner
iteration which can be a linear inner solver or a non-linear one. For linear
solvers, there is no need to computing $S^{-1}r$ accurately. For any
convergent method, say $\|I-\mathcal{I}_{\mathcal{Q}}^{-1}S\|\leq 1-\delta<1$,
we can bound the condition number $\kappa_{\mathcal{I}_{\mathcal{Q}}}(S)\leq
2/\delta-1.$ For example, if $\delta=2/3$, then
$\kappa_{\mathcal{I}_{\mathcal{Q}}}(S)\leq 2$. If $\delta$ is small, we may
use acceleration for the inner linear iterative methods.

The inner solver for $S^{-1}$ can be also a nonlinear one, e.g., the conjugate
gradient (CG) method. The bound (103), however, cannot be applied as now the
operator associated to CG is non-linear. Instead, the perturbation argument on
controlling the residual $\varepsilon_{\rm in}$ developed in Section 4.4 can
be applied. The inner iteration complexity will be
$\mathcal{O}(|\ln\varepsilon_{\rm
in}|\sqrt{\kappa_{\mathcal{I}_{\mathcal{Q}}}(S)})$ where a linear operator
$\mathcal{I}_{\mathcal{Q}}$ is used in the inner iteration as a preconditioner
of $S$. Therefore the total iteration complexity for the accelerated
transformed primal-dual method is

$\mathcal{O}\left(|\ln\epsilon_{\rm out}||\ln\varepsilon_{\rm
in}|\sqrt{\kappa_{\mathcal{I}_{\mathcal{V}}}(f)\,\kappa_{\mathcal{I}_{\mathcal{Q}}}(S)}\right),$

which achieves the optimal complexity bound for affinely constrained problems
salim2022optimal  provided $|\ln\varepsilon_{\rm in}|$ is not too large. Again
the convexity of $f$ can be relaxed to $f_{\beta}$ if we consider the
augmented Lagrangian.

For general strongly-convex-concave saddle point problems, as
$\lambda_{\min}=2/3$ and $\lambda_{\max}\leq L_{g_{S}}$, the bound (102) can
be bounded by

$\mathcal{O}\left(|\ln\epsilon_{\rm
out}|\sqrt{\kappa_{\mathcal{I}_{\mathcal{V}}}(f)}\,\kappa_{\mathcal{I}_{\mathcal{Q}}}(g_{S})\right).$
(104)

which relaxes the leading term $\kappa_{\mathcal{I}_{\mu}}(\mathcal{N})$ in
Theorem 5.1. The preconditioner $\mathcal{I}_{\mathcal{V}}$ is designed s.t.
$\kappa_{\mathcal{I}_{\mathcal{V}}}(f)$ is small and
$\mathcal{I}_{\mathcal{Q}}$ is for
$\kappa_{\mathcal{I}_{\mathcal{Q}}}(g_{S})$. Similar discussion on linear and
non-linear solvers still holds.

##  6 Conclusion

In this paper, we propose GSS methods and AGSS methods for solving a class of
strongly monotone operator equation achieving accelerated linear convergence
rates. The proof of strong Lyapunov property bridges the design of dynamical
system on the continuous level, the choice of Lyapunov function and the linear
convergence of the numerical schemes (as discretization of the continuous
flow). As direct applications, we derive optimal algorithms for strongly-
convex-strongly-concave saddle point systems with bilinear coupling. Combining
the transformed primal-dual methods in our recent work ChenWei2022Transformed
and augmented Lagrangian, accelerated linear convergence rates can be retained
for general convex-concave saddle point problems.

As it is well-known in convex minimization, the optimal mixed type convergence
rates is $\mathcal{O}(\min\\{1/k^{2},(1-1/\sqrt{\kappa(f)})^{-k})\\}$
chen2021unified , which enjoys a sublinear rate as long as $f$ is convex. With
more careful design of parameters, our framework may achieve sublinear
convergence rates for monotone operator equations. Then the iterative methods
with mixed type convergence rate can be smoother for deriving nonlinear
multigrid methods for solving saddle point problems, yielding iteration
complexity free of problem size and condition number.

###### Acknowledgements.

The authors would like to thank Dr. Hao Luo for fruitful discussion.

## Funding

L. Chen and J. Wei are supported by National Science Fundation DMS-2012465.

## Conflict of Interest

The authors have no conflicts of interest to declare that are relevant to the
content of this article.

## References

  * (1) Attouch, H., Peypouquet, J.: Convergence of inertial dynamics and proximal algorithms governed by maximally monotone operators.  Mathematical Programming 174(1), 391–432 (2019)
  * (2) Bai, Z.Z., Golub, G.H.: Accelerated Hermitian and skew-Hermitian splitting iteration methods for saddle-point problems.  IMA Journal of Numerical Analysis 27(1), 1–23 (2007)
  * (3) Bai, Z.Z., Golub, G.H., Ng, M.K.: Hermitian and skew-Hermitian splitting methods for non-Hermitian positive definite linear systems.  SIAM Journal on Matrix Analysis and Applications 24(3), 603–626 (2003)
  * (4) Bai, Z.Z., Golub, G.H., Ng, M.K.: On successive-overrelaxation acceleration of the Hermitian and skew-Hermitian splitting iterations.  Numerical Linear Algebra with Applications 14(4), 319–335 (2007)
  * (5) Bai, Z.Z., Golub, G.H., Ng, M.K.: On inexact Hermitian and skew-Hermitian splitting methods for non-Hermitian positive definite linear systems.  Linear Algebra and Its Applications 428(2-3), 413–440 (2008)
  * (6) Bai, Z.Z., Parlett, B.N., Wang, Z.Q.: On generalized successive overrelaxation methods for augmented linear systems.  Numerische Mathematik 102(1), 1–38 (2005)
  * (7) Benzi, M., Golub, G.H.: A preconditioner for generalized saddle point problems.  SIAM Journal on Matrix Analysis and Applications 26(1), 20–41 (2004)
  * (8) Benzi, M., Golub, G.H., Liesen, J.: Numerical solution of saddle point problems.  Acta Numer. 14, 1–137 (2005).  DOI 10.1017/S0962492904000212
  * (9) Browder, F.E., Petryshyn, W.V.: Construction of fixed points of nonlinear mappings in Hilbert space.  Journal of Mathematical Analysis and Applications 20(2), 197–228 (1967)
  * (10) Chambolle, A., Pock, T.: A first-order primal-dual algorithm for convex problems with applications to imaging.  Journal of mathematical imaging and vision 40(1), 120–145 (2011)
  * (11) Chambolle, A., Pock, T.: On the ergodic convergence rates of a first-order primal–dual algorithm.  Mathematical Programming 159(1), 253–287 (2016)
  * (12) Chen, G., Teboulle, M.: Convergence analysis of a proximal-like minimization algorithm using Bregman functions.  SIAM Journal on Optimization 3(3), 538–543 (1993)
  * (13) Chen, L., Luo, H.: A unified convergence analysis of first order convex optimization methods via strong Lyapunov functions.  arXiv preprint arXiv:2108.00132 (2021)
  * (14) Chen, L., Wei, J.: Transformed primal-dual methods for nonlinear saddle point systems.  arXiv preprint arXiv:2208.02444 (2022)
  * (15) Gabay, D., Mercier, B.: A dual algorithm for the solution of nonlinear variational problems via finite element approximation.  Computers & mathematics with applications 2(1), 17–40 (1976)
  * (16) Glowinski, R., Marroco, A.: Sur l’approximation, par éléments finis d’ordre un, et la résolution, par pénalisation-dualité d’une classe de problèmes de Dirichlet non linéaires.  Revue française d’automatique, informatique, recherche opérationnelle. Analyse numérique 9(R2), 41–76 (1975)
  * (17) Güler, O.: New proximal point algorithms for convex minimization.  SIAM Journal on Optimization 2(4), 649–664 (1992)
  * (18) Haddad, W.M., Chellaboina, V.: Nonlinear Dynamical Systems and Control: A Lyapunov-Based Approach.  Princeton University Press (2008).  URL http://www.jstor.org/stable/j.ctvcm4hws
  * (19) Hadjidimos, A.: Accelerated overrelaxation method.  Mathematics of Computation 32(141), 149–157 (1978)
  * (20) Hadjidimos, A.: Successive overrelaxation (SOR) and related methods.  Journal of Computational and Applied Mathematics 123(1-2), 177–199 (2000)
  * (21) He, B., Yuan, X.: Convergence analysis of primal-dual algorithms for a saddle-point problem: from contraction perspective.  SIAM Journal on Imaging Sciences 5(1), 119–149 (2012)
  * (22) Hestenes, M.R.: Multiplier and gradient methods.  Journal of optimization theory and applications 4(5), 303–320 (1969)
  * (23) Idema, R., Vuik, C.: A minimal residual method for shifted skew-symmetric systems.  Delft University of Technology (2007)
  * (24) Jiang, E.: Algorithm for solving shifted skew-symmetric linear system.  Frontiers of Mathematics in China 2(2), 227 (2007)
  * (25) Khalil, H.K.: Nonlinear systems; 3rd ed.  Prentice-Hall, Upper Saddle River, NJ (2002)
  * (26) Kim, D.: Accelerated proximal point method for maximally monotone operators.  Mathematical Programming 190(1), 57–87 (2021)
  * (27) Kovalev, D., Gasnikov, A., Richtárik, P.: Accelerated primal-dual gradient method for smooth and convex-concave saddle-point problems with bilinear coupling.  arXiv preprint arXiv:2112.15199 (2021)
  * (28) Liu, F., Nashed, M.Z.: Regularization of nonlinear ill-posed variational inequalities and convergence rates.  Set-Valued Analysis 6(4), 313–344 (1998)
  * (29) Luo, H., Chen, L.: From differential equation solvers to accelerated first-order methods for convex optimization.  Mathematical Programming (2021).  DOI 10.1007/s10107-021-01713-3
  * (30) Nesterov, Y.: Introductory lectures on convex optimization: A basic course, vol. 87.  Springer Science & Business Media (2003)
  * (31) Nesterov, Y.E.: A method of solving a convex programming problem with convergence rate ${O}\bigl{(}\frac{1}{k^{2}}\bigr{)}$.  In: Doklady Akademii Nauk, vol. 269, pp. 543–547. Russian Academy of Sciences (1983)
  * (32) Polyak, B.T.: Some methods of speeding up the convergence of iteration methods.  Ussr computational mathematics and mathematical physics 4(5), 1–17 (1964)
  * (33) Powell, M.J.: A method for nonlinear constraints in minimization problems.  Optimization pp. 283–298 (1969)
  * (34) Rockafellar, R.T.: Monotone operators and the proximal point algorithm.  SIAM journal on control and optimization 14(5), 877–898 (1976)
  * (35) Salim, A., Condat, L., Kovalev, D., Richtárik, P.: An optimal algorithm for strongly convex minimization under affine constraints.  In: International Conference on Artificial Intelligence and Statistics, pp. 4482–4498. PMLR (2022)
  * (36) Thekumparampil, K.K., He, N., Oh, S.: Lifted primal-dual method for bilinearly coupled smooth minimax optimization.  In: International Conference on Artificial Intelligence and Statistics, pp. 4281–4308. PMLR (2022)
  * (37) Zhang, J., Hong, M., Zhang, S.: On lower iteration complexity bounds for the convex concave saddle point problems.  Mathematical Programming 194(1), 901–935 (2022)
  * (38) Zhu, D.L., Marcotte, P.: Co-coercivity and its role in the convergence of iterative schemes for solving variational inequalities.  SIAM Journal on Optimization 6(3), 714–726 (1996)
  * (39) Zhu, M., Chan, T.: An efficient primal-dual hybrid gradient algorithm for total variation image restoration.  Ucla Cam Report 34, 8–34 (2008)