# Sublinear Convergence Rates of Extragradient-Type Methods: A Survey on
Classical and Recent Developments

Quoc Tran-Dinh

Department of Statistics and Operations Research

The University of North Carolina at Chapel Hill
318 Hanes Hall, UNC-Chapel Hill, NC 27599-3260.

Email<EMAIL_ADDRESS>
(Version 1 (March 30, 2023))

###### Abstract

The extragradient (EG), introduced by G. M. Korpelevich in 1976, is a well-
known method to approximate solutions of saddle-point problems and their
extensions such as variational inequalities and monotone inclusions. Over the
years, numerous variants of EG have been proposed and studied in the
literature. Recently, these methods have gained popularity due to new
applications in machine learning and robust optimization. In this work, we
survey the latest developments in the EG method and its variants for
approximating solutions of nonlinear equations and inclusions, with a focus on
the monotonicity and co-hypomonotonicity settings. We provide a unified
convergence analysis for different classes of algorithms, with an emphasis on
sublinear best-iterate and last-iterate convergence rates. We also discuss
recent accelerated variants of EG based on both Halpern fixed-point iteration
and Nesterov’s accelerated techniques. Our approach uses simple arguments and
basic mathematical tools to make the proofs as elementary as possible, while
maintaining generality to cover a broad range of problems.

##  1 Introduction

The generalized equation (also called the [non]linear inclusion) provides a
unified template to model various problems in computational mathematics and
related fields such as the optimality condition of optimization problems (in
both unconstrained and constrained settings), minimax optimization,
variational inequality, complementarity, two-person game, and fixed-point
problems, see, e.g., [11, 24, 50, 112, 116, 118, 120]. Theory and numerical
methods for this equation and its special cases have been extensively studied
for many decades, see, e.g., the following monographs and the references
quoted therein [11, 50, 94, 119]. At the same time, several applications of
this mathematical tool in operations research, economics, uncertainty
quantification, and transportations have been investigated [14, 52, 61, 50,
72]. In the last few years, there has been a surge of research in minimax
problems due to new applications in machine learning and robust optimization,
especially in generative adversarial networks (GANs), adversarial training,
and distributionally robust optimization, see, e.g., [4, 14, 55, 76, 84, 114]
as a few examples. Minimax problems have also found new applications in online
learning and reinforcement learning, among many others, see, e.g., [4, 9, 15,
55, 67, 76, 78, 84, 114, 139]. Such prominent applications have motivated the
research in minimax optimization and variational inequality problems (VIPs).
On the one hand, classical algorithms such as gradient descent-ascent,
extragradient, and primal-dual methods have been revisited, improved, and
extended. On the other hand, new variants such as accelerated extragradient
and accelerated operator splitting schemes have also been developed and
equipped with rigorous convergence guarantees and practical performance
evaluation. This new development motivates us to write this survey paper, with
the focus on sublinear convergence rate analysis.

Problem statements. Since there is a vast amount of literature on the
generalized equation, we will only present the recent developments on
sublinear convergence rates of the extragradient (EG) method and its variants
for approximating the solutions of the following generalized equation (also
known as a [composite] nonlinear inclusion) and its specific cases:

$\textrm{Find $x^{\star}\in\mathrm{dom}(\Phi)$ such that:}\quad 0\in\Phi
x^{\star}\equiv Fx^{\star}+Tx^{\star},$ (NI)

where $F:\mathbb{R}^{p}\to\mathbb{R}^{p}$ is a single-valued operator,
$T:\mathbb{R}^{p}\rightrightarrows 2^{\mathbb{R}^{p}}$ is a set-valued (or
multivalued) mapping from $\mathbb{R}^{p}$ to $2^{\mathbb{R}^{p}}$ (the set of
all subsets of $\mathbb{R}^{p}$), $\Phi:=F+T$, and
$\mathrm{dom}(\Phi):=\mathrm{dom}(F)\cap\mathrm{dom}(T)$ is the domain of
$\Phi$, which is the intersection of the domains of $F$ and $T$. In this
paper, we focus on the finite-dimensional Euclidean spaces $\mathbb{R}^{p}$
and $\mathbb{R}^{n}$ for ease of presentation. However, it is worth noting
that most of the results presented in this paper can be extended to Hilbert
spaces, as demonstrated in the existing literature.

Special cases. If $F=0$, then (NI) reduces to a generalized equation or a
[non]linear inclusion $0\in Tx^{\star}$. Alternatively, if $T=0$, then (NI)
reduces to a [non]linear equation:

$\textrm{Find $x^{\star}\in\mathrm{dom}(F)$ such that:}\quad Fx^{\star}=0.$
(NE)

If $T:=\partial{g}$, the subdifferential of a proper, closed, and convex
function $g:\mathbb{R}^{p}\to\mathbb{R}\cup\\{+\infty\\}$, then (NI) reduces a
mixed variational inequality problem (MVIP):

$\textrm{Find $x^{\star}\in\mathrm{dom}(\Phi)$ such that:}\quad\langle
Fx^{\star},x-x^{\star}\rangle+g(x)-g(x^{\star})\geq 0,\ \textrm{for all}\
x\in\mathrm{dom}(\Phi).$ (MVIP)

In particular, if $T=\mathcal{N}_{\mathcal{X}}$, the normal cone of a
nonempty, closed, and convex set $\mathcal{X}$ in $\mathbb{R}^{p}$ (i.e.
$g=\delta_{\mathcal{X}}$, the indicator of $\mathcal{X}$), then (MVIP) reduces
the classical (Stampacchia) variational inequality problem (VIP):

$\textrm{Find $x^{\star}\in\mathcal{X}$ such that:}\quad\langle
Fx^{\star},x-x^{\star}\rangle\geq 0,\ \textrm{for all}\ x\in\mathcal{X}.$
(VIP)

While (VIP) can be viewed as a primal VIP (or a strong VIP), its dual (or
weak) form can be written as

$\textrm{Find $x^{\star}\in\mathcal{X}$ such that:}\quad\langle
Fx,x-x^{\star}\rangle\geq 0,\ \textrm{for all}\ x\in\mathcal{X},$ (DVIP)

which is known as Minty’s variational inequality problem. If $F$ is monotone
(see the definition in Section 2) then both problems (VIP) and (MVIP) are
equivalent, i.e. their solution sets are identical, [50, 72]. One important
special case of (NI) or (VIP) is the optimality condition of minimax problems
of the form:

$\min_{u\in\mathbb{R}^{m}}\max_{v\in\mathbb{R}^{n}}\Big{\\{}\mathcal{L}(u,v):=\varphi(u)+\mathcal{H}(u,v)-\psi(v)\Big{\\}},$
(1)

where $\varphi:\mathbb{R}^{m}\to\mathbb{R}\cup\\{+\infty\\}$ and
$\psi:\mathbb{R}^{n}\to\mathbb{R}\cup\\{+\infty\\}$ are often proper, closed,
and convex functions, and
$\mathcal{H}:\mathbb{R}^{m}\times\mathbb{R}^{n}\to\mathbb{R}$ is a bifunction,
often assumed to be differentiable, but not necessarily convex-concave. If we
denote $x:=[u,v]$ as the concatenation of $x$ and $y$, and define
$T:=[\partial{\varphi},\partial{\psi}]$ and
$F:=[\nabla_{u}{\mathcal{H}}(u,v),-\nabla_{v}{\mathcal{H}}(u,v)]$, then the
optimality condition of (1) is exactly captured by (NI).

Related work. Extensive research has been conducted in the literature to
investigate the existence of solutions and theoretical properties of (NI) and
its special cases. This research has been conducted under various assumptions
of monotonicity and extensions, including quasi-monotone, pseudo-monotone, and
weakly monotone notions. Relevant literature references on this topic include
[11, 50, 72, 99, 144]. Moreover, solution methods for (NI) and its special
cases have been well-developed, particularly in the context of monotonicity
and related extensions such as quasi-monotone, pseudo-monotone, or star-
monotone notion. In addition, nonmonotone instances of (NI) have also received
extensive attention in the literature, with many theoretical results and
algorithms focusing on local properties. Additional information can be found
in references such as [12, 18, 17, 35, 44, 109, 115, 119].

Existing solution methods for (NI) and its special cases often rely on a
fundamental assumption: maximal monotonicity of $F$ and $T$, or of $\Phi$ to
guarantee global convergence. These methods generally generalize existing
optimization algorithms such as gradient, proximal-point, Newton, and
interior-point schemes to (NI) and its special cases [35, 50, 51, 92, 107,
118, 133, 137], while leveraging the splitting structure of (NI) to use
individual operators defined on $F$ and $T$. This approach leads to a class of
splitting algorithms for solving (NI) such as forward-backward splitting (FBS)
and Douglas-Rachford (DRS) splitting schemes, as seen in [11, 36, 41, 46, 80,
81]. Alternatively, other approaches rely on primal-dual, dual averaging, and
mirror descent techniques, with notable works including [31, 100, 104]. These
methods have also been further studied in many recent works such as [32, 33,
37, 39, 49, 63, 105, 131, 134, 145].

When it comes to convergence analysis for gradient-based/forward methods,
there is a fundamental challenge for generalized equation (NI) because an
objective function, which plays a central role in guaranteeing convergence for
optimization problems, does not exist. This creates a significant challenge,
particularly in nonmonotone settings. Additionally, unlike convex functions
where strong properties such as coerciveness and cyclic monotonicity hold for
their [sub]gradients beyond monotonicity, this is not the case for general
monotone and Lipschitz continuous operators. This lack of a strong property
results in gradient-based (or forward) methods being non-convergent, which
limits their practicality, see, e.g., [50]. To address this issue, the
extragradient (EG) method was introduced by G. M. Korpelevich in 1976 [74] and
also by A. S. Antipin in [3]. This method performs two sequential gradient
steps at each iteration, making it twice as expensive as the standard gradient
method, but converges under only the monotonicity and the Lipschitz continuity
of $F$. Since then, this method has been extended and modified in different
directions to reduce its per-iteration complexity, including in certain
nonmonotone settings, see, e.g., [1, 28, 29, 66, 69, 87, 89, 90, 97, 98, 113,
122, 123, 124, 135, 136]. Among these variants of EG, the past-extragradient
scheme in [113] and Teng’s forward-backward-forward splitting method in [136]
are the most notable ones. However, the results discussed here are only
applicable to the monotone setting of (NI) and its special cases.
Additionally, most of the convergence results discussed are asymptotic,
leading to sublinear “best-iterate” convergence rates of the residual norm
associated with (NI). Under stronger assumptions such as “strong
monotonicity”, linear convergence rates can be achieved. Such types of
convergence guarantees have been widely studied in the literature and are
beyond the scope of this paper, see, e.g., [11, 50, 72].

Motivated by recent applications in machine learning and robust optimization,
such as Generative Adversarial Networks (GANs), adversarial training,
distributionally robust optimization, reinforcement learning, and online
learning, several methods for solving minimax problems have become critically
important and attractive. This is particularly true in nonconvex-nonconcave,
large-scale, and stochastic settings, as evidenced in works such as [4, 9, 14,
15, 55, 67, 76, 78, 84, 114]. Several researchers have proposed and revisited
EG and its variants, including [16, 38, 43, 111]. A notable work is due to
[43], where the authors proposed an EG-plus (EG+) variant of EG, capable of
handling nonmonotone instances of (NE), known as weak-Minty solutions. In
[111], this method was further extended to (NI), while [16, 83] modified EG+
for Popov’s methods, as well as optimistic gradient variants.

In contrast to classical methods, there has been a significant focus on
developing accelerated methods for solving (NI) and its special cases under
both monotone and co-hypomonotone structures. Early works in this area relied
on dual averaging and mirror descent techniques such as those proposed in [37,
100, 104], which require the monotonicity or specific assumptions. Attouch et
al [8] proposed accelerated proximal-point methods for solving (NI) under the
maximal monotonicity of $\Phi$. Since then, numerous works have followed up
and explored Nesterov’s acceleration-type methods guided by dynamical systems,
utilizing momentum and correction terms for solving (NI) under monotone
assumptions, as demonstrated in works such as [5, 21, 22, 70, 86, 85].
Accelerated methods based on Halpern’s fixed-point iteration [60] have also
gained popularity. Although initially developed to approximate a fixed-point
of a nonexpansive operator, this method can be applied to solve (NE), (VIP),
and (NI) under monotonicity. In [77], it was shown that Halpern’s fixed-point
iteration can achieve $\mathcal{O}\big{(}1/k\big{)}$ last-iterate convergence
rates using a specific choice of parameters, where $k$ is the iteration
counter. The authors in [42] further exploited this approach to solve monotone
VIPs of the form (VIP). Yoon and Ryu extended Halpern’s fixed-point iteration
idea to EG methods to solve (NE) without the co-coerciveness assumption on $F$
in their pioneering work [143]. Lee and Kim [75] proposed a similar algorithm
for solving (NE) under the co-hypomonotonicity, further advancing [143]
without sacrificing the $\mathcal{O}\big{(}1/k\big{)}$-convergence rates. In
[132], the authors proposed a Halpern-type variant for the past-extragradient
method in [113] by adopting the technique from [143]. Recently, [25, 27]
extended [143] and [132] to (VIP) and (NI) under either monotonicity or co-
hypomonotonicity assumptions. New convergence analysis for these schemes can
also be found in [130]. Note that both Halpern’s fixed-point iteration and
Nesterov’s accelerated schemes for solving (NE) and (NI) are related to each
other, as shown in [129] for different methods, including EG. Nesterov’s
accelerated variants of EG can also be found in [130].

What does this paper survey? Our main objective is to provide a comprehensive
survey of both classical and recent sublinear convergence rate results for EG
and its variants for solving (NI) and its special cases, as summarized in
Table 1. Specifically, we survey the following results.

Table 1: Summary of the results surveyed in this paper and the most related
references

Methods | Assumptions | Add. Assumptions | Convergence Rates | Citations
---|---|---|---|---
For solving (NE)
EG/EG+/FBFS | wMs |  $F$ is chm |  $\mathcal{O}\big{(}1/\sqrt{k}\big{)}$ best and last | [43, 54, 56, 83]
PEG/OG/FRBS/RFBS/GR | wMs |  $F$ is chm |  $\mathcal{O}\big{(}1/\sqrt{k}\big{)}$ best and last | [16, 83]
EAG/FEG/AEG |  $F$ is chm | None |  $\mathcal{O}\big{(}1/k\big{)}$ last-iterate | [143, 70, 129]
PEAG/APEG |  $F$ is chm | None |  $\mathcal{O}\big{(}1/k\big{)}$ last-iterate | [132, 129]
For solving (NI), (MVIP), and (VIP)
EG/EG+ |  $\Phi$ is wMs |  $F$ is mono, $T$ is 3-cm |  $\mathcal{O}\big{(}1/\sqrt{k}\big{)}$ best and last | [26]
FBFS |  $\Phi$ is wMs | None |  $\mathcal{O}\big{(}1/\sqrt{k}\big{)}$ best-iterate | [111, 83]
OG/FRFS |  $\Phi$ is wMs | None |  $\mathcal{O}\big{(}1/\sqrt{k}\big{)}$ best-iterate | [83]
RFBS |  $F$ is mono |  $T$ is mono |  $\mathcal{O}\big{(}1/\sqrt{k}\big{)}$ best and last | [30, 87]
GR |  $F$ is mono |  $T$ is $3$-cm |  $\mathcal{O}\big{(}1/\sqrt{k}\big{)}$ best-iterate | [88]
EAG/FEG/AEG |  $\Phi$ is chm | None |  $\mathcal{O}\big{(}1/k\big{)}$ last-iterate | [25, 130]
PEAG/APEG |  $\Phi$ is chm | None |  $\mathcal{O}\big{(}1/k\big{)}$ last-iterate | [27, 130]

Abbreviations: EG $=$ extragradient; PEG $=$ past extragradient; FBFS $=$
forward-backward-forward splitting; OG $=$ optimistic gradient; FRBS $=$
forward-reflected-backward splitting; RFBS $=$ reflected-forward-backward
splitting; GR $=$ golden ratio; EAG $=$ extra-anchored gradient; FEG $=$ fast
extragradient; PEAG $=$ past extra-anchored gradient; AEG $=$ Nesterov’s
accelerated extragradient; and APEG $=$ Nesterov’s accelerated past
extragradient. In addition, wMs $=$ weak-Minty solution; mono $=$ monotone;
chm $=$ co-hypomonotone; and 3-cm $=$ $3$-cyclically monotone.

  * •

First, we present both the $\mathcal{O}\big{(}1/\sqrt{k}\big{)}$-best-iterate
and last-iterate sublinear convergence rate results of EG and its variants for
solving (NE). The best-iterate rate is classical for the monotone case, but
has recently been obtained under a weak Minty solution condition, see [43,
111] for EG and [16, 83] for past-EG in the non-composite case, i.e., for
solving (NE). The last-iterate convergence rates for EG and past-EG have been
recently proven in [54, 56] for the monotone equation (NE) and in [83] for the
co-hypomonotone case (see also [57]). In this paper, we provide a new and
unified proof that covers the results in [54, 56, 83]. Our results are stated
in a single theorem.

  * •

Second, we review the $\mathcal{O}\big{(}1/\sqrt{k}\big{)}$-sublinear best-
iterate convergence rates for EG and past-EG (also known as Popov’s method) to
solve (NI) under the monotonicity of $\Phi$. We unify the proof of both
methods in a single theorem and extend it to cover monotone inclusions of the
form (NI) instead of VIP or MVIP as in the literature. We also prove
$\mathcal{O}\big{(}1/\sqrt{k}\big{)}$ last-iterate convergence for the class
of EG-type schemes for solving (NI) under the monotonicity of $F$ and the
$3$-cyclical monotonicity of $T$ (in particular, for solving (MVIP)), which
covers the results in [26] as special cases. Next, we discuss the
$\mathcal{O}\big{(}1/\sqrt{k}\big{)}$-sublinear best-iterate convergence rates
of the FBFS scheme and its variant: optimistic gradient under the weak-Minty
solution notion, which was obtained in [83]. We again unify the proof in a
single theorem, and our analysis is also different from [83].

  * •

Third, we provide a new convergence analysis for both best-iterate and last-
iterate rates of the reflected forward-backward splitting (RFBS) methods for
solving (NI) under the monotonicity of $F$ and $T$. RFBS was proposed in [87]
to solve (VIP) and was extended to solve (NI) in [30]. The best-iterate rates
were proven in these works, and the last-iterate rate of RFBS for solving
(VIP) has recently been proven in [27]. Our result here is more general and
covers these works as special cases. In addition, we also review the best-
iterate convergence rate of the golden ration method in [88], but extend it to
the case $T$ is $3$-cyclically monotone, and extend the range of the golden-
ratio parameter $\omega$ to $1<\omega<1+\sqrt{3}$ instead of fixing it at
$\omega:=\frac{1+\sqrt{5}}{2}$ as in [88].

  * •

Fourth, we present a new analysis for the extra-anchored gradient (EAG) method
to solve (NI), which covers the results in [25, 143] as special cases. Our
result extends to $3$-cyclically monotone operator $T$.

  * •

Fifth, we summarize the convergence results of EAG and past EAG (also called
fast extragradient method [70]) for solving (NI) under both monotonicity and
co-hypomonotonicity of $\Phi$ from [130]. Note that EAG and past-EAG were
first proposed in [143] and [132], respectively, to solve monotone (NE). EAG
was extended to the co-hypomonotone case of (NE) in [70]. Recently, [25, 27]
extended EAG and past-EAG to solve (NI) under the co-hypomonotonicity of
$\Phi$.

  * •

Finally, we review two Nesterov’s accelerated extragradient methods presented
in [129, 130] for solving (NI) under the co-hypomonotonicity of $\Phi$, which
achieve the same last-iterate convergence rates as EAG. Note that Nesterov’s
accelerated extragradient methods have recently been studied in [21] for
solving (NE) via a dynamical system point of view.

What is not covered in this paper? The literature on EG and its variants is
extensive, and it is not feasible for us to cover it in detail in this paper.
First, there are various classical and recent variants of EG and past-EG, such
as those discussed in [28, 29, 66, 69, 89, 90, 123, 124, 122, 138], that are
not included in this paper. These methods are essentially rooted from EG with
the aim of improving the per-iteration complexity, theoretical aspects, or
practical performance. Second, we do not review results from methods such as
gradient/forward, forward-backward splitting, proximal-point and its variants,
inertial, dual averaging, mirror descent, and projective methods. The majority
of these methods is not immediately derived from EG, including recent
developments such as those in [19, 20, 23, 35, 34, 33, 47, 48, 79]. Third, we
do not cover stochastic and randomized methods, including recent works such as
[[2, 40, 53, 64, 68, 65, 110, 108, 128]. Fourth, we do not present adaptive
stepsizes/parameters and linesearch variants of EG-type methods. Fifth, we do
not disuss continuous view of EG-type methods via dynamical systems or
ordinary differential equations (ODEs), which is an emerging research topic in
recent years. Finally, we also do not cover specific applications to minimax
problems and other concrete applications.

Paper outline. This paper is organized as follows. Section 2 reviews basic
concepts and related results used in this paper. Section 3 covers the
convergence rate results of EG and its variants for solving (NE). Section 4
discusses the convergence rate results of EG and past-EG for solving (NI).
Section 5 provides a new convergence rate analysis of FBFS and OG for solving
(NI). Section 6 presents a new analysis for both the reflected-forward-
backward splitting and golden ratio methods for solving (NI). Section 7
focuses on the extra-anchored gradient method and its variants for solving
(NI). Finally, Section 8 presents Nesterov’s accelerated variants of EG for
solving (NI). We conclude this paper with some final remarks.

##  2 Background and Preliminary Results

To prepare for our survey, we will briefly review certain basic concepts and
properties of monotone operators and their extensions, as well as resolvents
and other related mathematical tools. These concepts and properties are well-
known and can be found in several monographs, including [11, 24, 50, 112, 116,
117, 118, 120].

###  2.1 Basic concepts, monotonicity, and Lipschitz continuity

We work with finite dimensional Euclidean spaces $\mathbb{R}^{p}$ and
$\mathbb{R}^{n}$ equipped with standard inner product
$\langle\cdot,\cdot\rangle$ and Euclidean norm $\|\cdot\|$. For a set-valued
or multivalued mapping $T:\mathbb{R}^{p}\rightrightarrows 2^{\mathbb{R}^{p}}$,
$\mathrm{dom}(T)=\left\\{x\in\mathbb{R}^{p}:Tx\not=\emptyset\right\\}$ denotes
its domain, $\mathrm{ran}(T):=\bigcup_{x\in\mathrm{dom}(T)}Tx$ is its range,
and $\mathrm{gra}(T)=\left\\{(x,y)\in\mathbb{R}^{p}\times\mathbb{R}^{p}:y\in
Tx\right\\}$ stands for its graph, where $2^{\mathbb{R}^{p}}$ is the set of
all subsets of $\mathbb{R}^{p}$. The inverse of $T$ is defined as
$T^{-1}y:=\\{x\in\mathbb{R}^{p}:y\in Tx\\}$. For a proper, closed, and convex
function $f:\mathbb{R}^{p}\to\mathbb{R}\cup\\{+\infty\\}$,
$\mathrm{dom}(f):=\\{x\in\mathbb{R}^{p}:f(x)<+\infty\\}$ denotes the domain of
$f$, $\partial{f}$ denotes the subdifferential of $f$, and $\nabla{f}$ stands
for the gradient of $f$. For any function $f$, which can be nonconvex, we call
$f^{*}(y):=\sup_{x\in\mathbb{R}^{p}}\\{\langle x,y\rangle-f(x)\\}$ the Fenchel
conjugate of $f$.

(a) Monotonicity. For a single-valued or multivalued mapping
$T:\mathbb{R}^{p}\rightrightarrows 2^{\mathbb{R}^{p}}$ and $\mu\in\mathbb{R}$,
we say that $T$ is $\mu$-monotone if $\langle
u-v,x-y\rangle\geq\mu\|x-y\|^{2}$ for all $(x,u),(y,v)\in\mathrm{gra}(T)$. If
$T$ is single-valued, then this condition reduces to $\langle Tx-
Ty,x-y\rangle\geq\mu\|x-y\|^{2}$ for all $x,y\in\mathrm{dom}(T)$. If $\mu=0$,
then we say that $T$ is monotone. If $\mu>0$, then $T$ is $\mu$-strongly
monotone (or sometimes called coercive), where $\mu>0$ is called a strong
monotonicity parameter. If $\mu<0$, then we say that $T$ is weakly monotone.
It is also called $-\mu$-hypomonotone, see [12]. If $T=\partial{g}$, the
subdifferential of a proper and convex function, then $T$ is also monotone. If
$g$ is $\mu$-strongly convex with $\mu>0$, then $T=\partial{g}$ is also
$\mu$-strongly monotone.

Alternatively, if there exists $\rho\in\mathbb{R}$ such that $\langle
u-v,x-y\rangle\geq\rho\|u-v\|^{2}$ for all $(x,u),(y,v)\in\mathrm{gra}(T)$,
then we say that $T$ is $\rho$-comonotone. If $\rho=0$, then this condition
reduces to the monotonicity of $T$. If $\rho>0$, then $T$ is called
$\rho$-cocoercive. In particular, if $\rho=1$, then $T$ is firmly
nonexpansive. If $\rho<0$, then $T$ is called $-\rho$-cohypomonotone, see,
e.g., [12, 35]. For a mapping $T$, we say that $T$ is pseudo-monotone if
$\langle u,y-x\rangle\geq 0$ implies $\langle v,y-x\rangle\geq 0$ for all
$(x,u),(y,v)\in\mathrm{gra}(T)$. Clearly, if $T$ is monotone, then it is also
pseudo-monotone, but the conversion is not true in general.

We say that $T$ is maximally $\mu$-monotone if $\mathrm{gra}(T)$ is not
properly contained in the graph of any other $\mu$-monotone operator. If
$\mu=0$, then we say that $T$ is maximally monotone. Note that $T$ is
maximally monotone, then $\eta T$ is also maximally monotone for any $\eta>0$,
and if $T$ and $U$ are maximally monotone, and
$\mathrm{dom}(T)\cap\mathrm{int}\left(\mathrm{dom}(U)\right)\not=\emptyset$,
then $T+U$ is maximally monotone. For a proper, closed, and convex function
$f:\mathbb{R}^{p}\to\mathbb{R}\cup\\{+\infty\\}$, the subdifferential
$\partial{f}$ of $f$ is maximally monotone.

For a given mapping $T$ such that
$\mathrm{zer}(T):=\left\\{x\in\mathrm{dom}(T):0\in Tx\right\\}\neq\emptyset$,
we say that $T$ is star-monotone (respectively, $\mu$-star-monotone or
$\rho$-star-comonotone (see [82])) if for some $x^{\star}\in\mathrm{zer}(T)$,
we have $\langle u,x-x^{\star}\rangle\geq 0$ (respectively, $\langle
u,x-x^{\star}\rangle\geq\mu\|x-x^{\star}\|^{2}$ or $\langle
u,x-x^{\star}\rangle\geq\rho\|u\|^{2}$) for all $(x,u)\in\mathrm{gra}(T)$.
Clearly, if $T$ is monotone (respectively, $\mu$-monotone or $\rho$-co-
monotone), then it is also star-monotone (respectively, $\mu$-star-monotone or
$\rho$-star-comonotone). However, the reverse statement does not hold in
general.

(b) Cyclic monotonicity. We also say that a mapping $T$ is $m$-cyclically
monotone ($m\geq 2$) if $\sum_{i=1}^{n}\langle u^{i},x^{i}-x^{i+1}\rangle\geq
0$ for all $(x^{i},u^{i})\in\mathrm{gra}(T)$ and $x_{1}=x_{m+1}$ (see [11]).
We say that $T$ is cyclically monotone if it is $m$-cyclically monotone for
every $m\geq 2$. If $T$ is $m$-cyclically monotone, then it is also
$\hat{m}$-cyclically monotone for any $2\leq\hat{m}\leq m$. Since a
$2$-cyclically monotone operator $T$ is monotone, any $m$-cyclically monotone
operator $T$ is $2$-cyclically monotone, and thus is also monotone. An
$m$-cyclically monotone operator $T$ is called maximally $m$-cyclically
monotone if $\mathrm{gra}(T)$ is not properly contained into the graph of any
other $m$-cyclically monotone operator. As proven in [11, Theorem 22.18] that
$T$ is maximally cyclically monotone iff $T=\partial{f}$, the subdifferential
of a proper, closed, and convex function $f$. However, there exist maximally
$m$-cyclically monotone operators (e.g., rotation linear operators) that are
not the subdifferential $\partial{f}$ of a proper, closed, and convex function
$f$, see, e.g., [11]. Furthermore, as indicated in [10, Example 2.16], there
exist maximally $3$-cyclically monotone operators that are not maximal
monotone.

(c) Lipschitz continuity and contraction. A single-valued or multivalued
mapping $T$ is said to be $L$-Lipschitz continuous if
$\sup\left\\{\|u-v\|:u\in Tx,\ v\in Ty\right\\}\leq L\|x-y\|$ for all
$x,y\in\mathrm{dom}(T)$, where $L\geq 0$ is a Lipschitz constant. If $T$ is
single-valued, then this condition becomes $\|Tx-Ty\|\leq L\|x-y\|$ for all
$x,y\in\mathrm{dom}(T)$. If $L=1$, then we say that $T$ is nonexpansive, while
if $L\in[0,1)$, then we say that $T$ is $L$-contractive, and $L$ is its
contraction factor. If $T$ is $\rho$-co-coercive with $\rho>0$, then $T$ is
also $L$-Lipschitz continuous with the Lipschitz constant $L:=\frac{1}{\rho}$.
However, the reverse statement is not true in general. For a continuously
differentiable function $f:\mathbb{R}^{p}\to\mathbb{R}$, we say that $f$ is
$L$-smooth if its gradient $\nabla{f}$ is $L$-Lipschitz continuous on
$\mathrm{dom}(f)$. If $f$ is convex and $L$-smooth, then $\nabla{f}$ is
$\frac{1}{L}$-co-coercive and vice versa, see, e.g., [102].

(d) Normal cone. Given a nonempty, closed, and convex set $\mathcal{X}$ in
$\mathbb{R}^{p}$, the normal cone of $\mathcal{X}$ is defined as
$\mathcal{N}_{\mathcal{X}}(x):=\\{w\in\mathbb{R}^{p}:\langle w,x-y\rangle\geq
0,\ \forall y\in\mathcal{X}\\}$ if $x\in\mathcal{X}$ and
$\mathcal{N}_{\mathcal{X}}(x)=\emptyset$, otherwise. The dual cone
$\mathcal{N}^{*}_{\mathcal{X}}(x):=\\{w\in\mathbb{R}^{p}:\langle
w,u\rangle\geq 0,\ \forall u\in\mathcal{N}_{\mathcal{X}}(x)\\}$ of
$\mathcal{N}_{\mathcal{X}}(x)$ at $x$ is $\mathcal{T}_{\mathcal{X}}(x)$, the
tangent cone of $\mathcal{X}$ at $x$. If $f:=\delta_{\mathcal{X}}$, the
indicator of $\mathcal{X}$, and $f^{*}$ is its Fenchel conjugate, then
$\partial{f}=\mathcal{N}_{\mathcal{X}}$, and
$\partial{f^{*}}=\mathcal{T}_{\mathcal{X}}$.

(e) Resolvent and proximal operators. The operator
$J_{T}x:=\left\\{y\in\mathbb{R}^{p}:x\in y+Ty\right\\}$ is called the
resolvent of $T$, denoted by $J_{T}x=(\mathbb{I}+T)^{-1}x$, where $\mathbb{I}$
is the identity mapping. If $T$ is $\rho$-monotone with $\rho>-1$, then
evaluating $J_{T}$ requires solving a strongly monotone inclusion $0\in
y-x+Ty$. If $T=\partial{f}$, the subdifferential of proper, closed, and convex
function $f$, then $J_{T}x$ reduces to the proximal operator of $f$, denoted
by $\mathrm{prox}_{f}$, which can be computed as
$\mathrm{prox}_{f}(x):=\mathrm{arg}\min_{y}\\{f(y)+(1/2)\|y-x\|^{2}\\}$. In
particular, if $T=\mathcal{N}_{\mathcal{X}}$, the normal cone of a closed and
convex set $\mathcal{X}$, then $J_{T}$ is the projection onto $\mathcal{X}$,
denoted by $\mathrm{proj}_{\mathcal{X}}$. If $T$ is maximally monotone, then
$\mathrm{ran}(\mathbb{I}+T)=\mathbb{R}^{p}$ (by Minty’s theorem) and $T$ is
firmly nonexpansive (and thus nonexpansive).

###  2.2 Best-iterate and last-iterate convergence rates

The results presented in this paper are related to two types of sublinear
convergence rates: the best-iterate and the last-iterate convergence rates. To
elaborate on these concepts, we assume that $D$ is a given metric (e.g.,
$\|Fx^{k}\|^{2}$ or $e(x^{k})^{2}$ defined by (2) below) defined on an iterate
sequence $\\{x^{k}\\}$ generated by the underlying algorithm for solving (NI)
or its special cases. For any $k\geq 0$ and a given order $\alpha>0$, if

$\min_{0\leq l\leq
k}D(x^{l})\leq\frac{1}{k+1}\sum_{l=0}^{k}D(x^{l})=\mathcal{O}\left(\frac{1}{k^{\alpha}}\right),$

then we say that $\left\\{x^{k}\right\\}$ has a
$\mathcal{O}\left(1/k^{\alpha}\right)$ best-iterate convergence rate. In this
case, we can take $\hat{x}_{k}:=x_{k_{\min}}$ with
$k_{\min}:=\mathrm{arg}\min_{0\leq l\leq k}D(x^{l})$ as the “best” output of
our algorithm. If we instead have
$D(x^{k})=\mathcal{O}\left(\frac{1}{k^{\alpha}}\right)$ with $x^{k}$ being the
$k$-th iterate, then we say that $\\{x^{k}\\}$ has a
$\mathcal{O}\left(1/k^{\alpha}\right)$ last-iterate convergence rate. We
emphasize that the convergence on the metric $D$ of $\\{x^{k}\\}$ does not
generally imply the convergence of $\\{x^{k}\\}$ itself, especially when
characterize the rate of convergence in different metrics.

###  2.3 Exact solutions and approximate solutions

There are different metrics to characterize exact and approximate solutions of
(NI). The most obvious one is the residual norm of $\Phi$, which is defined as

$e(x):=\min_{\xi\in Tx}\|Fx+\xi\|,\quad x\in\mathrm{dom}(\Phi).$ (2)

Clearly, if $e(x^{\star})=0$ for some $x^{\star}\in\mathrm{dom}(\Phi)$, then
$x^{\star}\in\mathrm{zer}(\Phi)$, a solution of (NI). If $T=0$, then
$e(x)=\|Fx\|$. However, if $e(\hat{x})\leq\epsilon$ for a given tolerance
$\epsilon>0$, then $\hat{x}$ can be considered as an $\epsilon$-approximate
solution of (NI). The algorithms presented in this paper use this metric as
one means to characterize approximate solutions.

Other metrics often used for monotone (VIP), a special case of (NI), are gap
functions and restricted gap functions [50, 72, 104], which are respectively
defined as

$\mathcal{G}(x):=\max_{y\in\mathcal{X}}\langle
Fy,y-x\rangle\quad\text{and}\quad\mathcal{G}_{\mathbb{B}}(x):=\max_{y\in\mathcal{X}\cap\mathbb{B}}\langle
Fy,y-x\rangle,$ (3)

where $\mathbb{B}$ is a given nonempty, closed, and bounded convex set. Note
that $\mathcal{G}(x)\geq 0$ for all $x\in\mathcal{X}$, and
$\mathcal{G}(x^{\star})=0$ iff $x^{\star}$ is a solution of (VIP). Therefore,
to characterize an $\epsilon$-approximate solution $\tilde{x}$ of (VIP), we
can impose $\mathcal{G}(\tilde{x})\leq\epsilon$. For the restricted gap
function $\mathcal{G}_{\mathbb{B}}$, if $x^{\star}$ is a solution of (VIP) and
$x^{\star}\in\mathbb{B}$, then $\mathcal{G}_{\mathbb{B}}(x^{\star})=0$.
Conversely, if $\mathcal{G}_{\mathbb{B}}(x^{\star})=0$ and
$x^{\star}\in\mathrm{int}\left(\mathbb{B}\right)$, then $x^{\star}$ is a
solution of (VIP) in $\mathbb{B}$ (see [104, Lemma 1]). For (DVIP), we can
also define similar dual gap functions and restricted dual gap functions
[104]. Gap functions have been widely used in the literature to characterize
approximate solutions generated by many numerical methods for solving (VIP) or
(DVIP), see, e.g., [33, 37, 50, 72, 100, 104].

If $J_{\eta T}$ is well-defined and single-valued for some $\eta>0$, and $F$
is single-valued, then we can use the following forward-backward splitting
residual operator:

$G_{\eta\Phi}x:=\tfrac{1}{\eta}\left(x-J_{\eta T}(x-\eta Fx)\right),$ (4)

to characterize solutions of (NI), where $F$ is single-valued and $J_{\eta T}$
is the resolvent of $\eta T$ for any $\eta>0$. It is clear that
$G_{\eta}x^{\star}=0$ iff $x^{\star}\in\mathrm{zer}(\Phi)$. In addition, if
$J_{\eta T}$ is firmly nonexpansive, then we also have

$\|G_{\eta\Phi}x\|\leq\|Fx+\xi\|,\quad(x,\xi)\in\mathrm{gra}(T).$ (5)

Hence, for a given tolerance $\epsilon>0$, if
$\|G_{\eta\Phi}\tilde{x}\|\leq\epsilon$, then we can say that $\tilde{x}$ is
an $\epsilon$-approximate solution of (NI). If $T:=\mathcal{N}_{\mathcal{X}}$,
i.e., (NI) reduces to (VIP), then, with $\eta=1$, $G_{\Phi}x$ reduces to the
classical natural map
$\Pi_{F,\mathcal{X}}x=x-\mathrm{proj}_{\mathcal{X}}(x-Fx)$ of (VIP), and
$r_{n}(x):=\|G_{\Phi}x\|=\|\Pi_{F,\mathcal{X}}x\|$ is the corresponding
natural residual at $x$. From (5), we have $r_{n}(x)\leq\|Fx+\xi\|$ for any
$\xi\in\mathcal{N}_{\mathcal{X}}(x)$.

###  2.4 Gradient/forward-type methods

Let us briefly recall the gradient/forward scheme for solving (NE) as follows.
Starting from $x^{0}\in\mathrm{dom}(F)$, at each iteration $k\geq 0$, we
update

$x^{k+1}:=x^{k}-\eta Fx^{k},$ (FW)

where $\eta>0$ is a given constant stepsize. If $F$ is $\rho$-co-coercive and
$0<\eta<\rho$, then $\\{x^{k}\\}$ converges to $x^{\star}\in\mathrm{zer}(F)$
(see, e.g., [50]). Otherwise, if $F$ is only monotone and $L$-Lipschitz
continuous, then there exist examples (e.g., $Fx=[x_{2},-x_{1}]$) showing that
(FW) is divergent for any choice of constant stepsize $\eta$.

To solve (NI), we can instead apply the forward-backward splitting method as
follows. Starting from $x^{0}\in\mathrm{dom}(F)$, at each iteration $k\geq 0$,
we update

$x^{k+1}:=J_{\eta T}(x^{k}-\eta Fx^{k}),$ (FBS)

where $\eta>0$ is a given constant stepsize. Similar to (FW), if $F$ is
$\rho$-co-coercive and $T$ is maximally monotone, then with $\eta\in(0,\rho)$,
$\\{x^{k}\\}$ generated by (FBS) converges to
$x^{\star}\in\mathrm{zer}(\Phi)$. If $F=\nabla{f}$, the gradient of a convex
and $L$-smooth function $f$, then $F$ is co-coercive. However, imposing the
co-coerciveness for a general mapping $F$ is often restrictive. Hence, both
(FW) and (FBS) are less practical.

##  3 Extragradient-Type Methods For Nonlinear Equations

As we have discussed before, the extragradient method was originally proposed
by G. M. Korpelevic in 1976 [73] and by A. S. Antipin around the same time [3]
to tackle saddle-point problems. Since then, this method has been extensively
studied in the literature, with numerous variants proposed (see, e.g., [50,
62, 71, 72, 91, 124, 126, 127, 142]). In recent years, the popularity of this
method has increased further due to its effectiveness in solving minimax
problems, including those in convex-concave and nonconvex-nonconcave settings,
which are common in machine learning and robust optimization. In this section,
we briefly survey both classical and recent works [43, 54, 56, 83] on the
extragradient method, as well as some closely related variants with minor
modifications. We unify the convergence analysis in one single theorem.

###  3.1 The class of extragradient methods for nonlinear equations

The class of extragradient methods for solving (NE) we discuss in this section
is presented as follows. Starting from an initial point
$x^{0}\in\mathrm{dom}(F)$, at each iteration $k\geq 0$, we update

$\left\\{\begin{array}[]{lcl}y^{k}&:=&x^{k}-\frac{\eta}{\beta}u^{k},\vspace{1ex}\\\
x^{k+1}&:=&x^{k}-\eta Fy^{k},\end{array}\right.$ (EG)

where $\eta>0$ is a given constant stepsize, $\beta\in(0,1]$ is a scaling
factor, and $u^{k}$ has two options as follows.

  * •

Option 1. If we set $u^{k}:=Fx^{k}$, then we obtain the extragradient scheme
[73] for (NE).

  * •

Option 2. If we choose $u^{k}:=Fy^{k-1}$, the we obtain the past-extragradient
method, also called Popov’s method [113], to solve (NE). This scheme is also
known as an optimistic gradient method in the literature, see also [38, 93,
96].

For Option 1 with $u^{k}:=Fx^{k}$, if $\beta=1$, then we obtain exactly the
classical extragradient method [73] for solving (NE). If $\beta<1$, then we
recover the extragradient-plus (EG$+$) scheme from [43] for solving (NE). If
we compute $x^{k}=y^{k}+\frac{\eta}{\beta}Fx^{k}$ from the first line of (EG)
and substitute it into the second line of (EG), then we get
$x^{k+1}=y^{k}-\eta(Fy^{k}-\frac{1}{\beta}Fx^{k})$. In this case, we obtain
from (EG) a forward-backward-forward splitting variant of Tseng’s method in
[136] as follows:

$\left\\{\begin{array}[]{lcl}y^{k}&:=&x^{k}-\frac{\eta}{\beta}Fx^{k},\vspace{1ex}\\\
x^{k+1}&:=&y^{k}-\eta(Fy^{k}-\frac{1}{\beta}Fx^{k}).\end{array}\right.$ (FBFS)

Clearly, if $\beta=1$, then we recover exactly Tseng’s method for solving
(NE).

For Option 2 with $u^{k}:=Fy^{k-1}$, we can show that it is equivalent to the
following variants. First, we can rewrite (EG) as

$\left\\{\begin{array}[]{lcl}x^{k+1}&:=&x^{k}-\eta Fy^{k}\vspace{1ex}\\\
y^{k+1}&:=&x^{k+1}-\tfrac{\eta}{\beta}Fy^{k}.\end{array}\right.$ (PEG)

This form shows us that (PEG) saves one evaluation $Fx^{k}$ of $F$ at each
iteration compared to Option 1. If $\beta=1$, then we obtain exactly the
Popov’s method in [113]. If we rotate the second line up and use $\beta=1$ as
$y^{k}=x^{k}-\eta Fy^{k-1}$, then we get the past-extragradient method.

Now, under this choice of $u^{k}$, from the first line of (EG), we have
$x^{k}=y^{k}+\frac{\eta}{\beta}u^{k}=y^{k}+\frac{\eta}{\beta}Fy^{k-1}$.
Substituting this expression into the first line of (PEG), we get
$x^{k+1}=y^{k}-\eta Fy^{k}+\frac{\eta}{\beta}Fy^{k-1}$. Substituting this
relation into the second line of (PEG), we can eliminate $x^{k+1}$ to get the
following variant:

$y^{k+1}:=y^{k}-\tfrac{\eta}{\beta}\big{(}(1+\beta)Fy^{k}-Fy^{k-1}\big{)}.$
(FRBS)

This scheme can be considered as a simplified variant of the forward-
reflected-backward splitting scheme in [89] for solving (NE) when we set
$\beta:=1$ as $y^{k+1}:=y^{k}-\eta(2Fy^{k}-Fy^{k-1})$.

Alternatively, from (PEG), we have $x^{k-1}-x^{k}=\eta Fy^{k-1}$ and
$\beta(x^{k}-y^{k})=\eta Fy^{k-1}$, leading to
$x^{k-1}-x^{k}=\beta(x^{k}-y^{k})$. Therefore, we get
$y^{k}=\frac{1}{\beta}((1+\beta)x^{k}-x^{k-1})$. Substituting this expression
into the first line of (PEG), we can show that

$x^{k+1}:=x^{k}-\eta F\big{(}\tfrac{1}{\beta}((1+\beta)x^{k}-x^{k-1})\big{)}.$
(RFB)

In particular, if $\beta=1$, then we obtain $x^{k+1}:=x^{k}-\eta
F(2x^{k}-x^{k-1})$, which turns out to be the reflected gradient method in
[87] or the reflected forward-backward splitting scheme in [30] for solving
(NE).

Using the relation $x^{k-1}-x^{k}=\beta(x^{k}-y^{k})$ above, we can compute
that
$x^{k}=\frac{\beta}{1+\beta}y^{k}+\frac{1}{1+\beta}x^{k-1}=\frac{(\omega-1)}{\omega}y^{k}+\frac{1}{\omega}x^{k-1}$,
where $\omega:=1+\beta$. Combining the two lines of (EG), we get
$y^{k+1}:=x^{k+1}-\tfrac{\eta}{\beta}Fy^{k}=x^{k}-\frac{\eta(1+\beta)}{\beta}Fy^{k}$.
Putting both expressions together, we get

$\left\\{\begin{array}[]{lcl}x^{k}&:=&\tfrac{(\omega-1)}{\omega}y^{k}+\tfrac{1}{\omega}x^{k-1},\vspace{1ex}\\\
y^{k+1}&:=&x^{k}-\tfrac{\eta(1+\beta)}{\beta}Fy^{k}.\end{array}\right.$ (GR)

This method is a simplified variant of the golden-ratio method in [88] for
solving (NE). Overall, the template (EG) covers a class of EG algorithms with
many common instances as discussed.

###  3.2 Convergence analysis

The results presented in this section were obtained in [83], but here we
provide a different proof and unify several methods in one. To analyze the
convergence of (EG), we first prove the following lemmas.

######  Lemma 3.1.

If $\left\\{(x^{k},y^{k})\right\\}$ is generated by (EG), then for any
$\gamma>0$ and any $\hat{x}\in\mathrm{dom}(F)$, we have

$\begin{array}[]{lcl}\|x^{k+1}-\hat{x}\|^{2}&\leq&\|x^{k}-\hat{x}\|^{2}-\beta\|y^{k}-x^{k}\|^{2}+\tfrac{\eta^{2}}{\gamma}\|Fy^{k}-u^{k}\|^{2}-2\eta\langle
Fy^{k},y^{k}-\hat{x}\rangle\vspace{1ex}\\\
&&-{~{}}(\beta-\gamma)\|x^{k+1}-y^{k}\|^{2}-(1-\beta)\|x^{k+1}-x^{k}\|^{2}.\end{array}$
(6)

###### Proof.

First, for any $\hat{x}\in\mathrm{dom}(F)$, using $x^{k+1}-x^{k}=-\eta Fy^{k}$
from the second line of (EG), we have

$\begin{array}[]{lcl}\|x^{k+1}-\hat{x}\|^{2}&=&\|x^{k}-\hat{x}\|^{2}+2\langle
x^{k+1}-x^{k},x^{k+1}-\hat{x}\rangle-\|x^{k+1}-x^{k}\|^{2}\vspace{1ex}\\\
&=&\|x^{k}-\hat{x}\|^{2}-2\eta\langle
Fy^{k},x^{k+1}-\hat{x}\rangle-\|x^{k+1}-x^{k}\|^{2}.\end{array}$

Next, using $\eta u^{k}=\beta(x^{k}-y^{k})$ from the first line of (EG), the
Cauchy-Schwarz inequality, the identity $2\langle
x^{k+1}-y^{k},x^{k}-y^{k}\rangle=\|x^{k}-y^{k}\|^{2}+\|x^{k+1}-y^{k}\|^{2}-\|x^{k+1}-x^{k}\|^{2}$,
and an elementary inequality $2wz\leq\gamma w^{2}+\frac{z^{2}}{\gamma}$ for
any $\gamma>0$ and $w,z\geq 0$, we can derive that

$\begin{array}[]{lcl}2\eta\langle Fy^{k},x^{k+1}-\hat{x}\rangle&=&2\eta\langle
Fy^{k},y^{k}-\hat{x}\rangle+2\eta\langle
Fy^{k}-u^{k},x^{k+1}-y^{k}\rangle+2\eta\langle
u^{k},x^{k+1}-y^{k}\rangle\vspace{1ex}\\\ &\geq&2\eta\langle
Fy^{k},y^{k}-\hat{x}\rangle-2\eta\|Fy^{k}-u^{k}\|\|x^{k+1}-y^{k}\|+2\beta\langle
x^{k+1}-y^{k},x^{k}-y^{k}\rangle\vspace{1ex}\\\ &\geq&2\eta\langle
Fy^{k},y^{k}-\hat{x}\rangle-\frac{\eta^{2}}{\gamma}\|Fy^{k}-u^{k}\|^{2}-\gamma\|x^{k+1}-y^{k}\|^{2}\vspace{1ex}\\\
&&+{~{}}\beta\big{[}\|x^{k}-y^{k}\|^{2}+\|x^{k+1}-y^{k}\|^{2}-\|x^{k+1}-x^{k}\|^{2}\big{]}\vspace{1ex}\\\
&=&2\eta\langle
Fy^{k},y^{k}-\hat{x}\rangle+\beta\|y^{k}-x^{k}\|^{2}-\frac{\eta^{2}}{\gamma}\|Fy^{k}-u^{k}\|^{2}\vspace{1ex}\\\
&&+{~{}}(\beta-\gamma)\|x^{k+1}-y^{k}\|^{2}-\beta\|x^{k+1}-x^{k}\|^{2}.\end{array}$

Finally, combining the last two expressions, we obtain (6). ∎

######  Lemma 3.2.

Let $F$ be $\rho$-co-hypomonotone, i.e. there exists $\rho\geq 0$ such that
$\langle Fx-Fy,x-y\rangle\geq-\rho\|Fx-Fy\|^{2}$ for all
$x,y\in\mathrm{dom}(F)$ and $L$-Lipschitz continuous. Let
$\left\\{(x^{k},y^{k})\right\\}$ be generated by (EG). Then, for any $c>0$ and
$\omega>0$, we have

$\begin{array}[]{lcl}\|Fx^{k+1}\|^{2}&\leq&\|Fx^{k}\|^{2}-\frac{[c\eta-2(1+c)\rho]}{c\eta}\|Fy^{k}-Fx^{k}\|^{2}+\frac{\left[\eta\omega+2(1+c)\rho\right]L^{2}\eta}{\beta^{2}}\|\beta
Fy^{k}-u^{k}\|^{2}\vspace{1ex}\\\
&&-{~{}}(\omega-1)\|Fx^{k+1}-Fy^{k}\|^{2}.\end{array}$ (7)

###### Proof.

Since $F$ is $\rho$-cohypomonotone, we have $\langle
Fx^{k+1}-Fx^{k},x^{k+1}-x^{k}\rangle+\rho\|Fx^{k+1}-Fx^{k}\|^{2}\geq 0$.
Substituting $x^{k+1}-x^{k}=-\eta Fy^{k}$ from the second line of (EG) into
this inequality, we can show that

$\begin{array}[]{lcl}0&\leq&2\langle Fx^{k},Fy^{k}\rangle-2\langle
Fx^{k+1},Fy^{k}\rangle+\frac{2\rho}{\eta}\|Fx^{k+1}-Fx^{k}\|^{2}\vspace{1ex}\\\
&\leq&\|Fx^{k}\|^{2}-\|Fy^{k}-Fx^{k}\|^{2}-\|Fx^{k+1}\|^{2}+\|Fx^{k+1}-Fy^{k}\|^{2}+\frac{2\rho}{\eta}\|Fx^{k+1}-Fx^{k}\|^{2}.\end{array}$

Now, by utilizing Young’s inequality, the $L$-Lipschitz continuity of $F$, and
$x^{k+1}-y^{k}=-\eta(Fy^{k}-\frac{1}{\beta}u^{k})$ from (EG), for any $c>0$
and $\omega\geq 1$, the last estimate leads to

$\begin{array}[]{lcl}\|Fx^{k+1}\|^{2}&\leq&\|Fx^{k}\|^{2}-\|Fy^{k}-Fx^{k}\|^{2}+\omega\|Fx^{k+1}-Fy^{k}\|^{2}+\frac{2\rho}{\eta}\|Fx^{k+1}-Fx^{k}\|^{2}-(\omega-1)\|Fx^{k+1}-Fy^{k}\|^{2}\vspace{1ex}\\\
&\leq&\|Fx^{k}\|^{2}-\frac{c\eta-2(1+c)\rho}{c\eta}\|Fy^{k}-Fx^{k}\|^{2}+\frac{[\eta\omega+2(1+c)\rho]L^{2}}{\eta}\|x^{k+1}-y^{k}\|^{2}-(\omega-1)\|Fx^{k+1}-Fy^{k}\|^{2}\vspace{1ex}\\\
&\leq&\|Fx^{k}\|^{2}-\frac{c\eta-2(1+c)\rho}{c\eta}\|Fy^{k}-Fx^{k}\|^{2}+\left[\eta\omega+2(1+c)\rho\right]L^{2}\eta\|Fy^{k}-\frac{1}{\beta}u^{k}\|^{2}\vspace{1ex}\\\
&&-{~{}}(\omega-1)\|Fx^{k+1}-Fy^{k}\|^{2},\end{array}$

which exactly proves (7). ∎

Now, we are ready to establish both the best-iterate and the last-iterate
convergence rates of (EG).

######  Theorem 3.1.

Suppose that $F$ in (NE) is $L$-Lipschitz continuous and
$\mathrm{zer}(F)\neq\emptyset$. Let $\left\\{(x^{k},y^{k})\right\\}$ be
generated by (EG) for solving (NE). Then, we have the following statements.

  * $\mathrm{(a)}$

$($Extragradient method$)$ Let us choose $u^{k}:=Fx^{k}$ and assume that there
exists $\rho\geq 0$ such that $\langle
Fx,x-x^{\star}\rangle\geq-\rho\|Fx\|^{2}$ for all $x\in\mathrm{dom}(F)$ and a
given $x^{\star}\in\mathrm{zer}(F)$ $($this condition holds if, in particular,
$F$ is $\rho$-co-hypomonotone on $\mathrm{dom}(F)$$)$. Then, if
$L\rho\leq\frac{3\sqrt{2}-2}{12}\approx 0.1869$, $\beta\in(0,1]$, and $\eta$
is chosen such that

$\begin{array}[]{l}0\leq\frac{\beta[1-\sqrt{1-24L\rho(3L\rho+1)}]}{2L(3L\rho+1)}<\eta<\frac{\beta[1+\sqrt{1-24L\rho(3L\rho+1)}]}{2L(3L\rho+1)}\leq\frac{\beta}{L},\end{array}$
(8)

then we have

$\min_{0\leq l\leq
k}\|Fx^{l}\|^{2}\leq\frac{1}{k+1}\sum_{l=0}^{k}\|Fx^{l}\|^{2}\leq\frac{C_{\rho}\|x^{0}-x^{\star}\|^{2}}{k+1},$
(9)

where
$C_{\rho}:=\frac{\beta^{2}}{\eta[\eta\beta-6\beta^{2}\rho-(3L\rho+1)L\eta^{2}]}>0$.
Consequently, $\\{\|x^{k}-x^{\star}\|\\}$ is nonincreasing and
$\lim_{k\to\infty}\|x^{k}-y^{k}\|=\lim_{k\to\infty}\|Fx^{k}\|=\lim_{k\to\infty}\|Fy^{k}\|=0$.
Moreover, we have $\min_{0\leq l\leq
k}\|Fx^{l}\|=\mathcal{O}\big{(}1/\sqrt{k}\big{)}$ showing the
$\mathcal{O}\big{(}1/\sqrt{k}\big{)}$ best-iterate convergence rate of
$\\{x^{k}\\}$.

In particular, if $\beta:=1$ and $F$ is $\rho$-co-hypomonotone on
$\mathrm{dom}(F)$ such that $L\rho\leq\frac{3\sqrt{2}-2}{12}$, then

$\|Fx^{k+1}\|^{2}\leq\|Fx^{k}\|^{2}-\psi\cdot\|Fy^{k}-Fx^{k}\|^{2}\quad\text{and}\quad\|Fx^{k}\|\leq\frac{\sqrt{C_{\rho}}\|x^{0}-x^{\star}\|}{\sqrt{k+1}},$
(10)

where $\psi:=1-\frac{4\rho}{\eta}-L^{2}\eta(\eta+4\rho)>0$. Hence, we have
$\|Fx^{k}\|=\mathcal{O}\big{(}1/\sqrt{k}\big{)}$ on the last-iterate $x^{k}$.

  * $\mathrm{(b)}$

$($Past-extragradient method$)$ Let us choose $u^{k}:=Fy^{k-1}$ and
$y^{-1}:=x^{0}$ and assume that there exists $\rho\geq 0$ such that $\langle
Fx,x-x^{\star}\rangle\geq-\rho\|Fx\|^{2}$ for all $x\in\mathrm{dom}(F)$ and a
given $x^{\star}\in\mathrm{zer}(F)$ $($in particular, if $F$ is $\rho$-co-
hypomonotone on $\mathrm{dom}(F)$$)$. Then, for fixed $\beta\in(0,1]$, if
$L\rho\leq\frac{\beta^{2}}{12}$ and $\eta$ is chosen such that

$\begin{array}[]{l}0\leq\frac{\beta-\sqrt{\beta^{2}-12L\rho}}{6L}<\eta<\frac{\beta+\sqrt{\beta^{2}-12L\rho}}{6L}\leq\frac{\beta}{3L},\end{array}$
(11)

then we have

$\min_{0\leq l\leq
k}\|Fx^{l}\|^{2}\leq\frac{1}{k+1}\sum_{l=0}^{k}\big{[}\|Fx^{l}\|^{2}+\kappa\|x^{l}-y^{l-1}\|^{2}\big{]}\leq\frac{\hat{C}_{\rho}\|x^{0}-x^{\star}\|^{2}}{k+1},$
(12)

where $\kappa:=\frac{(2L^{2}\eta^{2}+3)L}{2(\beta\eta-3L\eta^{2}-4\rho)}>0$
and
$\hat{C}_{\rho}:=\frac{2L^{2}\eta^{2}+3}{(\beta\eta-3L\eta^{2}-4\rho)\eta}>0$.
Consequently, we also have
$\lim_{k\to\infty}\|x^{k}-y^{k}\|=\lim_{k\to\infty}\|x^{k+1}-y^{k}\|=\lim_{k\to\infty}\|Fx^{k}\|=\lim_{k\to\infty}\|Fy^{k}\|=0$,
and $\min_{0\leq l\leq k}\|Fx^{l}\|=\mathcal{O}\big{(}1/\sqrt{k}\big{)}$
showing the $\mathcal{O}\big{(}1/\sqrt{k}\big{)}$ best-iterate convergence
rate of $\\{x^{k}\\}$.

In particular, if $\beta:=1$ and $F$ is $\rho$-co-hypomonotone on
$\mathrm{dom}(F)$ such that $12L\rho\leq 1$, then

$\|Fx^{k+1}\|^{2}+\hat{\kappa}\|Fx^{k+1}-Fy^{k}\|^{2}\leq\|Fx^{k}\|^{2}+\hat{\kappa}\|Fx^{k}-Fy^{k-1}\|^{2}\quad\text{and}\quad\|Fx^{k}\|\leq\frac{\sqrt{M_{\rho}}\|x^{0}-x^{\star}\|}{\sqrt{k+1}},$
(13)

where $\hat{\kappa}:=\frac{2(\eta+4\rho)L^{2}\eta}{1-2L^{2}\eta^{2}}$ and
$M_{\rho}:=\hat{C}_{\rho}\cdot\max\left\\{\frac{L^{2}\hat{\kappa}}{\kappa},1\right\\}$.
Moreover, we have the last-iterate convergence rate as
$\|Fx^{k}\|=\mathcal{O}\big{(}1/\sqrt{k}\big{)}$.

Before proving Theorem 3.1, we give the following remarks.

######  Remark 3.1.

If $\rho=0$ in Theorem 3.1, i.e. $F$ is star-monotone (and in particular,
monotone), then our condition on the stepsize $\eta$ reduces to
$0<\eta<\frac{\beta}{L}$ for the extragradient method and
$0<\eta<\frac{\beta}{3L}$ for the past-extragradient method. These choices are
standard and often seen in both methods. However, we have not yet optimized
the choice of these parameters for (EG) stated in Theorem 3.1. Here, the
stepsize $\eta$ can be carefully chosen so that we can possibly enlarge the
range of $L\rho$ (see [57] as an example).

######  Remark 3.2.

The results in Theorem 3.1 were proven in [83], and then were revised in [57].
The last-iterate convergence rates were proven in previous works such as [54]
for the monotone case but with an additional assumption. Note that the best-
iterate rates for the monotone or the star-monotone case are classical, which
can be found, e.g., in [50, 73]. The last-iterate convergence for the monotone
case can be found in recent works such as [54, 56]. The best-iterate rates for
the co-hypomonotone or the star-co-hypomonotone case can be found in [43],
while the last-iterate convergence rates were recently proven in [83].
Nevertheless, in this survey, we provide a unified analysis for all of these
variants of EG, which covers both the monotone and co-hypomonotone cases
altogether. Our analysis is also different from [83].

######  Remark 3.3.

One can easily modify the proof of Theorem 3.1 to handle the star-strongly-
monotone case of $F$. Indeed, if $F$ is $\mu$-star-strongly monotone, then we
have $\mu\|x^{k}-x^{\star}\|^{2}\leq\langle
Fx^{k},x^{k}-x^{\star}\rangle\leq\|Fx^{k}\|\|x^{k}-x^{\star}\|$. This leads to
$\|Fx^{k}\|^{2}\geq\mu^{2}\|x^{k}-x^{\star}\|^{2}$. Using this inequality,
$\hat{x}:=x^{\star}$, and $u^{k}:=Fx^{k}$ into (6), we obtain
$\|x^{k+1}-x^{\star}\|^{2}\leq(1-\eta^{2}(1-L^{2}\eta^{2})\mu^{2})\|x^{k}-x^{\star}\|^{2}$.
Clearly, if we choose $\eta\in\big{(}0,\frac{1}{L}\big{)}$, then
$\varphi:=1-\eta^{2}(1-L^{2}\eta^{2})\mu^{2}\in(0,1)$, and we obtain a linear
convergence rate of $\\{\|x^{k}-x^{\star}\|^{2}\\}$ with a contraction factor
$\varphi$. Note that since $\mu\leq L$, we have $\mu^{4}-4L^{2}\mu^{2}<0$.
Hence, $L^{2}\mu^{2}\eta^{4}-\mu^{2}\eta^{2}+1\geq 0$ always holds to
guarantee that $\varphi\in(0,1)$. Another proof for the monotone case can be
found, e.g., in [50].

###### Proof of Theorem 3.1.

(a) Extragradient method. Since $u^{k}:=Fx^{k}$, by the $L$-Lipschitz
continuity of $F$, we have $\|Fy^{k}-u^{k}\|=\|Fy^{k}-Fx^{k}\|\leq
L\|x^{k}-y^{k}\|$. Using this inequality,
$\hat{x}:=x^{\star}\in\mathrm{zer}(F)$, and $\langle
Fy^{k},y^{k}-x^{\star}\rangle\geq-\rho\|Fy^{k}\|^{2}$ into (6), we have

$\begin{array}[]{lcl}\|x^{k+1}-x^{\star}\|^{2}&\leq&\|x^{k}-x^{\star}\|^{2}-\left(\beta-\tfrac{L^{2}\eta^{2}}{\gamma}\right)\|y^{k}-x^{k}\|^{2}+2\eta\rho\|Fy^{k}\|^{2}\vspace{1ex}\\\
&&-{~{}}(\beta-\gamma)\|x^{k+1}-y^{k}\|^{2}-(1-\beta)\|x^{k+1}-x^{k}\|^{2}.\end{array}$
(14)

Now, using the first line of (EG) with $u^{k}:=Fx^{k}$ as
$\beta(x^{k}-y^{k})=\eta Fx^{k}$, the $L$-Lipschitz continuity of $F$, and
Young’s inequality, we have

$\begin{array}[]{lcl}\|Fy^{k}\|^{2}&\leq&\frac{3}{2}\|Fy^{k}-Fx^{k}\|^{2}+3\|Fx^{k}\|^{2}\leq\frac{3(2\beta^{2}+L^{2}\eta^{2})}{2\eta^{2}}\|x^{k}-y^{k}\|^{2}.\end{array}$
(15)

Substituting this inequality, $\gamma:=L\eta$, and
$\beta^{2}\|y^{k}-x^{k}\|^{2}=\eta^{2}\|Fx^{k}\|^{2}$ into (14), we obtain

$\begin{array}[]{lcl}\|x^{k+1}-x^{\star}\|^{2}&\leq&\|x^{k}-x^{\star}\|^{2}-(\beta-L\eta)\|x^{k+1}-y^{k}\|^{2}-(1-\beta)\|x^{k+1}-x^{k}\|^{2}\vspace{1ex}\\\
&&-{~{}}\frac{\eta[\eta\beta-6\beta^{2}\rho-(3L\rho+1)L\eta^{2}]}{\beta^{2}}\|Fx^{k}\|^{2}.\end{array}$
(16)

Let us choose $\eta>0$ such that
$\eta\beta-6\beta^{2}\rho-(3L\rho+1)L\eta^{2}>0$, which holds if $\eta$
satisfies

$\begin{array}[]{lcl}0\leq\frac{\beta\big{[}1-\sqrt{1-24L\rho(3L\rho+1)}\big{]}}{2L(3L\rho+1)}<\eta<\frac{\beta\big{[}1+\sqrt{1-24L\rho(3L\rho+1)}\big{]}}{2L(3L\rho+1)}\leq\frac{\beta}{L},\end{array}$

provided that $L\rho\leq\frac{3\sqrt{2}-2}{12}\approx 0.1869$. This condition
is exactly (8). In this case, we also have $\beta-L\eta\geq 0$, and (16)
implies (9). The next statement is a consequence of (16) combining with (15)
using standard arguments.

Next, if $F$ is $\rho$-co-hypomonotone, then with $\beta:=1$, using $\|\beta
Fy^{k}-u^{k}\|=\|Fy^{k}-Fx^{k}\|$ and $c=\omega=1$ into (7), we obtain

$\begin{array}[]{lcl}\|Fx^{k+1}\|^{2}\leq\|Fx^{k}\|^{2}-\frac{\left[\eta-4\rho-L^{2}\eta^{2}(\eta+4\rho)\right]}{\eta}\|Fy^{k}-Fx^{k}\|^{2}.\end{array}$

However, by the choice of $\eta$ as in (8), one has
$\eta-6\rho-(3L\rho+1)L\eta^{2}>0$. It is obvious to check that
$\psi:=\eta-4\rho-L^{2}\eta^{2}(\eta+4\rho)\geq\eta-6\rho-(3L\rho+1)L\eta^{2}>0$.
This condition leads to the first part of (10). The second part of (10) is a
consequence of the first part and (9).

(b) Past-extragradient method. If we choose $u^{k}:=Fy^{k-1}$, then by Young’s
inequality and the $L$-Lipschitz continuity of $F$, we have

$\begin{array}[]{lcl}\|Fy^{k}-u^{k}\|^{2}&=&\|Fy^{k}-Fy^{k-1}\|^{2}\leq
3\|Fy^{k}-Fx^{k}\|^{2}+\frac{3}{2}\|Fx^{k}-Fy^{k-1}\|^{2}\vspace{1ex}\\\
&\leq&3L^{2}\|x^{k}-y^{k}\|^{2}+\frac{3L^{2}}{2}\|x^{k}-y^{k-1}\|^{2}.\end{array}$
(17)

Moreover, since $\langle Fx,x-x^{\star}\rangle\geq-\rho\|Fx\|^{2}$ for all
$x\in\mathrm{dom}(F)$, using this condition, $\eta Fy^{k}=x^{k+1}-x^{k}$ from
(EG), and Young’s inequality, we can lower bound that

$\begin{array}[]{lcl}2\eta\langle
Fy^{k},y^{k}-x^{\star}\rangle&\geq&-2\rho\eta\|Fy^{k}\|^{2}=-\frac{2\rho}{\eta}\|x^{k+1}-x^{k}\|^{2}\geq-\frac{4\rho}{\eta}\big{[}\|x^{k+1}-y^{k}\|^{2}+\|y^{k}-x^{k}\|^{2}\big{]}.\end{array}$
(18)

Substituting (17), (18), and $\hat{x}:=x^{\star}\in\mathrm{zer}(F)$ into (6),
we obtain

$\begin{array}[]{lcl}\|x^{k+1}-x^{\star}\|^{2}&+&\left(\frac{3L^{2}\eta^{2}}{\gamma}-\gamma\right)\|x^{k+1}-y^{k}\|^{2}\leq\|x^{k}-x^{\star}\|^{2}+\left(\tfrac{3L^{2}\eta^{2}}{\gamma}-\gamma\right)\|x^{k}-y^{k-1}\|^{2}\vspace{1ex}\\\
&&-{~{}}\left(\beta-\frac{3L^{2}\eta^{2}}{\gamma}-\frac{4\rho}{\eta}\right)\left[\|x^{k+1}-y^{k}\|^{2}+\|y^{k}-x^{k}\|^{2}\right]-\left(\tfrac{3L^{2}\eta^{2}}{2\gamma}-\gamma\right)\|x^{k}-y^{k-1}\|^{2}\vspace{1ex}\\\
&&-{~{}}(1-\beta)\|x^{k+1}-x^{k}\|^{2}.\end{array}$ (19)

Next, by the second line of (EG) and Young’s inequality, we can easily show
that
$\|Fy^{k}\|^{2}=\frac{1}{\eta^{2}}\|x^{k+1}-x^{k}\|^{2}\leq\frac{3}{2\eta^{2}}\|x^{k}-y^{k}\|^{2}+\frac{3}{\eta^{2}}\|x^{k+1}-y^{k}\|^{2}$.
Alternatively, by Young’s inequality and the $L$-Lipschitz continuity of $F$,
we also have
$\|Fx^{k}\|^{2}\leq\left(L^{2}+\frac{3}{2\eta^{2}}\right)\|x^{k}-y^{k}\|^{2}+\left(1+\frac{2L^{2}\eta^{2}}{3}\right)\|Fy^{k}\|^{2}$.
Combining both inequalities, we get

$\begin{array}[]{lcl}\|Fx^{k}\|^{2}\leq\frac{2L^{2}\eta^{2}+3}{\eta^{2}}\left[\|x^{k}-y^{k}\|^{2}+\|x^{k+1}-y^{k}\|^{2}\right].\end{array}$

Now, let us choose $\gamma:=L\eta$ and using the last inequality into (19), we
can show that

$\begin{array}[]{lcl}\|x^{k+1}-x^{\star}\|^{2}+2L\eta\|x^{k+1}-y^{k}\|^{2}&\leq&\|x^{k}-x^{\star}\|^{2}+2L\eta\|x^{k}-y^{k-1}\|^{2}-(1-\beta)\|x^{k+1}-x^{k}\|^{2}\vspace{1ex}\\\
&&-{~{}}\frac{\eta(\beta\eta-3L\eta^{2}-4\rho)}{2L^{2}\eta^{2}+3}\|Fx^{k}\|^{2}-\tfrac{L\eta}{2}\|x^{k}-y^{k-1}\|^{2}.\end{array}$
(20)

If we choose $\eta>0$ such that $\beta\eta-3L\eta^{2}-4\rho>0$, then (20)
implies (12). The next statement of Part (b) follows from (12) and the fact
that
$\|Fy^{k-1}\|^{2}\leq\frac{L^{2}+\kappa}{\kappa}\big{[}\|Fx^{k}\|^{2}+\kappa\|x^{k}-y^{k-1}\|^{2}\big{]}$.
Note that, the condition $\beta\eta-3L\eta^{2}-4\rho>0$ holds if
$0\leq\frac{\beta-\sqrt{\beta^{2}-12L\rho}}{6L}<\eta<\frac{\beta+\sqrt{\beta^{2}-12L\rho}}{6L}\leq\frac{\beta}{3L}$,
provided that $12L\rho\leq\beta^{2}$ as stated in Theorem 3.1.

To prove the last-iterate rate, let us choose
$\omega:=\frac{1+8L^{2}\rho\eta}{1-2L^{2}\eta^{2}}>1$, provided that
$\sqrt{2}L\eta<1$. Moreover, for $\beta=1$ and $u^{k}:=Fy^{k-1}$, we have
$\|\beta Fy^{k}-u^{k}\|^{2}=\|Fy^{k}-Fy^{k-1}\|^{2}\leq
2\|Fy^{k}-Fx^{k}\|^{2}+2\|Fx^{k}-Fy^{k-1}\|^{2}$. Substituting this
inequality, $c:=1$, and $\omega:=\frac{1+8L^{2}\rho\eta}{1-2L^{2}\eta^{2}}$
into (7), we obtain

$\begin{array}[]{lcl}\|Fx^{k+1}\|^{2}+(\omega-1)\|Fx^{k+1}-Fy^{k}\|^{2}&\leq&\|Fx^{k}\|^{2}+(\omega-1)\|Fx^{k}-Fy^{k-1}\|^{2}\vspace{1ex}\\\
&&-{~{}}\left[1-\frac{4\rho}{\eta}-\frac{2L^{2}\eta(\eta+4\rho)}{1-2L^{2}\eta^{2}}\right]\|Fy^{k}-Fx^{k}\|^{2}.\end{array}$
(21)

It is obvious to show that the conditions $1-\frac{4\rho}{\eta}-3L\eta>0$ and
$\sqrt{2}L\eta<1$ guarantee that
$1-\frac{4\rho}{\eta}-\frac{2L^{2}\eta(\eta+4\rho)}{1-2L^{2}\eta^{2}}\geq 0$.
Hence, if we define
$\hat{\kappa}:=\omega-1=\frac{2(\eta+4\rho)L^{2}\eta}{1-2L^{2}\eta^{2}}$, then
(21) reduces to

$\begin{array}[]{lcl}\|Fx^{k+1}\|^{2}+\hat{\kappa}\|Fx^{k+1}-Fy^{k}\|^{2}&\leq&\|Fx^{k}\|^{2}+\hat{\kappa}\|Fx^{k}-Fy^{k-1}\|^{2}.\end{array}$
(22)

For $C_{0}:=\max\left\\{\frac{L^{2}\hat{\kappa}}{\kappa},1\right\\}$, we have
$\|Fx^{k}\|^{2}+\hat{\kappa}\|Fx^{k}-Fy^{k-1}\|^{2}\leq
C_{0}\left[\|Fx^{k}\|^{2}+\frac{\kappa}{L^{2}}\|Fx^{k}-Fy^{k-1}\|^{2}\right]\leq
C_{0}\left[\|Fx^{k}\|^{2}+\kappa\|x^{k}-y^{k-1}\|^{2}\right]$. Combining this
inequality and (12), we get

$\frac{1}{k+1}\sum_{l=0}^{k}\big{[}\|Fx^{l}\|^{2}+\hat{\kappa}\|Fx^{l}-Fy^{l-1}\|^{2}\big{]}\leq\frac{C_{0}}{k+1}\sum_{l=0}^{k}\big{[}\|Fx^{l}\|^{2}+\kappa\|x^{l}-y^{l-1}\|^{2}\big{]}\leq\frac{C_{0}\hat{C}_{\rho}\|x^{0}-x^{\star}\|^{2}}{k+1},$

Using (22) into the last bound, we obtain (13). Note that the condition
$3L\eta<1$ guarantees that $\sqrt{2}L\eta<1$. The remaining statement of (b)
in Theorem 3.1 is a direct consequence of (13) and (18). ∎

##  4 Extragradient-Type Methods for Monotone Inclusions

In this section, we go beyond (NE) to survey recent results on both best-
iterate and last-iterate convergence rates of the EG method and its variants
for solving (NI). Again, we provide a unified analysis that covers a wide
class of EG variants of the monotone instances of (NI) as can be seen below.

###  4.1 The class of extragradient methods

The class of EG methods for solving (NI) we consider in this section can be
described as follows. Starting from an initial point
$x^{0}\in\mathrm{dom}(\Phi)$, at each iteration $k\geq 0$, we update

$\left\\{\begin{array}[]{lcl}y^{k}&:=&J_{\frac{\eta}{\beta}T}(x^{k}-\frac{\eta}{\beta}u^{k}),\vspace{1ex}\\\
x^{k+1}&:=&J_{\eta T}(x^{k}-\eta Fy^{k}),\end{array}\right.$ (EG2)

where $J_{\eta T}$ is the resolvent of $\eta T$, $\eta>0$ is a given stepsize,
and $\beta>0$ is a scaling factor. Here, we consider two different choices of
$u^{k}$ as follows:

  * •

Option 1. If $u^{k}:=Fx^{k}$, then we obtain
$y^{k}:=J_{\frac{\eta}{\beta}T}(x^{k}-\frac{\eta}{\beta}Fx^{k})$, leading to
the well-known extragradient method (or extragradient-plus – EG+ in [43] if
$\beta\in(0,1)$) for solving (NI).

  * •

Option 2. If $u^{k}:=Fy^{k-1}$, then we obtain
$y^{k}:=J_{\frac{\eta}{\beta}T}(x^{k}-\frac{\eta}{\beta}Fy^{k-1})$, leading to
the past-extragradient method (or equivalently, Popov’s method [113]) for
solving (NI).

Clearly, when $T=\mathcal{N}_{\mathcal{X}}$, the normal cone of a nonempty,
closed, and convex set $\mathcal{X}$, then $J_{\gamma
T}=\mathrm{proj}_{\mathcal{X}}$, the projection onto $\mathcal{X}$ and hence,
(EG2) reduces to the extragradient variant for solving (VIP) widely studied in
the literature [50, 72]. In terms of computational complexity, (EG2) requires
two evaluations of $F$ at $x^{k}$ and $y^{k}$, and two evaluations of the
resolvent $J_{\eta T}$ at each iteration. It costs as twice as one iteration
of the forward-backward splitting method (FBS). However, its does not require
the co-coerciveness of $F$ to guarantee convergence. Again, we use a scaling
factor $\beta$ as in (EG), which covers EG+ in [43] as a special case.

Now, for given $\zeta^{k}\in Ty^{k}$ and $\xi^{k+1}\in Tx^{k+1}$, we denote
$\tilde{w}^{k}:=Fx^{k}+\zeta^{k}$, and $\hat{w}^{k+1}:=Fy^{k}+\xi^{k+1}$.
Then, we can rewrite (EG2) equivalently to

$\left\\{\begin{array}[]{lclclcll}y^{k}&:=&x^{k}-\frac{\eta}{\beta}(u^{k}+\zeta^{k})&=&x^{k}-\frac{\eta}{\beta}(\tilde{w}^{k}+u^{k}-Fx^{k}),&\zeta^{k}&\in&Ty^{k},\vspace{1ex}\\\
x^{k+1}&:=&x^{k}-\eta(Fy^{k}+\xi^{k+1})&=&x^{k}-\eta\hat{w}^{k+1},&\xi^{k+1}&\in&Tx^{k+1}.\end{array}\right.$
(23)

This representation makes (EG2) looks like (EG), and it is a key step for our
convergence analysis.

###  4.2 One-iteration analysis

We establish both the best-iterate and last-iterate convergence rates of (EG2)
under the assumption that $F$ is monotone and $T$ is maximally $3$-cyclically
monotone. Note that if $T$ is maximally cyclically monotone, then
$T=\partial{g}$, the subdifferential of a proper, closed, and convex function
due to [11, Theorem 22.18]. In this case, (NI) reduces to (MVIP). However, we
do not require $T$ to be maximally cyclically monotone, but only $3$-maximally
cyclically monotone, which may not be necessarily identical to $\partial{g}$.
Therefore, our result below is more general than existing variants in the
recent literature, including [26].

To analyze the convergence of (EG2), we also define

$w^{k}:=Fx^{k}+\xi^{k}\quad\text{for some}\quad\xi^{k}\in Tx^{k}.$ (24)

The following lemma provides key estimates to establish convergence of (EG2).

######  Lemma 4.1.

Suppose that $\\{(x^{k},y^{k})\\}$ is generated by (EG2), $w^{k}$ is defined
by (24) and $T$ is maximally $3$-cyclically monotone. Then, for any
$\gamma>0$, any $x^{\star}\in\mathrm{zer}(\Phi)$, we have

$\begin{array}[]{lcl}\|x^{k+1}-x^{\star}\|^{2}&\leq&\|x^{k}-x^{\star}\|^{2}-(1-\beta)\|x^{k+1}-x^{k}\|^{2}-(\beta-\gamma)\|x^{k+1}-y^{k}\|^{2}\vspace{1ex}\\\
&&-{~{}}\beta\|x^{k}-y^{k}\|^{2}+\frac{\eta^{2}}{\gamma}\|Fy^{k}-u^{k}\|^{2}-2\eta\langle
Fy^{k}-Fx^{\star},y^{k}-x^{\star}\rangle.\end{array}$ (25)

If, in addition, $F$ is monotone and $\beta:=1$, then for $\omega\geq 1$,
$\gamma>0$, and $t>0$, we have

$\begin{array}[]{lcl}\|w^{k+1}\|^{2}+(\omega-1)\|w^{k+1}-\hat{w}^{k+1}\|^{2}&\leq&\|w^{k}\|^{2}-(1-\gamma)\|w^{k}-\tilde{w}^{k}\|^{2}+\left[\frac{1}{\gamma}+\frac{\omega(1+t)L^{2}\eta^{2}}{t}\right]\|Fx^{k}-u^{k}\|^{2}\vspace{1ex}\\\
&&-{~{}}\big{[}1-\omega(1+t)L^{2}\eta^{2}\big{]}\|\hat{w}^{k+1}-\tilde{w}^{k}\|^{2}.\end{array}$
(26)

###### Proof.

Firstly, since $\xi^{k+1}\in Tx^{k+1}$, $\zeta^{k}\in Ty^{k}$, and
$\xi^{\star}=-Fx^{\star}\in Tx^{\star}$, by the maximally $3$-cyclic
monotonicity of $T$, we have
$\langle\xi^{k+1},x^{k+1}-x^{\star}\rangle+\langle\xi^{\star},x^{\star}-y^{k}\rangle+\langle\zeta^{k},y^{k}-x^{k+1}\rangle\geq
0$, leading to
$\langle\xi^{k+1}-\zeta^{k},x^{k+1}-x^{\star}\rangle\geq\langle\zeta^{k}-\xi^{\star},x^{\star}-y^{k}\rangle=-\langle
Fx^{\star}+\zeta^{k},y^{k}-x^{\star}\rangle$. Utilizing this inequality and
the second line $x^{k}-x^{k+1}=\eta(Fy^{k}+\xi^{k+1})$ of (23), for any
$x^{\star}\in\mathrm{zer}(\Phi)$, we can derive that

$\begin{array}[]{lcl}\|x^{k+1}-x^{\star}\|^{2}&=&\|x^{k}-x^{\star}\|^{2}-2\langle
x^{k}-x^{k+1},x^{k+1}-x^{\star}\rangle-\|x^{k+1}-x^{k}\|^{2}\vspace{1ex}\\\
&=&\|x^{k}-x^{\star}\|^{2}-2\eta\langle
Fy^{k}+\xi^{k+1},x^{k+1}-x^{\star}\rangle-\|x^{k+1}-x^{k}\|^{2}\vspace{1ex}\\\
&=&\|x^{k}-x^{\star}\|^{2}-2\eta\langle
Fy^{k}+\zeta^{k},x^{k+1}-x^{\star}\rangle-\|x^{k+1}-x^{k}\|^{2}-2\eta\langle\xi^{k+1}-\zeta^{k},x^{k+1}-x^{\star}\rangle\vspace{1ex}\\\
&\leq&\|x^{k}-x^{\star}\|^{2}-\|x^{k+1}-x^{k}\|^{2}-2\eta\langle
Fy^{k}+\zeta^{k},x^{k+1}-y^{k}\rangle-2\eta\langle
Fy^{k}-Fx^{\star},y^{k}-x^{\star}\rangle.\end{array}$

Next, from the first line of (23), we have
$\eta(Fy^{k}+\zeta^{k})=\beta(x^{k}-y^{k})+\eta(Fy^{k}-u^{k})$. Therefore, by
the Cauchy-Schwarz inequality and Young’s inequality, for any $\gamma>0$, we
can derive that

$\begin{array}[]{lcl}2\eta\langle
Fy^{k}+\zeta^{k},x^{k+1}-y^{k}\rangle&=&2\beta\langle
x^{k}-y^{k},x^{k+1}-y^{k}\rangle+2\eta\langle
Fy^{k}-u^{k},x^{k+1}-y^{k}\rangle\vspace{1ex}\\\
&\geq&\beta\left[\|x^{k}-y^{k}\|^{2}+\|x^{k+1}-y^{k}\|^{2}-\|x^{k+1}-x^{k}\|^{2}\right]-2\eta\|Fy^{k}-u^{k}\|\|x^{k+1}-y^{k}\|\vspace{1ex}\\\
&\geq&\beta\|x^{k}-y^{k}\|^{2}+(\beta-\gamma)\|x^{k+1}-y^{k}\|^{2}-\beta\|x^{k+1}-x^{k}\|^{2}-\frac{\eta^{2}}{\gamma}\|Fy^{k}-u^{k}\|^{2}.\end{array}$

Finally, substituting this inequality into the above estimate, we obtain (25).

To prove (30), we process as follows. Using again the $3$-cyclic monotonicity
of $T$ but with $\xi^{k}\in Tx^{k}$, we have
$\langle\xi^{k+1},x^{k+1}-x^{k}\rangle+\langle\xi^{k},x^{k}-y^{k}\rangle+\langle\zeta^{k},y^{k}-x^{k+1}\rangle\geq
0$. By the monotonicity of $F$, we get $\langle
Fx^{k+1}-Fx^{k},x^{k+1}-x^{k}\rangle\geq 0$. Summing up these inequalities and
using $w^{k}=Fx^{k}+\xi^{k}$ and $\tilde{w}^{k}:=Fx^{k}+\zeta^{k}$, we have

$\langle w^{k+1}-\tilde{w}^{k},x^{k+1}-x^{k}\rangle+\langle
w^{k}-\tilde{w}^{k},x^{k}-y^{k}\rangle\geq 0.$ (27)

From the second line of (EG2), we have
$x^{k+1}-x^{k}=-\eta(Fy^{k}+\xi^{k+1})=-\eta\hat{w}^{k+1}$. From the first
line of (EG2) and $\beta=1$, we also have
$x^{k}-y^{k}=\frac{\eta}{\beta}(\tilde{w}^{k}+u^{k}-Fx^{k})=\eta\tilde{w}^{k}+\eta(u^{k}-Fx^{k})$.
Substituting these expressions into (27), and using an elementary inequality
$2\langle z,s\rangle\leq\gamma\|s\|^{2}+\frac{\|z\|^{2}}{\gamma}$ for any
$\gamma>0$, we have

$\begin{array}[]{lcl}0&\leq&2\langle\tilde{w}^{k},\hat{w}^{k+1}\rangle-2\langle
w^{k+1},\hat{w}^{k+1}\rangle+2\langle
w^{k},\tilde{w}^{k}\rangle-2\|\tilde{w}^{k}\|^{2}+2\langle
w^{k}-\tilde{w}^{k},u^{k}-Fx^{k}\rangle\vspace{1ex}\\\
&=&\|w^{k}\|^{2}-\|w^{k+1}\|^{2}+\|w^{k+1}-\hat{w}^{k+1}\|^{2}-\|\hat{w}^{k+1}-\tilde{w}^{k}\|^{2}-\|w^{k}-\tilde{w}^{k}\|^{2}+2\langle
w^{k}-\tilde{w}^{k},Fx^{k}-u^{k}\rangle\vspace{1ex}\\\
&\leq&\|w^{k}\|^{2}-\|w^{k+1}\|^{2}+\|w^{k+1}-\hat{w}^{k+1}\|^{2}-\|\hat{w}^{k+1}-\tilde{w}^{k}\|^{2}-(1-\gamma)\|w^{k}-\tilde{w}^{k}\|^{2}+\frac{1}{\gamma}\|Fx^{k}-u^{k}\|^{2}.\end{array}$

This inequality leads to

$\begin{array}[]{lcl}\|w^{k+1}\|^{2}&\leq&\|w^{k}\|^{2}+\|w^{k+1}-\hat{w}^{k+1}\|^{2}+\frac{1}{\gamma}\|Fx^{k}-u^{k}\|^{2}-(1-\gamma)\|w^{k}-\tilde{w}^{k}\|^{2}-\|\hat{w}^{k+1}-\tilde{w}^{k}\|^{2}.\end{array}$

Now, by the $L$-Lipschitz continuity of $F$, (EG2), and Young’s inequality,
for any $t>0$, we have

$\begin{array}[]{lcl}\|w^{k+1}-\hat{w}^{k+1}\|^{2}&=&\|Fx^{k+1}-Fy^{k}\|^{2}\leq
L^{2}\|x^{k+1}-y^{k}\|^{2}=L^{2}\eta^{2}\|\hat{w}^{k+1}-\tilde{w}^{k}+Fx^{k}-u^{k}\|^{2}\vspace{1ex}\\\
&\leq&(1+t)L^{2}\eta^{2}\|\hat{w}^{k+1}-\tilde{w}^{k}\|^{2}+\frac{(1+t)L^{2}\eta^{2}}{t}\|Fx^{k}-u^{k}\|^{2}.\end{array}$

Multiplying this inequality by $\omega\geq 1$ and adding the last inequality,
we obtain (30). ∎

###  4.3 Unified convergence analysis

The following theorem proves the best-iterate and the last-iterate convergence
of (EG2).

######  Theorem 4.1.

Suppose that $\mathrm{zer}(\Phi)\neq\emptyset$, $F$ in (NI) is $L$-Lipschitz
continuous and satisfies $\langle Fx-Fx^{\star},x-x^{\star}\rangle\geq 0$ for
all $x\in\mathrm{dom}(F)$ and some $x^{\star}\in\mathrm{zer}(\Phi)$, and $T$
is maximally $3$-cyclically monotone. Let $\\{(x^{k},y^{k})\\}$ be generated
by (EG2). Then, the following statements hold.

  * $\mathrm{(a)}$

$($EG method$)$ If we choose $u^{k}:=Fx^{k}$ and $0<\eta<\frac{\beta}{L}$,
then we have

$\min_{1\leq l\leq
k+1}\|Fx^{l}+\xi^{l}\|^{2}\leq\frac{1}{k+1}\sum_{l=1}^{k+1}\|Fx^{l}+\xi^{l}\|^{2}\leq\frac{C_{0}\|x^{0}-x^{\star}\|^{2}}{k+1},\quad\xi^{l}\in
Tx^{l},$ (28)

where $C_{0}:=\frac{3+2L^{2}}{\eta^{2}(\beta-L\eta)}>0$. As a consequence,
$\\{\|x^{k}-x^{\star}\|\\}$ is nonincreasing and

$\lim_{k\to\infty}\|x^{k}-y^{k}\|=\lim_{k\to\infty}\|Fy^{k}+\zeta^{k}\|=\lim_{k\to\infty}\|Fx^{k}+\xi^{k}\|=0.$
(29)

Moreover, $\\{x^{k}\\}$ converges to $x^{\star}$, a solution of (NI).

  * $($Last-iterate convergence rate of EG$)$ If, in addition, $F$ is monotone, then we have

$\begin{array}[]{lcl}\|Fx^{k+1}+\xi^{k+1}\|^{2}\leq\|Fx^{k}+\xi^{k}\|^{2}\quad\text{and}\quad\|Fx^{k}+\xi^{k}\|\leq\frac{\sqrt{C_{0}}\|x^{0}-x^{\star}\|}{\sqrt{k}}.\end{array}$
(30)

Hence, we have a last-iterate convergence rate
$\|Fx^{k}+\xi^{k}\|=\mathcal{O}\big{(}1/\sqrt{k}\big{)}$ of
$\|Fx^{k}+\xi^{k}\|$, where $\xi^{k}\in Tx^{k}$.

  * $\mathrm{(b)}$

$($Past-EG method$)$ If we choose $u^{k}:=Fy^{k-1}$ with $y^{-1}:=x^{0}$ and
$0<\eta<\frac{\beta}{3L}$, then

$\min_{1\leq l\leq
k+1}\|Fx^{l}+\xi^{l}\|^{2}\leq\frac{1}{k+1}\sum_{l=1}^{k+1}\left[\|Fx^{l}+\xi^{l}\|^{2}+\psi\cdot\|x^{l}-y^{l-1}\|^{2}\right]\leq\frac{\hat{C}_{0}\|x^{0}-x^{\star}\|^{2}}{k+1},$
(31)

where $\hat{C}_{0}:=\frac{3+2L^{2}}{\eta^{2}(\beta-3L\eta)}>0$ and
$\psi:=\frac{L(3+2L^{2})}{2\eta(\beta-3L\eta)}$. Moreover, $\\{x^{k}\\}$
converges to $x^{\star}$ and (29) still holds.

  * $($Last-iterate convergence rate of Past-EG$)$ If, in addition, $F$ is monotone, then we have

$\begin{array}[]{ll}&\|Fx^{k+1}+\xi^{k+1}\|^{2}+\kappa\|Fx^{k+1}-Fy^{k}\|^{2}\leq\|Fx^{k}+\xi^{k}\|^{2}+\kappa\|Fx^{k}-Fy^{k-1}\|^{2}\vspace{1ex}\\\
\text{and}&\|Fx^{k}+\xi^{k}\|^{2}\leq\|Fx^{k}+\xi^{k}\|^{2}+\kappa\|Fx^{k}-Fy^{k-1}\|^{2}\leq\frac{\hat{C}_{0}\|x^{0}-x^{\star}\|^{2}}{m_{0}k},\end{array}$
(32)

where $\kappa:=\frac{1+9L^{2}\eta^{2}}{1-9L^{2}\eta^{2}}>0$ and
$m_{0}:=\max\left\\{\frac{\kappa}{\psi},1\right\\}$. Therefore, we have
$\|Fx^{k}+\xi^{k}\|=\mathcal{O}\big{(}1/\sqrt{k}\big{)}$.

In both cases $\mathrm{(a)}$ and $\mathrm{(b)}$, if $T$ is maximally monotone,
then we also have

$\min_{1\leq l\leq
k}\|G_{\eta\Phi}x^{l}\|=\mathcal{O}\left(\frac{1}{\sqrt{k}}\right)\quad\text{and}\quad\min_{1\leq
l\leq k}\|G_{\eta\Phi}y^{l}\|=\mathcal{O}\left(\frac{1}{\sqrt{k}}\right),$
(33)

where $G_{\eta\Phi}x:=\frac{1}{\eta}(x-J_{\eta T}(x-\eta Fx))$ is given by
(4). Moreover,
$\lim_{k\to\infty}\|G_{\eta\Phi}y^{k}\|=\lim_{k\to\infty}\|G_{\eta\Phi}x^{k}\|=0$.

######  Remark 4.1.

If $\beta=1$, then the condition $0<\eta<\frac{1}{L}$ in Part (a) is the same
as in classical EG methods [50]. Similarly, when $\beta=1$, the condition
$0<\eta\leq\frac{1}{3L}$ in Part (b) is the same as the one in [113]. It
remains open to establish both best-iterate and last-iterate convergence rates
of (EG2) under weak-Minty solution assumption or the co-hypomonotonicity of
$\Phi$.

###### Proof of Theorem 4.1.

(a) (EG method) First, since $u^{k}:=Fx^{k}$, by the $L$-Lipschitz continuity
of $F$, we have $\|Fy^{k}-u^{k}\|=\|Fy^{k}-Fx^{k}\|\leq L\|x^{k}-y^{k}\|$.
Using this inequality and $\langle
Fy^{k}-Fx^{\star},y^{k}-x^{\star}\rangle\geq 0$ from our assumption into (25),
we have

$\begin{array}[]{lcl}\|x^{k+1}-x^{\star}\|^{2}&\leq&\|x^{k}-x^{\star}\|^{2}-(1-\beta)\|x^{k+1}-x^{k}\|^{2}-(\beta-\gamma)\|x^{k+1}-y^{k}\|^{2}-\left(\beta-\frac{L^{2}\eta^{2}}{\gamma}\right)\|x^{k}-y^{k}\|^{2}.\end{array}$
(34)

Next, using the second line of (23), Young’s inequality, and the $L$-Lipschitz
continuity of $F$, we have

$\begin{array}[]{lcl}\eta^{2}\|Fx^{k+1}+\xi^{k+1}\|^{2}&\overset{\tiny\eqref{eq:EG4NI_reform}}{=}&\|x^{k+1}-x^{k}+\eta(Fx^{k+1}-Fy^{k})\|^{2}\vspace{1ex}\\\
&\leq&(1+\frac{2L^{2}}{3})\|x^{k+1}-x^{k}\|^{2}+(1+\frac{3}{2L^{2}})\|Fx^{k+1}-Fy^{k}\|^{2}\vspace{1ex}\\\
&\leq&\frac{3}{2}(1+\frac{2L^{2}}{3})\|x^{k+1}-y^{k}\|^{2}+(1+\frac{3}{2L^{2}})L^{2}\|x^{k+1}-y^{k}\|^{2}+3(1+\frac{2L^{2}}{3})\|x^{k}-y^{k}\|^{2}\vspace{1ex}\\\
&=&(3+2L^{2})[\|x^{k+1}-y^{k}\|^{2}+\|x^{k}-y^{k}\|^{2}].\end{array}$ (35)

Substituting (35) into (34) with $\gamma:=L\eta$, and assuming that
$L\eta<\beta$, we get

$\begin{array}[]{lcl}\|x^{k+1}-x^{\star}\|^{2}&\leq&\|x^{k}-x^{\star}\|^{2}-(\beta-L\eta)\big{[}\|x^{k+1}-y^{k}\|^{2}+\|x^{k}-y^{k}\|^{2}\big{]}\vspace{1ex}\\\
&\leq&\|x^{k}-x^{\star}\|^{2}-\frac{(\beta-L\eta)\eta^{2}}{3+2L^{2}}\|Fx^{k+1}+\xi^{k+1}\|^{2}.\end{array}$

Now, using this estimate, we can easily prove (28) in Theorem 4.1. Note that,
by (5), Young’s inequality, the $L$-Lipschitz continuity of $F$, and
$x^{k}-y^{k}=\frac{\eta}{\beta}(u^{k}+\zeta^{k})=\frac{\eta}{\beta}(Fx^{k}+\zeta^{k})$
from the first line of (23), we have

$\begin{array}[]{lcl}\|Fy^{k}+\zeta^{k}\|^{2}&\leq&\frac{3}{2}\|Fy^{k}-Fx^{k}\|^{2}+2\|Fx^{k}+\zeta^{k}\|^{2}\leq\frac{3(L^{2}\eta^{2}+2\beta^{2})}{2\eta^{2}}\|x^{k}-y^{k}\|^{2}.\end{array}$
(36)

Therefore, the remaining statements are direct consequences of (28), (36), and
(35) using standard arguments.

Finally, since $\beta=1$, using $u^{k}:=Fx^{k}$, $\omega:=1$, and $t=\gamma\to
0^{+}$ into (30), and noting that $L\eta\leq 1$, we get

$\begin{array}[]{lcl}\|w^{k+1}\|^{2}&\leq&\|w^{k}\|^{2}-(1-L^{2}\eta^{2})\|\hat{w}^{k+1}-\tilde{w}^{k}\|^{2}-\|w^{k}-\tilde{w}^{k}\|^{2}\leq\|w^{k}\|^{2}.\end{array}$

This shows that $\left\\{\|w^{k}\|\right\\}$ is monotonically nonincreasing.
Combining this property and (28), we obtain (30).

(b) (Past-extragradient method) If we choose $u^{k}:=Fy^{k-1}$, then by
Young’s inequality and the $L$-Lipschitz continuity of $F$, similar to the
proof of (17), we have

$\begin{array}[]{lcl}\|Fy^{k}-u^{k}\|^{2}=\|Fy^{k}-Fy^{k-1}\|^{2}\leq
3L^{2}\|x^{k}-y^{k}\|^{2}+\frac{3L^{2}}{2}\|x^{k}-y^{k-1}\|^{2}.\end{array}$

Substituting this expression into (25), using $\langle
Fy^{k}-Fx^{\star},y^{k}-x^{\star}\rangle\geq 0$, and choosing $\gamma:=L\eta$,
we obtain

$\begin{array}[]{lcl}\|x^{k+1}-x^{\star}\|^{2}+\frac{3L\eta}{2}\|x^{k+1}-y^{k}\|^{2}&\leq&\|x^{k}-x^{\star}\|^{2}+\frac{3L\eta}{2}\|x^{k}-y^{k-1}\|^{2}-(1-\beta)\|x^{k+1}-x^{k}\|^{2}\vspace{1ex}\\\
&&-{~{}}\frac{L\eta}{2}\|x^{k+1}-y^{k}\|^{2}-(\beta-3L\eta)\big{[}\|x^{k+1}-y^{k}\|^{2}+\|x^{k}-y^{k}\|^{2}\big{]}.\end{array}$
(37)

Assuming that $3L\eta<\beta$. Then, combining (37) and (35), we can easily
prove (31) in Theorem 4.1. Moreover, (37) also implies
$\lim_{k\to\infty}\|x^{k}-y^{k}\|=\lim_{k\to\infty}\|x^{k+1}-y^{k}\|=\lim_{k\to\infty}\|Fx^{k}+\xi^{k}\|=0$.
By the second line of (23), we have $\eta\|Fy^{k}+\zeta^{k}\|\leq
L\eta\|x^{k}-y^{k}\|+L\eta\|x^{k}-y^{k-1}\|+\|\eta(Fy^{k-1}+\zeta^{k}\|=(L\eta+1)\|x^{k}-y^{k}\|+L\eta\|x^{k}-y^{k-1}\|$.
Using this relation and
$\lim_{k\to\infty}\|x^{k}-y^{k}\|=\lim_{k\to\infty}\|x^{k+1}-y^{k}\|=0$, we
obtain $\lim_{k\to\infty}\|Fy^{k}+\zeta^{k}\|=0$. The convergence of
$\\{x^{k}\\}$ to $x^{\star}$ follows from standard arguments.

Next, assume that $3L\eta<1$ and $u^{k}:=Fy^{k-1}$. Then, substituting
$\gamma:=1$, $t:=\frac{1}{8}$, $\omega:=\frac{2}{1-9L^{2}\eta^{2}}>1$ into
(30), and using $\|Fx^{k}-u^{k}\|=\|Fx^{k}-Fy^{k-1}\|=\|w^{k}-\hat{w}^{k}\|$,
we obtain

$\begin{array}[]{lcl}\|w^{k+1}\|^{2}+\kappa\|w^{k+1}-\hat{w}^{k+1}\|^{2}&\leq&\|w^{k}\|^{2}+\kappa\|w^{k}-\hat{w}^{k}\|^{2}-\frac{4-27L^{2}\eta^{2}}{4(1-9L^{2}\eta^{2})}\|\hat{w}^{k+1}-\tilde{w}^{k}\|^{2}.\end{array}$

This is exactly first line of (32). Since $3L\eta<1$, we have
$\kappa:=\frac{1+9L^{2}\eta^{2}}{1-9L^{2}\eta^{2}}>0$ and
$\frac{4-27L^{2}\eta^{2}}{4(1-9L^{2}\eta^{2})}>0$. Let
$m_{0}:=\max\\{\frac{\kappa}{\psi},1\\}$, where $\psi$ is given in (31). Then,
we have

$\begin{array}[]{lcl}\|Fx^{k}+\xi^{k}\|^{2}+\kappa\|Fx^{k}-Fy^{k-1}\|^{2}&\leq&\|Fx^{k}+\xi^{k}\|^{2}+\kappa
L^{2}\|x^{k}-y^{k-1}\|^{2}\leq
m_{0}\big{[}\|Fx^{k}+\xi^{k}\|^{2}+\psi\cdot\|x^{k}-y^{k-1}\|^{2}\big{]}.\end{array}$

Combining this inequality and (31), we obtain
$\frac{1}{k+1}\sum_{l=1}^{k+1}\left[\|Fx^{l}+\xi^{l}\|^{2}+\kappa\|Fx^{l}-Fy^{l-1}\|^{2}\right]\leq\frac{\hat{C}_{0}\|x^{0}-x^{\star}\|^{2}}{m_{0}(k+1)}$.
This bound together with the first line of (32) imply the second line of (32).

Finally, since $T$ is maximally monotone, $J_{\eta T}$ is single-valued and
nonexpansive. By using $\|G_{\Phi}x^{k}\|\leq\|Fx^{k}+\xi^{k}\|$ from (5) and
either (28) or (31), we obtain (33). Using again (5) and the limits (29), we
obtain
$\lim_{k\to\infty}\|G_{\Phi}x^{k}\|\leq\lim_{k\to\infty}\|Fx^{k}+\xi^{k}\|=0$
and
$\lim_{k\to\infty}\|G_{\Phi}y^{k}\|\leq\lim_{k\to\infty}\|Fy^{k}+\zeta^{k}\|=0$.
∎

##  5 Forward-Backward-Forward Splitting-Type Methods for (NI)

Alternative to (EG2), we now survey recent results on the best-iterate
convergence rates of the FBFS method and its variants for solving (NI). As
before, we provide a unified analysis that covers a wide class of FBFS
variants, which can also solve (NI) under a weak-Minty solution and
particularly, the co-hypomonotonicity.

###  5.1 The class of forward-backward-forward splitting methods

The forward-backward-forward splitting (FBFS) method was proposed by P. Tseng
in [136] for solving (NI), which is originally called a modified forward-
backward splitting method. This method was developed to solve (NI) with
additional constraints. Instead of presenting the original scheme in [136], we
modify it using the idea of EG$+$ in [43] and combine two variants in one.
Starting from an initial point $x^{0}\in\mathrm{dom}(\Phi)$, at each iteration
$k\geq 0$, we update

$\left\\{\begin{array}[]{lcl}y^{k}&\in&J_{\frac{\eta}{\beta}T}(x^{k}-\frac{\eta}{\beta}u^{k}),\vspace{1ex}\\\
x^{k+1}&:=&\beta y^{k}+(1-\beta)x^{k}-\eta(Fy^{k}-u^{k}),\end{array}\right.$
(FBFS2)

where $\eta>0$ is a given stepsize, $\beta\in(0,1]$ is a scaling factor, and
$u^{k}$ is one of the following choices:

  * •

Option 1. If we choose $u^{k}:=Fx^{k}$, then we obtain a variant of Tseng’s
FBFS method. In particular, if $\beta=1$, then we get exactly Tseng’s FBFS
method in [136] for solving (NI). Note that one can extend (FBFS2) to cover
the case $\mathrm{zer}(\Phi)\cap\mathcal{C}\neq\emptyset$ for some subset
$\mathcal{C}$ of $\mathbb{R}^{p}$ as presented in [136]. Nevertheless, for
simplicity, we assume that $\mathcal{X}=\mathbb{R}^{p}$. As shown in (FBFS),
if $T=0$, then (FBFS2) reduces to the classical extragradient method (EG).
However, if $T\neq 0$, then (FBFS2) is different from (EG2).

  * •

Option 2. If we choose $u^{k}:=Fy^{k-1}$, where $y^{-1}:=x^{0}$, then we
obtain a past-FBFS variant. This variant can also be referred to as a
generalized variant of the optimistic gradient (OG) method, see, e.g., [38,
95, 96]. If $\beta=1$, then $y^{k+1}\in J_{\eta T}(x^{k+1}-\eta Fy^{k})$ and
$x^{k+1}=y^{k}-\eta(Fy^{k}-Fy^{k-1})$. Combining these two expressions,
(FBFS2) reduces to

$\begin{array}[]{lcl}y^{k+1}&\in&J_{\eta
T}\left(y^{k}-\eta(2Fy^{k}-Fy^{k-1})\right).\end{array}$ (FRBS2)

This is exactly the forward-reflected-backward splitting (FRBS) method in
[89].

Compared to (EG2), we do not require $T$ to be monotone in (FBFS2). However,
to guarantee the well-definedness of $\\{(x^{k},y^{k})\\}$, we need
$y^{k}\in\mathrm{ran}(J_{\eta T})$ and $y^{k}\in\mathrm{dom}(F)$. Hence, we
can assume that $\mathrm{ran}(J_{\eta
T})\subseteq\mathrm{dom}(F)=\mathbb{R}^{p}$ and $\mathrm{dom}(J_{\eta
T})=\mathbb{R}^{p}$. This requirement makes (FBFS2) cover a broader class of
problems than (EG2), and it obviously holds if $T$ is maximally monotone and
$F$ is monotone and Lipschitz continuous as in (EG2). In addition, (FBFS2)
only requires one evaluation of $J_{\eta T}$ instead of two as in (EG2),
reducing the per-iteration complexity when $J_{\eta T}$ is expensive to
evaluate.

Similar to (23), we can rewrite (FBFS2) equivalently to

$\left\\{\begin{array}[]{lcl}y^{k}&:=&x^{k}-\frac{\eta}{\beta}(u^{k}+\zeta^{k}),\quad\zeta^{k}\in
Ty^{k},\vspace{1ex}\\\
x^{k+1}&:=&x^{k}+\beta(y^{k}-x^{k})-\eta(Fy^{k}-u^{k}).\end{array}\right.$
(38)

This representation is an important step for our convergence analysis below.

###  5.2 One-iteration analysis

The following lemma provides a key estimate to establish convergence of
(FBFS2).

######  Lemma 5.1.

Suppose that $\\{(x^{k},y^{k})\\}$ is generated by (FBFS2) and $T$ is not
necessary monotone, but $\mathrm{ran}(J_{\eta
T})\subseteq\mathrm{dom}(F)=\mathbb{R}^{p}$ and $\mathrm{dom}(J_{\eta
T})=\mathbb{R}^{p}$. Then, for any $\gamma>0$, any
$x^{\star}\in\mathrm{zer}(\Phi)$, we have

$\begin{array}[]{lcl}\|x^{k+1}-x^{\star}\|^{2}&\leq&\|x^{k}-x^{\star}\|^{2}-(1-\beta)\|x^{k+1}-x^{k}\|^{2}-(\beta-\gamma)\|x^{k+1}-y^{k}\|^{2}\vspace{1ex}\\\
&&-{~{}}\beta\|x^{k}-y^{k}\|^{2}+\frac{\eta^{2}}{\gamma}\|Fy^{k}-u^{k}\|^{2}-2\eta\langle
Fy^{k}+\zeta^{k},y^{k}-x^{\star}\rangle.\end{array}$ (39)

###### Proof.

First, combining the first and second lines of (38), we obtain
$x^{k+1}=x^{k}-\beta(x^{k}-y^{k})+\eta(u^{k}-Fy^{k})=x^{k}-\eta(Fy^{k}+\zeta^{k})$.
Using this relation, we have

$\begin{array}[]{lcl}\|x^{k+1}-x^{\star}\|^{2}&=&\|x^{k}-x^{\star}\|^{2}-2\eta\langle
Fy^{k}+\zeta^{k},x^{k+1}-y^{k}\rangle-\|x^{k+1}-x^{k}\|^{2}-2\eta\langle
Fy^{k}+\zeta^{k},y^{k}-x^{\star}\rangle.\end{array}$

Next, from the first line of (38), we have
$\eta(Fy^{k}+\zeta^{k})=\beta(x^{k}-y^{k})+\eta(Fy^{k}-u^{k})$. Therefore, by
the Cauchy-Schwarz inequality and Young’s inequality, we can derive that

$\begin{array}[]{lcl}2\eta\langle
Fy^{k}+\zeta^{k},x^{k+1}-y^{k}\rangle&=&2\beta\langle
x^{k}-y^{k},x^{k+1}-y^{k}\rangle+2\eta\langle
Fy^{k}-u^{k},x^{k+1}-y^{k}\rangle\vspace{1ex}\\\
&\geq&\beta\left[\|x^{k}-y^{k}\|^{2}+\|x^{k+1}-y^{k}\|^{2}-\|x^{k+1}-x^{k}\|^{2}\right]-2\eta\|Fy^{k}-u^{k}\|\|x^{k+1}-y^{k}\|\vspace{1ex}\\\
&\geq&\beta\|x^{k}-y^{k}\|^{2}+(\beta-\gamma)\|x^{k+1}-y^{k}\|^{2}-\beta\|x^{k+1}-x^{k}\|^{2}-\frac{\eta^{2}}{\gamma}\|Fy^{k}-u^{k}\|^{2}.\end{array}$

Finally, substituting this inequality into the above estimate, we obtain (39).
∎

###  5.3 Unified convergence analysis

The following theorem establishes the best-iterate convergence rates of
(FBFS2) under star-co-hypomonotonicity.

######  Theorem 5.1.

Suppose that $F$ in (NI) is $L$-Lipschitz continuous, and
$\mathrm{zer}(\Phi)\neq\emptyset$. Suppose additionally that $\Phi$ is
$\rho$-star co-hypomonotone $($i.e. $\langle
u,x-x^{\star}\rangle\geq-\rho\|u\|^{2}$ for all $(x,u)\in\mathrm{gra}(\Phi)$
and for any $x^{\star}\in\mathrm{zer}(\Phi)$, where $\rho\geq 0$, $T$ is not
necessarily monotone, but $\mathrm{ran}(J_{\eta
T})\subseteq\mathrm{dom}(F)=\mathbb{R}^{p}$ and $\mathrm{dom}(J_{\eta
T})=\mathbb{R}^{p}$. Let $\\{(x^{k},y^{k})\\}$ be generated by (FBFS2). Then,
the following statements hold.

  * $\mathrm{(a)}$

$($FBFS method$)$ Let us choose $u^{k}:=Fx^{k}$ and assume that
$L\rho\leq\frac{3\sqrt{2}-2}{12}\approx 0.1869$. Then, for any
$\beta\in(0,1]$, if $\eta$ is chosen such that

$\begin{array}[]{lcl}0\leq\frac{\beta\big{[}1-\sqrt{1-24L\rho(3L\rho+1)}\big{]}}{2L(3L\rho+1)}<\eta<\frac{\beta\big{[}1+\sqrt{1-24L\rho(3L\rho+1)}\big{]}}{2L(3L\rho+1)}\leq\frac{\beta}{L},\end{array}$
(40)

then we have

$\min_{1\leq l\leq
k+1}\|Fx^{l}+\xi^{l}\|^{2}\leq\frac{1}{k+1}\sum_{l=1}^{k+1}\|Fx^{l}+\xi^{l}\|^{2}\leq\frac{C_{\rho}\|x^{0}-x^{\star}\|^{2}}{k+1},\quad\xi^{l}\in
Tx^{l},$ (41)

where
$C_{\rho}:=\frac{3+2L^{2}}{\eta(\beta\eta-6\beta^{2}\rho-(3L\rho+1)L\eta^{2})}>0$.
As a consequence, $\\{\|x^{k}-x^{\star}\|\\}$ is nonincreasing and

$\lim_{k\to\infty}\|x^{k}-y^{k}\|=\lim_{k\to\infty}\|Fy^{k}+\zeta^{k}\|=\lim_{k\to\infty}\|Fx^{k}+\xi^{k}\|=0.$
(42)

Moreover, $\\{x^{k}\\}$ converges to $x^{\star}\in\mathrm{zer}(\Phi)$, a
solution of (NI).

  * $\mathrm{(b)}$

$($Past-FBFS/OG method$)$ Let us choose $u^{k}:=Fy^{k-1}$ with $y^{-1}:=x^{0}$
and assume that $L\rho\leq\frac{2\sqrt{3}-3}{24}\approx 0.01934$. Then, for
any $\beta\in(0,1]$, if $\eta$ is chosen such that

$\begin{array}[]{lcl}0\leq\frac{\beta\big{[}1-\sqrt{1-48L\rho(4L\rho+1)}\big{]}}{6L(4L\rho+1)}<\eta<\frac{\beta\big{[}1+\sqrt{1-48L\rho(4L\rho+1)}\big{]}}{6L(4L\rho+1)}\leq\frac{\beta}{3L},\end{array}$
(43)

then we have

$\min_{1\leq l\leq
k+1}\|Fx^{l}+\xi^{l}\|^{2}\leq\frac{1}{k+1}\sum_{l=1}^{k+1}\left[\|Fx^{l}+\xi^{l}\|^{2}+\frac{L\eta^{2}+8\beta^{2}\rho}{2\eta}\|x^{l}-y^{l-1}\|^{2}\right]\leq\frac{\hat{C}_{\rho}\|x^{0}-x^{\star}\|^{2}}{k+1},$
(44)

where
$\hat{C}_{\rho}:=\frac{3+2L^{2}}{\eta(\beta\eta-4\beta^{2}\rho-3(4L\rho+1)L\eta^{2})}>0$.
Moreover, $\\{x^{k}\\}$ converges to $x^{\star}\in\mathrm{zer}(\Phi)$ and (29)
still holds.

In both cases $\mathrm{(a)}$ and $\mathrm{(b)}$, if $J_{\eta T}$ is single-
valued and nonexpansive, then

$\min_{1\leq l\leq
k}\|G_{\eta\Phi}x^{l}\|=\mathcal{O}\left(\frac{1}{\sqrt{k}}\right)\quad\text{and}\quad\min_{1\leq
l\leq k}\|G_{\eta\Phi}y^{l}\|=\mathcal{O}\left(\frac{1}{\sqrt{k}}\right),$
(45)

where $G_{\eta\Phi}x:=\frac{1}{\eta}(x-J_{\eta T}(x-\eta Fx))$ is given by
(4). Moreover,
$\lim_{k\to\infty}\|G_{\eta\Phi}y^{k}\|=\lim_{k\to\infty}\|G_{\eta\Phi}x^{k}\|=0$.

###### Proof.

(a) (FBFS method) For $u^{k}:=Fx^{k}$, by the $L$-Lipschitz continuity of $F$,
we have $\|Fy^{k}-u^{k}\|=\|Fy^{k}-Fx^{k}\|\leq L\|y^{k}-x^{k}\|$. Using this
inequality into (39), we have

$\begin{array}[]{lcl}\|x^{k+1}-x^{\star}\|^{2}&\leq&\|x^{k}-x^{\star}\|^{2}-(1-\beta)\|x^{k+1}-x^{k}\|^{2}-(\beta-\gamma)\|x^{k+1}-y^{k}\|^{2}\vspace{1ex}\\\
&&-{~{}}\left(\beta-\frac{L^{2}\eta^{2}}{\gamma}\right)\|x^{k}-y^{k}\|^{2}-2\eta\langle
Fy^{k}+\zeta^{k},y^{k}-x^{\star}\rangle.\end{array}$ (46)

Now, by Young’s inequality, the $L$-Lipschitz continuity of $F$, and
$x^{k}-y^{k}=\frac{\eta}{\beta}(u^{k}+\zeta^{k})=\frac{\eta}{\beta}(Fx^{k}+\zeta^{k})$
from the first line of (38), we have

$\begin{array}[]{lcl}\|Fy^{k}+\zeta^{k}\|^{2}&\leq&\frac{3}{2}\|Fy^{k}-Fx^{k}\|^{2}+2\|Fx^{k}+\zeta^{k}\|^{2}\leq\frac{3(L^{2}\eta^{2}+2\beta^{2})}{2\eta^{2}}\|x^{k}-y^{k}\|^{2}.\end{array}$
(47)

Similarly, using the second line of (38), Young’s inequality, and the
$L$-Lipschitz continuity of $F$, we have

$\begin{array}[]{lcl}\eta^{2}\|Fx^{k+1}+\xi^{k+1}\|^{2}&=&\|x^{k+1}-x^{k}+\eta(Fx^{k+1}-Fy^{k})\|^{2}\vspace{1ex}\\\
&\leq&(1+\frac{2L^{2}}{3})\|x^{k+1}-x^{k}\|^{2}+(1+\frac{3}{2L^{2}})\|Fx^{k+1}-Fy^{k}\|^{2}\vspace{1ex}\\\
&\leq&\frac{3}{2}(1+\frac{2L^{2}}{3})\|x^{k+1}-y^{k}\|^{2}+(1+\frac{3}{2L^{2}})L^{2}\|x^{k+1}-y^{k}\|^{2}+3(1+\frac{2L^{2}}{3})\|x^{k}-y^{k}\|^{2}\vspace{1ex}\\\
&=&(3+2L^{2})[\|x^{k+1}-y^{k}\|^{2}+\|x^{k}-y^{k}\|^{2}].\end{array}$ (48)

Next, since $\langle u,x-x^{\star}\rangle\geq-\rho\|u\|^{2}$ for any
$(x,u)\in\mathrm{gra}(\Phi)$, we have $\langle
Fy^{k}+\zeta^{k},y^{k}-x^{\star}\rangle\geq-\rho\|Fy^{k}+\zeta^{k}\|^{2}$
since $(y^{k},\zeta^{k})\in\mathrm{gra}(T)$. Using this relation and (47) into
(46) with $\gamma:=L\eta$, we get

$\begin{array}[]{lcl}\|x^{k+1}-x^{\star}\|^{2}&\leq&\|x^{k}-x^{\star}\|^{2}-(\beta-L\eta)\|x^{k}-y^{k}\|^{2}-(\beta-L\eta)\|x^{k+1}-y^{k}\|^{2}+2\eta\rho\|Fy^{k}+\zeta^{k}\|^{2}\vspace{1ex}\\\
&\leq&\|x^{k}-x^{\star}\|^{2}-(\beta-L\eta)\|x^{k+1}-y^{k}\|^{2}-\tfrac{\beta\eta-6\beta^{2}\rho-(3L\rho+1)L\eta^{2}}{\eta}\|x^{k}-y^{k}\|^{2}.\end{array}$
(49)

Let us impose $\beta\eta-6\beta^{2}\rho-(3L\rho+1)L\eta^{2}>0$, which holds if

$\begin{array}[]{lcl}0\leq\frac{\beta\big{[}1-\sqrt{1-24L\rho(3L\rho+1)}\big{]}}{2L(3L\rho+1)}<\eta<\frac{\beta\big{[}1+\sqrt{1-24L\rho(3L\rho+1)}\big{]}}{2L(3L\rho+1)}\leq\frac{\beta}{L},\end{array}$

provided that $L\rho\leq\frac{3\sqrt{2}-2}{12}\approx 0.1869$. This choice of
$\eta$ is exactly (43). Let
$\psi:=\tfrac{\beta\eta-6\beta^{2}\rho-(3L\rho+1)L\eta^{2}}{\eta}=\beta-L\eta-\frac{3\rho}{\eta}(2\beta^{2}+L^{2}\eta^{2})>0$.
Then, we have $\beta-L\eta-\psi\geq 0$. In this case, (49) becomes

$\begin{array}[]{lcl}\|x^{k+1}-x^{\star}\|^{2}&\leq&\|x^{k}-x^{\star}\|^{2}-\psi[\|x^{k}-y^{k}\|^{2}+\|x^{k+1}-y^{k}\|^{2}]-(\beta-L\eta-\psi)\|x^{k+1}-y^{k}\|^{2}\vspace{1ex}\\\
&\leq&\|x^{k}-x^{\star}\|^{2}-\psi[\|x^{k}-y^{k}\|^{2}+\|x^{k+1}-y^{k}\|^{2}].\end{array}$

Finally, using this estimate, (48), we can easily prove (41) in Theorem 5.1.
The remaining statements are direct consequences of (41), (47), and (48) using
standard arguments.

(b) (Past-FBFS/OG method) For $u^{k}:=Fy^{k-1}$, using Young’s inequality and
the $L$-Lipschitz continuity of $F$, we can derive

$\begin{array}[]{lcl}\|Fy^{k}-u^{k}\|^{2}=\|Fy^{k}-Fy^{k-1}\|^{2}\leq
3L^{2}\|x^{k}-y^{k}\|^{2}+\frac{3L^{2}}{2}\|x^{k}-y^{k-1}\|^{2}.\end{array}$

Next, since $\langle u,x-x^{\star}\rangle\geq-\rho\|u\|^{2}$ for all
$(x,u)\in\mathrm{gra}(\Phi)$, we have $\langle
Fy^{k}+\zeta^{k},y^{k}-x^{\star}\rangle\geq-\rho\|Fy^{k}+\zeta^{k}\|^{2}$.
Using this inequality and the first line of (38) as
$\eta(Fy^{k}+\zeta^{k})=\beta(x^{k}-y^{k})+\eta(Fy^{k}-Fy^{k-1})$.

$\begin{array}[]{lcl}\langle
Fy^{k}+\zeta^{k},y^{k}-x^{\star}\rangle&\geq&-\rho\|Fy^{k}+\zeta^{k}\|^{2}=-\frac{\rho}{\eta^{2}}\|\beta(x^{k}-y^{k})+\eta(Fy^{k}-Fy^{k-1})\|^{2}\vspace{1ex}\\\
&\geq&-\frac{2\rho\beta^{2}}{\eta^{2}}\|x^{k}-y^{k}\|^{2}-4\rho\|Fy^{k}-Fx^{k}\|^{2}-4\rho\|Fx^{k}-Fy^{k-1}\|^{2}\vspace{1ex}\\\
&\geq&-\frac{2\rho(\beta^{2}+2L^{2}\eta^{2})}{\eta^{2}}\|x^{k}-y^{k}\|^{2}-4\rho
L^{2}\|x^{k}-y^{k-1}\|^{2}.\end{array}$

Substituting the last two expressions into (39) and choosing $\gamma:=L\eta$,
we obtain

$\begin{array}[]{lcl}\mathcal{V}_{k+1}&\leq&\mathcal{V}_{k}-(1-\beta)\|x^{k+1}-x^{k}\|^{2}-\left(\beta-\frac{5L\eta}{2}-8\rho
L^{2}\eta\right)\|x^{k+1}-y^{k}\|^{2}\vspace{1ex}\\\
&&-{~{}}\left[\beta-3L\eta-\frac{4\rho(\beta^{2}+2L^{2}\eta^{2})}{\eta}\right]\|x^{k}-y^{k}\|^{2},\end{array}$
(50)

where $\mathcal{V}_{k}:=\|x^{k}-x^{\star}\|^{2}+L\eta\left(\frac{3}{2}+8\rho
L\right)\|x^{k}-y^{k-1}\|^{2}$.

Let us define
$\hat{\psi}:=\beta-3L\eta-\frac{4\rho(\beta^{2}+2L^{2}\eta^{2})}{\eta}$ and
$\hat{\varphi}:=\beta-\frac{5L\eta}{2}-8\rho L^{2}\eta$ and impose the
condition that $\hat{\psi}>0$. Then, it is clear that
$\hat{\varphi}-\hat{\psi}=\frac{L\eta}{2}+\frac{4\beta^{2}\rho}{\eta}\geq 0$.
Moreover, (50) becomes
$\mathcal{V}_{k+1}\leq\mathcal{V}_{k}-\hat{\psi}\big{[}\|x^{k+1}-y^{k}\|^{2}+\|x^{k}-y^{k}\|^{2}\big{]}-(\hat{\varphi}-\hat{\psi})\|x^{k+1}-y^{k}\|^{2}$.
Combining this estimate and (48), we can easily prove (44) in Theorem 5.1.
Note that the condition $\hat{\psi}>0$ holds if $\eta$ is chosen as in (43)
provided that $L\rho\leq\frac{2\sqrt{3}-3}{24}\approx 0.01934$. Using (44), we
can easily prove that $\\{x^{k}\\}$ converges to $x^{\star}$. Moreover, (44)
also implies
$\lim_{k\to\infty}\|x^{k}-y^{k}\|=\lim_{k\to\infty}\|x^{k+1}-y^{k}\|=\lim_{k\to\infty}\|Fx^{k}+\xi^{k}\|=0$.
By the second line of (38), we have

$\begin{array}[]{lcl}\eta\|Fy^{k}+\zeta^{k}\|\leq
L\eta\|x^{k}-y^{k}\|+L\eta\|x^{k}-y^{k-1}\|+\|\eta(Fy^{k-1}+\zeta^{k}\|=(L\eta+1)\|x^{k}-y^{k}\|+L\eta\|x^{k}-y^{k-1}\|.\end{array}$

Combing this inequality and
$\lim_{k\to\infty}\|x^{k}-y^{k}\|=\lim_{k\to\infty}\|x^{k+1}-y^{k}\|=0$, we
obtain $\lim_{k\to\infty}\eta\|Fy^{k}+\zeta^{k}\|=0$.

Finally, if $J_{\eta T}$ is single-valued and nonexpansive, then by using
$\|G_{\Phi}x^{k}\|\leq\|Fx^{k}+\xi^{k}\|$ from (5) and either (41) or (44), we
obtain (45). Using again (5) and the limits (42), we obtain
$\lim_{k\to\infty}\|G_{\Phi}x^{k}\|\leq\lim_{k\to\infty}\|Fx^{k}+\xi^{k}\|=0$
and
$\lim_{k\to\infty}\|G_{\Phi}y^{k}\|\leq\lim_{k\to\infty}\|Fy^{k}+\zeta^{k}\|=0$.
∎

######  Remark 5.1.

The $\rho$-star co-hypomonotone condition in Theorem 5.1(b) trivially holds if
$\Phi:=F+T$ is $\rho$-co-hypomonotone. Hence, the star co-hypomonotone
condition is generally weaker than the $\rho$-co-hypomonotonicity. Note that
we have not tried to optimize the parameters in Theorem 5.1, and generally in
the whole paper. By careful tightening bounds and selecting parameters in our
analysis (e.g., where Young’s inequality is used), we can improve the range of
parameters. Note that the convergence analysis for the monotone case is very
classical, which can be found, e.g., in [50, 73]. However, the convergence
rates for $\beta<1$, and for the star co-hypomonotone case are recent results.
The best-iterate convergence rate of (FBFS2) was proven in [83] for the star
co-hypomonotone case, but using a potential function. Here, we provide a
different proof using classical results in [50] combining with the star co-
hypomonotonicity of $\Phi$.

##  6 Two Other Variants of The Extragradient Method

In this section, we review two additional methods: the reflected forward-
backward splitting (RFBS) algorithm [30, 87] and the golden ratio (GR) scheme
[88]. The last-iterate analysis for RFBS scheme was recently given in [27],
but only for (VIP). Here, we provide a new analysis for both the best-iterate
and the last-iterate rates for RFBS to solve (NI), which is more general than
(VIP). The best-iterate convergence analysis for GR is modified the proof from
[88] to expand the range of parameters. Nevertheless, the last-iterate
convergence rate analysis of GR is still open.

###  6.1 Reflected forward-backward splitting method

The reflected forward-backward splitting method was proposed by Malitsky in
[87] to solve (VIP) and it is called the projected reflected gradient method.
It was generalized to solve monotone (NI) in [30], which is called the
reflected forward-backward splitting (RFBS) scheme. The last iterate
convergence rate of the projected reflected gradient method for (VIP) was
recently proven in [27]. In this subsection, we survey this method for solving
(NI). We provide a new best-iterate convergence rate analysis compared to
[30]. We also present an elementary proof for the last-iterate convergence
rate of RFBS for solving monotone (NI).

The reflected forward-backward splitting (RFBS) method to approximate a
solution of (NI) is described as follows. Starting from
$x^{0}\in\mathrm{dom}(\Phi)$, we choose $x^{-1}:=x^{0}$ and at each iteration
$k\geq 0$, we update

$\left\\{\begin{array}[]{lcl}y^{k}&:=&2x^{k}-x^{k-1},\vspace{1ex}\\\
x^{k+1}&:=&J_{\eta T}(x^{k}-\eta Fy^{k}),\end{array}\right.$ (RFBS2)

where $\eta>0$ is a given step-size, determined later. Clearly, if we
eliminate $y^{k}$, then (RFBS2) can be written as

$x^{k+1}:=J_{\eta T}(x^{k}-\eta F(2x^{k}-x^{k-1})).$

From the second line of (RFBS2), we have $\xi^{k+1}:=\frac{1}{\eta}(x^{k}-\eta
Fy^{k}-x^{k+1})\in Tx^{k+1}$. As before, if we denote

$w^{k}:=Fx^{k}+\xi^{k}\quad\text{and}\quad\hat{w}^{k}:=Fy^{k-1}+\xi^{k},$ (51)

then we can rewrite (RFBS2) equivalently to

$\left\\{\begin{array}[]{lclcl}y^{k}&:=&x^{k}+x^{k}-x^{k-1}&=&x^{k}-\eta\hat{w}^{k},\vspace{1ex}\\\
x^{k+1}&:=&x^{k}-\eta\hat{w}^{k+1}.\end{array}\right.$ (52)

This expression leads to $x^{k+1}-y^{k}=-\eta(\hat{w}^{k+1}-\hat{w}^{k})$.
Next, we prove the following lemmas for our analysis.

######  Lemma 6.1.

Assume that $T$ in (NI) is maximally monotone and $F$ in (NI) is $L$-Lipschitz
continuous and satisfies $\langle Fx-Fx^{\star},x-x^{\star}\rangle\geq 0$ for
all $x\in\mathrm{dom}(\Phi)$ and some $x^{\star}\in\mathrm{zer}(\Phi)$. Let
$\left\\{(x^{k},y^{k})\right\\}$ be generated by (RFBS2) using $\eta>0$ and
$\mathcal{V}_{k}$ be defined as

$\begin{array}[]{lcl}\mathcal{V}_{k}&:=&\|x^{k}-x^{\star}\|^{2}+2\|x^{k}-x^{k-1}\|^{2}+\big{(}1-\sqrt{2}L\eta\big{)}\|x^{k}-y^{k-1}\|^{2}+2\eta\langle
Fy^{k-1}-Fx^{\star},x^{k}-x^{k-1}\rangle.\end{array}$ (53)

Then, we have

$\begin{array}[]{lcl}\mathcal{V}_{k}&\geq&\mathcal{V}_{k+1}+\big{[}1-(1+\sqrt{2})L\eta\big{]}\left[\|y^{k}-x^{k}\|^{2}+\|x^{k}-y^{k-1}\|^{2}\right],\vspace{1.5ex}\\\
\mathcal{V}_{k}&\geq&(1-L\eta)\|x^{k}-x^{\star}\|^{2}+(1-(1+\sqrt{2})L\eta)\|x^{k}-y^{k-1}\|^{2}+2(1-L\eta)\|x^{k}-x^{k-1}\|^{2}.\end{array}$
(54)

###### Proof.

First, since (RFBS2) is equivalent to (52), we have
$x^{k+1}-x^{k}=-\eta\hat{w}^{k+1}$ from the second line of (52). Using this
expression, for any $x^{\star}\in\mathrm{zer}(\Phi)$, we have

$\begin{array}[]{lcl}\|x^{k+1}-x^{\star}\|^{2}&=&\|x^{k}-x^{\star}\|^{2}-2\eta\langle\hat{w}^{k+1},x^{k+1}-x^{\star}\rangle-\|x^{k+1}-x^{k}\|^{2}.\end{array}$
(55)

Next, since $Fx^{\star}+\xi^{\star}=0$ from the fact that
$x^{\star}\in\mathrm{zer}(\Phi)$ and $T$ is monotone, we have
$\langle\xi^{k+1},x^{k+1}-x^{\star}\rangle\geq\langle\xi^{\star},x^{k+1}-x^{\star}\rangle=-\langle
Fx^{\star},x^{k+1}-x^{\star}\rangle$, where $\xi^{\star}\in Tx^{\star}$. Using
this relation, we can prove that

$\begin{array}[]{lcl}\langle\hat{w}^{k+1},x^{k+1}-x^{\star}\rangle&=&\langle
Fy^{k},x^{k+1}-x^{\star}\rangle+\langle\xi^{k+1},x^{k+1}-x^{\star}\rangle\vspace{1ex}\\\
&\geq&\langle Fy^{k}-Fx^{\star},x^{k+1}-y^{k}\rangle+\langle
Fy^{k}-Fx^{\star},y^{k}-x^{\star}\rangle.\end{array}$ (56)

Utilizing $y^{k}-x^{k}=x^{k}-x^{k-1}$ from the first line of (RFBS2), we can
further expand

$\begin{array}[]{lcl}\langle Fy^{k}-Fx^{\star},y^{k}-x^{k+1}\rangle&=&\langle
Fy^{k-1}-Fx^{\star},y^{k}-x^{k}\rangle-\langle
Fy^{k}-Fx^{\star},x^{k+1}-x^{k}\rangle+\langle
Fy^{k}-Fy^{k-1},y^{k}-x^{k}\rangle\vspace{1ex}\\\ &=&\langle
Fy^{k-1}-Fx^{\star},x^{k}-x^{k-1}\rangle-\langle
Fy^{k}-Fx^{\star},x^{k+1}-x^{k}\rangle\vspace{1ex}\\\ &&+{~{}}\langle
Fy^{k}-Fy^{k-1},y^{k}-x^{k}\rangle.\end{array}$ (57)

Now, from the second line of (52), we have $\eta\xi^{k+1}=x^{k}-x^{k+1}-\eta
Fy^{k}$, leading to
$\eta(\xi^{k+1}-\xi^{k})=2x^{k}-x^{k-1}-x^{k+1}-\eta(Fy^{k}-Fy^{k-1})=y^{k}-x^{k+1}-\eta(Fy^{k}-Fy^{k-1})$.
By the monotonicity of $T$, we have $\langle
y^{k}-x^{k+1}-\eta(Fy^{k}-Fy^{k-1}),x^{k+1}-x^{k}\rangle=\eta\langle\xi^{k+1}-\xi^{k},x^{k+1}-x^{k}\rangle\geq
0$, leading to

$\begin{array}[]{lcl}2\eta\langle
Fy^{k}-Fy^{k-1},x^{k+1}-x^{k}\rangle&\leq&2\langle
y^{k}-x^{k+1},x^{k+1}-x^{k}\rangle\vspace{1ex}\\\
&=&\|y^{k}-x^{k}\|^{2}-\|x^{k+1}-x^{k}\|^{2}-\|x^{k+1}-y^{k}\|^{2}.\end{array}$

By the Cauchy-Schwarz inequality, the Lipschitz continuity of $F$, and Young’s
inequality, we can show that

$\begin{array}[]{lcl}2\eta\langle
Fy^{k}-Fy^{k-1},y^{k}-x^{k+1}\rangle&\leq&2\eta
L\|y^{k}-y^{k-1}\|\|x^{k+1}-y^{k}\|\vspace{1ex}\\\ &\leq&\eta
L(\sqrt{2}+1)\|y^{k}-x^{k}\|^{2}+\eta L\|x^{k}-y^{k-1}\|^{2}+\sqrt{2}\eta
L\|x^{k+1}-y^{k}\|^{2}.\end{array}$

Summing up the last two inequalities, we get

$\begin{array}[]{lcl}2\eta\langle
Fy^{k}-Fy^{k-1},y^{k}-x^{k}\rangle&\leq&\left[1+\eta
L(\sqrt{2}+1)\right]\|y^{k}-x^{k}\|^{2}+\eta
L\|x^{k}-y^{k-1}\|^{2}\vspace{1ex}\\\
&&-{~{}}(1-\sqrt{2}L\eta)\|x^{k+1}-y^{k}\|^{2}-\|x^{k+1}-x^{k}\|^{2}.\end{array}$

Combining this inequality, (56), and (57) with (55), we can eventually prove
that

$\begin{array}[]{lcl}\|x^{k+1}-x^{\star}\|^{2}&=&\|x^{k}-x^{\star}\|^{2}+2\eta\langle\hat{w}^{k+1},x^{\star}-x^{k+1}\rangle-\|x^{k+1}-x^{k}\|^{2}\vspace{1ex}\\\
&\leq&\|x^{k}-x^{\star}\|^{2}+2\eta\langle
Fy^{k-1}-Fx^{\star},x^{k}-x^{k-1}\rangle-2\eta\langle
Fy^{k}-Fx^{\star},x^{k+1}-x^{k}\rangle\vspace{1ex}\\\
&&+{~{}}\left[1+L\eta(\sqrt{2}+1)\right]\|y^{k}-x^{k}\|^{2}-2\|x^{k+1}-x^{k}\|^{2}\vspace{1ex}\\\
&&+{~{}}L\eta\|x^{k}-y^{k-1}\|^{2}-(1-\sqrt{2}L\eta)\|x^{k+1}-y^{k}\|^{2}-2\eta\langle
Fy^{k}-Fx^{\star},y^{k}-x^{\star}\rangle.\end{array}$

Using the definition (53) of $\mathcal{V}_{k}$, $\langle
Fy^{k}-Fx^{\star},y^{k}-x^{\star}\rangle\geq 0$, and
$y^{k}-x^{k}=x^{k}-x^{k-1}$, this estimate becomes

$\begin{array}[]{lcl}\mathcal{V}_{k+1}&:=&\|x^{k+1}-x^{\star}\|^{2}+2\|x^{k+1}-x^{k}\|^{2}+(1-\sqrt{2}L\eta)\|x^{k+1}-y^{k}\|^{2}+2\eta\langle
Fy^{k}-Fx^{\star},x^{k+1}-x^{k}\rangle\vspace{1ex}\\\
&\leq&\|x^{k}-x^{\star}\|^{2}+\left[1+(\sqrt{2}+1)L\eta\right]\|y^{k}-x^{k}\|^{2}+L\eta\|x^{k}-y^{k-1}\|^{2}+2\eta\langle
Fy^{k-1}-Fx^{\star},x^{k}-x^{k-1}\rangle\vspace{1ex}\\\
&=&\mathcal{V}_{k}-\left[1-(\sqrt{2}+1)L\eta\right]\left(\|y^{k}-x^{k}\|^{2}+\|x^{k}-y^{k-1}\|^{2}\right),\end{array}$

which proves the first inequality of (54).

Next, using Young’s inequality twice, we can show that

$\begin{array}[]{lcl}2\eta\langle
Fy^{k-1}-Fx^{\star},x^{k}-x^{k-1}\rangle&\geq&-\frac{L\eta}{2}\|y^{k-1}-x^{\star}\|^{2}-2L\eta\|x^{k}-x^{k-1}\|^{2}\vspace{1ex}\\\
&\geq&-L\eta\|x^{k}-x^{\star}\|^{2}-L\eta\|x^{k}-y^{k-1}\|^{2}-2L\eta\|x^{k}-x^{k-1}\|^{2}.\end{array}$

Substituting this estimate into (53), we get

$\begin{array}[]{lcl}\mathcal{V}_{k}&:=&\|x^{k}-x^{\star}\|^{2}+2\|x^{k}-x^{k-1}\|^{2}+\left(1-\sqrt{2}L\eta\right)\|x^{k}-y^{k-1}\|^{2}+2\eta\langle
Fy^{k-1}-Fx^{\star},x^{k}-x^{k-1}\rangle\vspace{1ex}\\\
&\geq&(1-L\eta)\|x^{k}-x^{\star}\|^{2}+(1-(1+\sqrt{2})L\eta)\|x^{k}-y^{k-1}\|^{2}+2(1-L\eta)\|x^{k}-x^{k-1}\|^{2},\end{array}$

which proves the second line of (54). ∎

######  Lemma 6.2.

Assume that $F+T$ in (NI) is maximally monotone and $F$ in (NI) is
$L$-Lipschitz continuous. Let $\left\\{(x^{k},y^{k})\right\\}$ be generated by
(RFBS2) using $\eta>0$. Then, we have

$\begin{array}[]{lcl}&\|Fx^{k}+\xi^{k}\|^{2}\leq\frac{5L^{2}\eta^{2}+3}{3\eta^{2}}\|x^{k}-y^{k}\|^{2}+\frac{5L^{2}\eta^{2}+3}{5\eta^{2}}\|x^{k}-y^{k-1}\|^{2}.\end{array}$
(58)

Moreover, if $\sqrt{2}L\eta<1$, then with
$\kappa:=\frac{2L^{2}\eta^{2}}{1-2L^{2}\eta^{2}}>0$, we also have

$\begin{array}[]{lcl}\|Fx^{k+1}+\xi^{k+1}\|^{2}+\kappa\|Fx^{k+1}-Fy^{k}\|^{2}&\leq&\|Fx^{k}+\xi^{k}\|^{2}+\kappa\|Fx^{k}-Fy^{k-1}\|^{2}\vspace{1ex}\\\
&&-{~{}}\left(\frac{1-4L^{2}\eta^{2}}{1-2L^{2}\eta^{2}}\right)\|Fy^{k}-Fx^{k}+\xi^{k+1}-\xi^{k}\|^{2}.\end{array}$
(59)

###### Proof.

First, by Young’s inequality, $x^{k}-y^{k}=-\eta(Fy^{k-1}+\xi^{k})$ from (52),
and the $L$-Lipschitz continuity of $F$, we have

$\begin{array}[]{lcl}\|Fx^{k}+\xi^{k}\|^{2}&\leq&\big{(}1+\frac{5L^{2}\eta^{2}}{3}\big{)}\|Fy^{k-1}+\xi^{k}\|^{2}+\left(1+\frac{3}{5L^{2}\eta^{2}}\right)\|Fx^{k}-Fy^{k-1}\|^{2}\vspace{1ex}\\\
&\leq&\frac{5L^{2}\eta^{2}+3}{3\eta^{2}}\|x^{k}-y^{k}\|^{2}+\frac{5L^{2}\eta^{2}+3}{5\eta^{2}}\|x^{k}-y^{k-1}\|^{2}.\end{array}$

This estimate is exactly the first line of (58).

Next, by the monotonicity of $F+T$, we have $\langle
w^{k+1}-w^{k},x^{k+1}-x^{k}\rangle\geq 0$. Substituting
$x^{k+1}-x^{k}=-\eta\hat{w}^{k+1}$ into this inequality, we get

$\begin{array}[]{lcl}0&\leq&2\langle w^{k},\hat{w}^{k+1}\rangle-2\langle
w^{k+1},\hat{w}^{k+1}\rangle=\|w^{k}\|^{2}-\|w^{k+1}\|+\|w^{k+1}-\hat{w}^{k+1}\|^{2}-\|\hat{w}^{k+1}-w^{k}\|^{2}.\end{array}$

This inequality implies that

$\begin{array}[]{lcl}\|w^{k+1}\|^{2}&\leq&\|w^{k}\|^{2}+\|w^{k+1}-\hat{w}^{k+1}\|^{2}-\|\hat{w}^{k+1}-w^{k}\|^{2}.\end{array}$

On the other hand, by the Lipschitz continuity of $F$ and
$x^{k+1}-y^{k}=-\eta(\hat{w}^{k+1}-\hat{w}^{k})$, we have
$\|w^{k+1}-\hat{w}^{k+1}\|^{2}=\|Fx^{k+1}-Fy^{k}\|^{2}\leq
L^{2}\|x^{k+1}-y^{k}\|^{2}=L^{2}\eta^{2}\|\hat{w}^{k+1}-\hat{w}^{k}\|^{2}\leq
2L^{2}\eta^{2}\|\hat{w}^{k+1}-w^{k}\|^{2}+2L^{2}\eta^{2}\|w^{k}-\hat{w}^{k}\|^{2}$.
Multiplying this inequality by $\omega\geq 1$, and adding the result to the
last inequality, we get

$\begin{array}[]{lcl}\|w^{k+1}\|^{2}+(\omega-1)\|w^{k+1}-\hat{w}^{k+1}\|^{2}&\leq&\|w^{k}\|^{2}+2\omega
L^{2}\eta^{2}\|w^{k}-\hat{w}^{k}\|^{2}-(1-2\omega
L^{2}\eta^{2})\|\hat{w}^{k+1}-w^{k}\|^{2}.\end{array}$

Finally, let us choose $\omega\geq 1$ such that $\omega-1=2\omega
L^{2}\eta^{2}$. If $2L^{2}\eta^{2}<1$, then
$\omega:=\frac{1}{1-2L^{2}\eta^{2}}$ satisfies $\omega-1=2\omega
L^{2}\eta^{2}$. Consequently, the last estimate leads to (58). ∎

Now, we are ready to establish the convergence of (RFBS2) in the following
theorem.

######  Theorem 6.1.

Assume that $\mathrm{zer}(\Phi)\neq\emptyset$, $T$ in (NI) is maximally
monotone, and $F$ in (NI) is $L$-Lipschitz continuous and satisfies $\langle
Fx-Fx^{\star},x-x^{\star}\rangle\geq 0$ for all $x\in\mathrm{dom}(\Phi)$ and
some $x^{\star}\in\mathrm{zer}(\Phi)$. Let $\left\\{(x^{k},y^{k})\right\\}$ be
generated by (RFBS2) using $\eta\in\left(0,\frac{\sqrt{2}-1}{L}\right)$. Then,
we have the following statements.

  * •

The $\mathcal{O}\big{(}1/\sqrt{k}\big{)}$ best-iterate convergence rate. The
following bound holds:

$\frac{1}{k+1}\sum_{l=0}^{k}\|Fx^{l}+\xi^{l}\|^{2}\leq\frac{1}{k+1}\sum_{l=0}^{k}\left[\|Fx^{l}+\xi^{l}\|^{2}+\kappa\|Fx^{l}-Fy^{l-1}\|^{2}\right]\leq\frac{C_{0}\|x^{0}-x^{\star}\|^{2}}{k+1},$
(60)

where $\kappa:=\frac{2L^{2}\eta^{2}}{1-L^{2}\eta^{2}}>0$ and
$C_{0}:=\frac{5L^{2}\eta^{2}+3}{3\eta^{2}\left[1-(1+\sqrt{2})L\eta\right]}>0$.

  * •

The $\mathcal{O}\big{(}1/\sqrt{k}\big{)}$ last-iterate convergence rate. If
$\Phi$ is additionally monotone, then we also have

$\|Fx^{k}+\xi^{k}\|^{2}\leq\|Fx^{k}+\xi^{k}\|^{2}+\kappa\|Fx^{k}-Fy^{k-1}\|^{2}\leq\frac{C_{0}\|x^{0}-x^{\star}\|^{2}}{k+1}.$
(61)

As a consequence, we have the last-iterate convergence rate
$\mathcal{O}\left(\frac{1}{\sqrt{k}}\right)$ of the residual norm
$\|Fx^{k}+\xi^{k}\|$.

###### Proof.

First, since $0<\eta<\frac{\sqrt{2}-1}{L}$, we have $1-(1+\sqrt{2})L\eta>0$
and $\kappa:=\frac{2L^{2}\eta^{2}}{1-2L^{2}\eta^{2}}<\frac{2}{3}$. From (58),
we have

$\begin{array}[]{lcl}\|Fx^{k}+\xi^{k}\|^{2}+\kappa\|Fx^{k}-Fy^{k-1}\|^{2}&\leq&\frac{5L^{2}\eta^{2}+3}{3\eta^{2}}\|x^{k}-y^{k}\|^{2}+\left(\frac{5L^{2}\eta^{2}+3}{5\eta^{2}}+\frac{2L^{2}}{3}\right)\|x^{k}-y^{k-1}\|^{2}\vspace{1ex}\\\
&\leq&\frac{5L^{2}\eta^{2}+3}{3\eta^{2}}\left[\|x^{k}-y^{k}\|^{2}+\|x^{k}-y^{k-1}\|^{2}\right].\end{array}$

Combining this estimate and (54), we get

$\begin{array}[]{lcl}\frac{3\eta^{2}\left[1-(1+\sqrt{2})L\eta\right]}{5L^{2}\eta^{2}+3}\left[\|Fx^{k}+\xi^{k}\|^{2}+\kappa\|Fx^{k}-Fy^{k-1}\|^{2}\right]&\leq&(1-(1+\sqrt{2})L\eta)\left[\|x^{k}-y^{k}\|^{2}+\|x^{k}-y^{k-1}\|^{2}\right]\vspace{1ex}\\\
&\leq&\mathcal{V}_{k}-\mathcal{V}_{k+1}.\end{array}$

Summing up this inequality from $l=0$ to $l=k$, and using $y^{-1}:=x^{0}$, we
get

$\begin{array}[]{lcl}\frac{3\eta^{2}\left[1-(1+\sqrt{2})L\eta\right]}{5L^{2}\eta^{2}+3}\sum_{l=0}^{k}\left[\|Fx^{l}+\xi^{l}\|^{2}+\kappa\|Fx^{l}-Fy^{l-1}\|^{2}\right]&\leq&\mathcal{V}_{0}-\mathcal{V}_{k+1}\leq\mathcal{V}_{0}=\|x^{0}-x^{\star}\|^{2}.\end{array}$

This inequality implies (60). Finally, combining (60) and (59), we obtain
(61). ∎

######  Remark 6.1.

We can modify (RFBS2) to capture adaptive parameters as
$y^{k}:=x^{k}+\beta_{k}(x^{k}-x^{k-1})$ and
$x^{k+1}:=J_{\eta_{k}T}(x^{k}-\eta_{k}Fy^{k})$, where
$\eta_{k}:=\beta_{k}\eta_{k-1}$ for some $\beta_{k}>0$. Then, by imposing
appropriate bounds on $\eta_{k}$, we can still prove the convergence of this
variant by modifying the proof of Theorem 6.1. We also note that our best-
iterate convergence analysis of (RFBS2) in this paper is relatively different
from [30], while the last-iterate convergence rate analysis is new and very
simple.

###  6.2 The golden ratio method for (NI)

The golden ratio (GR) method for solving (NI) is presented as follows.
Starting from $x^{0}\in\mathrm{dom}(\Phi)$, we set $y^{-1}:=x^{0}$, and at
each iteration $k\geq 0$, we update

$\left\\{\begin{array}[]{lcl}y^{k}&:=&\frac{\omega-1}{\omega}x^{k}+\frac{1}{\omega}y^{k-1},\vspace{1ex}\\\
x^{k+1}&:=&J_{\eta T}(y^{k}-\eta Fx^{k}),\end{array}\right.$ (GR2)

where $J_{\eta T}$ is the resolvent of $\eta T$, $\omega>1$ is given, and
$\eta\in(0,\frac{\omega}{2L})$.

This method was proposed by Malitsky in [88] to solve monotone (MVIP), where
$\omega$ is chosen as $\omega:=\frac{\sqrt{5}+1}{2}$, leading to the name:
golden ratio. We now extend it to solve (NI) for the case $F$ is monotone and
$L$-Lipschitz continuous, and $T$ is maximally $3$-cyclically monotone.
Moreover, we extend our analysis for any $\omega\in(1,1+\sqrt{3})$ instead of
fixing $\omega:=\frac{\sqrt{5}+1}{2}$. We call this extension the GR2$+$
scheme.

Let us denote $\breve{w}^{k}:=Fx^{k-1}+\xi^{k}$ for $\xi^{k}\in Tx^{k}$. Then,
we can rewrite the second line of (GR2) as
$x^{k+1}:=y^{k}-\eta(Fx^{k}+\xi^{k+1})=y^{k}-\eta\breve{w}^{k+1}$ for
$\xi^{k+1}\in Tx^{k+1}$. In this case, we have
$x^{k}=y^{k-1}-\eta(Fx^{k-1}+\xi^{k})$, leading to
$y^{k-1}=x^{k}+\eta\breve{w}^{k}$. Combining this expression and the first
line of (GR2), we have
$y^{k}=\frac{\omega-1}{\omega}x^{k}+\frac{1}{\omega}(x^{k}+\eta\breve{w}^{k})=x^{k}+\frac{\eta}{\omega}\breve{w}^{k}$.
Consequently, we can rewrite (GR2) equivalently as follows:

$\left\\{\begin{array}[]{lcllcl}y^{k}&:=&x^{k}+\frac{\eta}{\omega}\breve{w}^{k},\vspace{1ex}\\\
x^{k+1}&:=&y^{k}-\eta(Fx^{k}+\xi^{k+1})&=&y^{k}-\eta\breve{w}^{k+1}.\end{array}\right.$
(62)

If we eliminate $y^{k}$, then we obtain

$\left\\{\begin{array}[]{lcl}x^{k+1}&:=&J_{\eta
T}\left(x^{k}-\eta\left(Fx^{k}-\frac{1}{\omega}(Fx^{k-1}+\xi^{k})\right)\right)=x^{k}-\eta\breve{w}^{k+1}+\frac{\eta}{\omega}\breve{w}^{k},\vspace{1ex}\\\
\xi^{k+1}&:=&\frac{1}{\eta}(x^{k}-x^{k+1})-\left(Fx^{k}-\frac{1}{\omega}(Fx^{k-1}+\xi^{k})\right).\end{array}\right.$
(63)

The convergence of (GR2) is established based on the following key lemma.

######  Lemma 6.3.

Suppose that $\mathrm{zer}(\Phi)\neq\emptyset$, $T$ in (NI) is maximally
$3$-cyclically monotone, and $F$ is $L$-Lipschitz continuous. Let
$\left\\{(x^{k},y^{k})\right\\}$ be generated by (GR2) with $\omega>1$. Then,
for any $x^{\star}\in\mathrm{zer}(\Phi)$, we have

$\begin{array}[]{lcl}\omega\|y^{k+1}-x^{\star}\|^{2}+(\omega-1)(\omega-\gamma)\|x^{k+1}-x^{k}\|^{2}&\leq&\omega\|y^{k}-x^{\star}\|^{2}+\frac{(\omega-1)L^{2}\eta^{2}}{\gamma}\|x^{k}-x^{k-1}\|^{2}\vspace{1ex}\\\
&&-{~{}}\omega(\omega-1)\|x^{k}-y^{k}\|^{2}-\frac{(\omega-1)(1-\omega^{2}+\omega)}{\omega}\|x^{k+1}-y^{k}\|^{2}\vspace{1ex}\\\
&&-{~{}}2\eta(\omega-1)\langle
Fx^{k}-Fx^{\star},x^{k}-x^{\star}\rangle.\end{array}$ (64)

###### Proof.

Since $T$ is $3$-cyclically monotone, for $\xi^{k+1}\in Tx^{k+1}$, $\xi^{k}\in
Tx^{k}$, and $x^{\star}\in Tx^{\star}$, we have

$\begin{array}[]{lcll}\langle\xi^{k+1},x^{k+1}-x^{\star}\rangle+\langle\xi^{\star},x^{\star}-x^{k}\rangle+\langle\xi^{k},x^{k}-x^{k+1}\rangle\geq
0.\end{array}$ (65)

From (62), we have $\eta\xi^{k+1}=y^{k}-x^{k+1}-\eta Fx^{k}$ and
$\eta\xi^{k}=y^{k-1}-x^{k}-\eta Fx^{k-1}$. Moreover, since
$x^{\star}\in\mathrm{zer}(\Phi)$, we have $Fx^{\star}+\xi^{\star}=0$, leading
to $\eta\xi^{\star}=-\eta Fx^{\star}$. Substituting these expressions into
(65), we have

$\begin{array}[]{lcll}\langle y^{k}-x^{k+1}-\eta
Fx^{k},x^{k+1}-x^{\star}\rangle+\langle y^{k-1}-x^{k}-\eta
Fx^{k-1},x^{k}-x^{k+1}\rangle-\langle Fx^{\star},x^{\star}-x^{k}\rangle\geq
0.\end{array}$

However, since $y^{k}:=\frac{\omega-1}{\omega}x^{k}+\frac{1}{\omega}y^{k-1}$
from the first line of (GR2), we get $y^{k-1}-x^{k}=\omega(y^{k}-x^{k})$.
Substituting this relation into the last inequality, rearranging the result,
we obtain

$\langle y^{k}-x^{k+1},x^{k+1}-x^{\star}\rangle+\omega\langle
x^{k}-y^{k},x^{k+1}-x^{k}\rangle+\eta\langle
Fx^{k-1}-Fx^{k},x^{k+1}-x^{k}\rangle-\eta\langle
Fx^{k}-Fx^{\star},x^{k}-x^{\star}\rangle\geq 0.$ (66)

By Young’s inequality and the $L$-Lipschitz continuity of $F$, for any
$\gamma>0$, we get the first line of the following:

$\begin{array}[]{lcl}2\eta\langle
Fx^{k-1}-Fx^{k},x^{k+1}-x^{k}\rangle&\leq&\frac{L^{2}\eta^{2}}{\gamma}\|x^{k}-x^{k-1}\|^{2}+\gamma\|x^{k+1}-x^{k}\|^{2},\vspace{1ex}\\\
2\langle
y^{k}-x^{k+1},x^{k+1}-x^{\star}\rangle&=&\|y^{k}-x^{\star}\|^{2}-\|x^{k+1}-x^{\star}\|^{2}-\|x^{k+1}-y^{k}\|^{2},\vspace{1ex}\\\
2\langle
x^{k}-y^{k},x^{k+1}-x^{k}\rangle&=&\|x^{k+1}-y^{k}\|^{2}-\|x^{k}-y^{k}\|^{2}-\|x^{k+1}-x^{k}\|^{2}.\end{array}$

Substituting these expressions into (66), and rearranging the result, we can
show that

$\begin{array}[]{lcl}\|x^{k+1}-x^{\star}\|^{2}&\leq&\|y^{k}-x^{\star}\|^{2}+(\omega-1)\|x^{k+1}-y^{k}\|^{2}-\omega\|x^{k}-y^{k}\|^{2}-(\omega-\gamma)\|x^{k+1}-x^{k}\|^{2}\vspace{1ex}\\\
&&+{~{}}\frac{L^{2}\eta^{2}}{\gamma}\|x^{k}-x^{k-1}\|^{2}-2\eta\langle
Fx^{k}-Fx^{\star},x^{k}-x^{\star}\rangle.\end{array}$ (67)

Now, using $(\omega-1)x^{k+1}=\omega y^{k+1}-y^{k}$ and
$\omega(y^{k+1}-y^{k})=(\omega-1)(x^{k+1}-y^{k})$ from the first line of
(GR2), we can derive that

$\begin{array}[]{lcl}(\omega-1)^{2}\|x^{k+1}-x^{\star}\|^{2}&=&\omega(\omega-1)\|y^{k+1}-x^{\star}\|^{2}-(\omega-1)\|y^{k}-x^{\star}\|^{2}+\omega\|y^{k+1}-y^{k}\|^{2}\vspace{1ex}\\\
&=&\omega(\omega-1)\|y^{k+1}-x^{\star}\|^{2}-(\omega-1)\|y^{k}-x^{\star}\|^{2}+\frac{(\omega-1)^{2}}{\omega}\|x^{k+1}-y^{k}\|^{2}.\end{array}$

Simplifying this expression to get
$(\omega-1)\|x^{k+1}-x^{\star}\|^{2}=\omega\|y^{k+1}-x^{\star}\|^{2}-\|y^{k}-x^{\star}\|^{2}+\frac{(\omega-1)}{\omega}\|x^{k+1}-y^{k}\|^{2}$.
Combining it and (67), and rearranging the result, we obtain (64). ∎

Now, we are ready to state the convergence of (GR2) in the following theorem.

######  Theorem 6.2.

Assume that $\mathrm{zer}(\Phi)\neq\emptyset$, $T$ in (NI) is maximally
$3$-cyclically monotone, and $F$ in (NI) is $L$-Lipschitz continuous and
satisfies $\langle Fx-Fx^{\star},x-x^{\star}\rangle\geq 0$ for all
$x\in\mathrm{dom}(\Phi)$ and some $x^{\star}\in\mathrm{zer}(\Phi)$. Let
$\left\\{(x^{k},y^{k})\right\\}$ be generated by (GR2). Then, the following
statements hold.

  * •

The best-iterate rate of GR2. If $1<\omega\leq\frac{1+\sqrt{5}}{2}$ and
$\eta\in\left(0,\frac{\omega}{2L}\right)$, then

$\frac{1}{k+1}\sum_{l=0}^{k}\|Fx^{l}+\xi^{l}\|^{2}\leq\frac{1}{k+1}\sum_{l=0}^{k}(\omega-1)\left[\omega\|x^{l}-y^{l}\|^{2}+\varphi\cdot\|x^{l}-x^{l-1}\|^{2}\right]\leq\frac{C_{0}\|x^{0}-x^{\star}\|^{2}}{k+1},$
(68)

where $\varphi:=\frac{\omega^{2}-4L^{2}\eta^{2}}{2\omega}>0$ and
$C_{0}:=\frac{(\omega^{2}-2L^{2}\eta^{2})\omega^{2}}{(\omega^{2}-4L^{2}\eta^{2})\eta^{2}(\omega-1)}>0$.

  * •

The best-iterate rate for GR2$+$. If $\frac{1+\sqrt{5}}{2}<\omega<1+\sqrt{3}$
and $0<\eta<\frac{\psi}{2L}$, then

$\frac{1}{k+1}\sum_{l=0}^{k}\|Fx^{l}+\xi^{l}\|^{2}\leq\frac{1}{k+1}\sum_{l=0}^{k}(\omega-1)\left[\psi\cdot\|x^{l}-y^{l}\|^{2}+\kappa\cdot\|x^{l}-y^{l-1}\|^{2}\right]\leq\frac{\hat{C}_{0}\|x^{0}-x^{\star}\|^{2}}{k+1},$
(69)

where $\psi:=\frac{2\omega+2-\omega^{2}}{\omega}>0$,
$\kappa:=\frac{\psi^{2}-4L^{2}\eta^{2}}{2\psi}$, and
$\hat{C}_{0}:=\frac{[\psi^{2}-2L^{2}\eta^{2}(2\omega^{2}-\psi^{2})]\omega}{(\omega-1)(\psi^{2}-4L^{2}\eta^{2})\eta^{2}\psi}>0$.

###### Proof.

First, to guarantee that $1+\omega-\omega^{2}\geq 0$ and $\omega>1$, we need
to choose $1<\omega\leq\frac{\sqrt{5}+1}{2}$. If $0<\eta<\frac{\omega}{2L}$,
then by choosing $\gamma:=\frac{\omega}{2}$, we have
$\psi:=\frac{(\omega-1)(\omega\gamma-\gamma^{2}-L^{2}\eta^{2})}{\gamma}=\frac{(\omega-1)(\omega^{2}-4L^{2}\eta^{2})}{2\omega}>0$.
Using this relation and $\langle Fx^{k}-Fx^{\star},x^{k}-x^{\star}\rangle\geq
0$, if we define
$\mathcal{V}_{k}:=\omega\|y^{k}-x^{\star}\|^{2}+\frac{\omega(\omega-1)}{2}\|x^{k}-x^{k-1}\|^{2}\geq
0$, then we can deduce from (64) that

$\begin{array}[]{lcl}\mathcal{V}_{k+1}&\leq&\mathcal{V}_{k}-\psi\cdot\|x^{k}-x^{k-1}\|^{2}-\omega(\omega-1)\|x^{k}-y^{k}\|^{2}.\end{array}$
(70)

Next, using $y^{k}-x^{k}=\frac{\eta}{\omega}\breve{w}^{k}$ and
$\breve{w}^{k}=Fx^{k-1}+\xi^{k}$, by Young’s inequality, we have

$\begin{array}[]{lcl}\|w^{k}\|^{2}&=&\|Fx^{k}+\xi^{k}\|^{2}\leq\left(1+\frac{\psi\omega}{L^{2}\eta^{2}(\omega-1)}\right)\|Fx^{k}-Fx^{k-1}\|^{2}+\left(1+\frac{L^{2}\eta^{2}(\omega-1)}{\psi\omega}\right)\|\breve{w}^{k}\|^{2}\vspace{1ex}\\\
&\leq&\left(1+\frac{\psi\omega}{L^{2}\eta^{2}(\omega-1)}\right)L^{2}\|x^{k}-x^{k-1}\|^{2}+\left(1+\frac{L^{2}\eta^{2}(\omega-1)}{\psi\omega}\right)\frac{\omega^{2}}{\eta^{2}}\|x^{k}-y^{k}\|^{2}\vspace{1ex}\\\
&=&\frac{(\omega-1)}{\psi\eta^{2}(\omega-1)}\left[\psi\cdot\|x^{k}-x^{k-1}\|^{2}+\omega(\omega-1)\|x^{k}-y^{k}\|^{2}\right].\end{array}$
(71)

Combining this estimate and (70), and noting that $\mathcal{V}_{k}\geq 0$, we
can show that

$\begin{array}[]{lcl}\sum_{l=0}^{k}\|w^{l}\|^{2}&\leq&\frac{L^{2}\eta^{2}(\omega-1)+\psi\omega}{\psi\eta^{2}(\omega-1)}\sum_{l=0}^{k}\left[\psi\cdot\|x^{l}-x^{l-1}\|^{2}+\omega(\omega-1)\|x^{l}-y^{l}\|^{2}\right]\vspace{1ex}\\\
&\leq&\frac{L^{2}\eta^{2}(\omega-1)+\psi\omega}{\psi\eta^{2}(\omega-1)}\left[\mathcal{V}_{0}-\mathcal{V}_{k+1}\right]\leq\frac{L^{2}\eta^{2}(\omega-1)+\psi\omega}{\psi\eta^{2}(\omega-1)}\cdot\mathcal{V}_{0}\vspace{1ex}\\\
&=&\frac{(\omega^{2}-2L^{2}\eta^{2})\omega^{2}}{(\omega^{2}-4L^{2}\eta^{2})\eta^{2}(\omega-1)}\cdot\|x^{0}-x^{\star}\|^{2},\end{array}$

which is exactly (68), where we have used
$\mathcal{V}_{0}:=\omega\|y^{0}-x^{\star}\|^{2}+\frac{\omega(\omega-1)}{2}\|x^{0}-x^{-1}\|^{2}=\omega\|x^{0}-x^{\star}\|^{2}$
due to $x^{-1}=y^{0}=x^{0}$.

Next, if $1.6180\approx\frac{1+\sqrt{5}}{2}<\omega<1+\sqrt{3}\approx 2.7321$,
then we have $\omega^{2}-\omega-1>0$ and
$\psi:=\omega-\frac{2(\omega^{2}-\omega-1)}{\omega}>0$. In this case, using
$\|x^{k+1}-y^{k}\|^{2}\leq 2\|x^{k+1}-x^{k}\|^{2}+2\|y^{k}-x^{k}\|^{2}$ and
$\langle Fx^{k}-Fx^{\star},x^{k}-x^{\star}\rangle\geq 0$ into (64),
rearranging the result, and using $\gamma:=\frac{\psi}{2}$, we get

$\begin{array}[]{lcl}\omega\|y^{k+1}-x^{\star}\|^{2}+\frac{\psi(\omega-1)}{2}\|x^{k+1}-x^{k}\|^{2}&\leq&\omega\|y^{k}-x^{\star}\|^{2}+\frac{\psi(\omega-1)}{2}\|x^{k}-x^{k-1}\|^{2}-\psi(\omega-1)\|x^{k}-y^{k}\|^{2}\vspace{1ex}\\\
&&-{~{}}\frac{(\omega-1)(\psi^{2}-4L^{2}\eta^{2})}{2\psi}\|x^{k}-x^{k-1}\|^{2}.\end{array}$
(72)

Similar to the proof of (71), we have
$\|w^{k}\|^{2}\leq\frac{\psi^{2}-2L^{2}\eta^{2}(2\omega^{2}-\psi^{2})}{(\psi^{2}-4L^{2}\eta^{2})\eta^{2}\psi}\big{[}\frac{\psi^{2}-4L^{2}\eta^{2}}{2\psi}\|x^{k}-x^{k-1}\|^{2}+\psi\|x^{k}-y^{k}\|^{2}\big{]}$.
Combining this inequality and (72), with same argument as in the proof of
(68), we obtain (69). ∎

##  7 Accelerated Extragradient Methods for Nonlinear Inclusions

Introduction. The convergence rate on the residual norm $\|Fx^{k}+\xi^{k}\|$
of the EG method and its variants using constant stepsize discussed so far is
$\mathcal{O}\big{(}1/\sqrt{k}\big{)}$, which is unimprovable for standard and
constant stepsize EG-type methods as shown in [54]. In this section, we survey
recent development on accelerated methods that can theoretically achieve a
$\mathcal{O}\left(1/k\right)$ last-iterate convergence rate on
$\|Fx^{k}+\xi^{k}\|$ using variable stepsizes. We will present two different
approaches to develop accelerated methods for solving (NE) and (NI) without
using averaging sequences. The first one relies on Halpern’s fixed-point
iteration [60], and the second approach leverages Nesterov’s accelerated
techniques. Halpern’s fixed-point iteration is a classical method to
approximate a fixed-point of a nonexpansive operator, or equivalently, to find
a root of a co-coercive operator. This method has been intensively studied in
fixed-point theory, but the first work showing a $\mathcal{O}\left(1/k\right)$
last-iterate convergence rate is due to F. Lieder in [77]. This method was
then extended and intensively studied in [42] for root-finding problems and
VIPs. In a pioneering work [143], Yoon and Ryu extended Halpern’s fixed-point
iteration to the EG method for solving (NE), which is called extra-anchored
gradient (EAG) method. This new method still achieves
$\mathcal{O}\left(1/k\right)$ last-iterate convergence but only requires the
monotonicity and Lipschitz continuity of $F$. Lee and Kim in [75] further
advanced [143] to the co-hypomonotone setting of (NE) and still achieved the
same rates. The authors in [132] exploited the technique in [143] and applied
it to the past-extragradient method in [113] and obtained a past extra-
anchored gradient (PEAG) method that has the same
$\mathcal{O}\left(1/k\right)$-rates (up to a constant factor). Recently, [25]
and [27] expanded the results in [75, 132, 143] to develop methods for solving
(VIP) and (NI), and preserved the same $\mathcal{O}\left(1/k\right)$ last-
iterate convergence rates on $\|Fx^{k}+\xi^{k}\|$.

In this section, we summarize the above mentioned results and provide a
unified convergence analysis obtained from a recent work [130] that covers all
the results from [25, 27, 75, 132, 143] in a unified fashion.

###  7.1 The extra-anchored gradient method for (NI)

The algorithm. The extra-anchored gradient (EAG) for solving (NI) we discuss
here is presented as follows. Starting from $x^{0}\in\mathrm{dom}(\Phi)$, at
each iteration $k\geq 0$, we update

$\left\\{\begin{array}[]{lcl}y^{k}&:=&J_{\hat{\eta}_{k}T}\left(\tau_{k}x^{0}+(1-\tau_{k})x^{k}-\hat{\eta}_{k}Fx^{k}\right),\vspace{1ex}\\\
x^{k+1}&:=&J_{\eta_{k}T}\left(\tau_{k}x^{0}+(1-\tau_{k})x^{k}-\eta_{k}Fy^{k}\right),\end{array}\right.$
(EAG2)

where $\tau_{k}\in(0,1)$, $\hat{\eta}_{k}>0$, and $\eta_{k}>0$ are given,
which will be determined later. Here, we assume that $T$ is maximally
$3$-cyclically monotone, and hence covers the special cases
$T=\mathcal{N}_{\mathcal{X}}$, the normal cone of $\mathcal{X}$ and
$T=\partial{g}$, the subdifferential of a convex function $g$. In fact, [25]
considers the special case (VIP) of (NI) when $T:=\mathcal{N}_{\mathcal{X}}$
is normal cone of a nonempty, closed, and convex set $\mathcal{X}$, and
$\hat{\eta}_{k}:=\eta_{k}$. As we can see, the scheme (EAG2) purely extends
the extra-anchored gradient (EAG) scheme from [143] for (NE) to (NI), when $F$
is monotone and Lipschitz continuous, and $T$ is maximally $3$-cyclically
monotone.

Convergence analysis. For simplicity of analysis, we recall the following
quantities defined ealier:

$w^{k}:=Fx^{k}+\xi^{k},\quad\hat{w}^{k}:=Fy^{k-1}+\xi^{k},\quad\text{and}\quad\tilde{w}^{k}:=Fx^{k}+\zeta^{k},$
(73)

where $\xi^{k}\in Tx^{k}$ and $\zeta^{k}\in Ty^{k}$. Then, we can equivalently
rewrite (EAG2) as follows:

$\left\\{\begin{array}[]{lcl}y^{k}&:=&\tau_{k}x^{0}+(1-\tau_{k})x^{k}-\hat{\eta}_{k}\tilde{w}^{k},\vspace{1ex}\\\
x^{k+1}&:=&\tau_{k}x^{0}+(1-\tau_{k})x^{k}-\eta_{k}\hat{w}^{k+1}.\end{array}\right.$
(74)

To establish the convergence of (EAG2), we use the following potential
function as in [25, 75, 130, 143]:

$\mathcal{V}_{k}:=a_{k}\|w^{k}\|^{2}+b_{k}\langle w^{k},x^{k}-x^{0}\rangle,$
(75)

where $a_{k}>0$ and $b_{k}>0$ are given parameters. Let us prove the
convergence of (EAG2).

######  Theorem 7.1.

For (NI), assume that $\mathrm{zer}(\Phi)\neq\emptyset$, $F$ is $L$-Lipschitz
continuous and monotone, and $T$ is maximally $3$-cyclically monotone. Let
$\\{(x^{k},y^{k})\\}$ be generated by (EAG2) using

$\tau_{k}:=\frac{1}{k+2},\quad\eta_{k}:=\eta\in\left(0,\frac{1}{L}\right],\quad\text{and}\quad\hat{\eta}_{k}:=(1-\tau_{k})\eta.$
(76)

Then, for all $k\geq 0$ and any $x^{\star}\in\mathrm{zer}(\Phi)$, the
following result holds:

$\|Fx^{k}+\xi^{k}\|^{2}\leq\frac{4\|x^{0}-x^{\star}\|^{2}+2\eta^{2}\|Fx^{0}+\xi^{0}\|^{2}}{\eta^{2}(k+1)^{2}},\quad\text{where}\quad\xi^{k}\in
Tx^{k}.$ (77)

Consequently, we have the last-iterate convergence rate
$\mathcal{O}\left(1/k\right)$ of the residual norm $\|Fx^{k}+\xi^{k}\|$.

###### Proof.

Since $T$ is maximally $3$-cyclically monotone, $\xi^{k+1}\in Tx^{k+1}$,
$\xi^{k}\in Tx^{k}$, and $\zeta^{k}\in Ty^{k}$, we have

$\begin{array}[]{lcl}\langle\xi^{k+1},x^{k+1}-x^{k}\rangle+\langle\xi^{k},x^{k}-y^{k}\rangle+\langle\zeta^{k},y^{k}-x^{k+1}\rangle\geq
0.\end{array}$

By the monotonicity of $F$, we also have $\langle
Fx^{k+1}-Fx^{k},x^{k+1}-x^{k}\rangle\geq 0$. Summing up this inequality and
the last one, and then using the fact that $w^{k+1}=Fx^{k+1}+\xi^{k+1}$,
$w^{k}:=Fx^{k}+\xi^{k}$, and $\tilde{w}^{k}:=Fx^{k}+\zeta^{k}$, we obtain

$\begin{array}[]{lcl}\langle
w^{k+1},x^{k+1}-x^{k}\rangle-\langle\tilde{w}^{k},x^{k+1}-x^{k}\rangle+\langle
w^{k}-\tilde{w}^{k},x^{k}-y^{k}\rangle\geq 0.\end{array}$ (78)

Now, from (74), we have

$\begin{array}[]{lcl}x^{k+1}-x^{k}&=&-\frac{\tau_{k}}{1-\tau_{k}}(x^{k+1}-x^{0})-\frac{\eta_{k}}{1-\tau_{k}}\hat{w}^{k+1},\vspace{1ex}\\\
x^{k+1}-x^{k}&=&-\tau_{k}(x^{k}-x^{0})-\eta_{k}\hat{w}^{k+1},\vspace{1ex}\\\
x^{k}-y^{k}&=&\tau_{k}(x^{k}-x^{0})+\hat{\eta}_{k}\tilde{w}^{k}.\end{array}$

Substituting these relations into (78), and rearranging terms, we arrive at

$\begin{array}[]{lcl}\tau_{k}\langle
w^{k},x^{k}-x^{0}\rangle-\frac{\tau_{k}}{1-\tau_{k}}\langle
w^{k+1},x^{k+1}-x^{0}\rangle&\geq&\frac{\eta_{k}}{1-\tau_{k}}\langle
w^{k+1},\hat{w}^{k+1}\rangle-\eta_{k}\langle\tilde{w}^{k},\hat{w}^{k+1}\rangle-\hat{\eta}_{k}\langle
w^{k},\tilde{w}^{k}\rangle+\hat{\eta}_{k}\|\tilde{w}^{k}\|^{2}.\end{array}$

Multiplying this inequality by $\frac{b_{k}}{\tau_{k}}$ and assume that
$b_{k+1}=\frac{b_{k}}{1-\tau_{k}}$, and then using (75), we can show that

$\begin{array}[]{lcl}\mathcal{V}_{k}-\mathcal{V}_{k+1}&=&a_{k}\|w^{k}\|^{2}-a_{k+1}\|w^{k+1}\|^{2}+b_{k}\langle
w^{k},x^{k}-x^{0}\rangle-\frac{b_{k}}{1-\tau_{k}}\langle
w^{k+1},x^{k+1}-x^{0}\rangle\vspace{1ex}\\\
&\geq&\frac{b_{k+1}\eta_{k}}{\tau_{k}}\langle
w^{k+1}-\tilde{w}^{k},\hat{w}^{k+1}\rangle+b_{k+1}\eta_{k}\langle\tilde{w}^{k},\hat{w}^{k+1}\rangle-\frac{b_{k}\hat{\eta}_{k}}{\tau_{k}}\langle
w^{k},\tilde{w}^{k}\rangle+\frac{b_{k}\hat{\eta}_{k}}{\tau_{k}}\|\tilde{w}^{k}\|^{2}\vspace{1ex}\\\
&&+{~{}}a_{k}\|w^{k}\|^{2}-a_{k+1}\|w^{k+1}\|^{2}.\end{array}$ (79)

Next, from (74), we have
$x^{k+1}-y^{k}=-\eta_{k}\hat{w}^{k+1}+\hat{\eta}_{k}\tilde{w}^{k}$. Using this
expression and the $L$-Lipschitz continuity of $F$, we have
$\|w^{k+1}-\hat{w}^{k+1}\|^{2}=\|Fx^{k+1}-Fy^{k}\|^{2}\leq
L^{2}\|x^{k+1}-y^{k}\|^{2}=L^{2}\|\eta_{k}\hat{w}^{k+1}-\hat{\eta}_{k}\tilde{w}^{k}\|^{2}$,
leading to

$\begin{array}[]{lcl}\|w^{k+1}\|^{2}+(1-L^{2}\eta_{k}^{2})\|\hat{w}^{k+1}\|^{2}-2\langle
w^{k+1}-\tilde{w}^{k},\hat{w}^{k+1}\rangle-2(1-L^{2}\eta_{k}\hat{\eta}_{k})\langle\hat{w}^{k+1},\tilde{w}^{k}\rangle-L^{2}\hat{\eta}_{k}^{2}\|\tilde{w}^{k}\|^{2}\leq
0.\end{array}$

Multiplying this inequality by $\frac{b_{k+1}\eta_{k}}{2\tau_{k}}$, adding the
result to (78), and using $\hat{\eta}_{k}=(1-\tau_{k})\eta_{k}$, we obtain

$\begin{array}[]{lcl}\mathcal{V}_{k}-\mathcal{V}_{k+1}&\geq&\left(\frac{b_{k+1}\eta_{k}}{2\tau_{k}}-a_{k+1}\right)\|w^{k+1}\|^{2}+\frac{b_{k+1}\eta_{k}(1-L^{2}\eta_{k}^{2})}{2\tau_{k}}\|\hat{w}^{k+1}\|^{2}+a_{k}\|w^{k}\|^{2}+\frac{b_{k+1}\eta_{k}(1-\tau_{k})^{2}(2-L^{2}\eta_{k}^{2})}{2\tau_{k}}\|\tilde{w}^{k}\|^{2}\vspace{1ex}\\\
&&-{~{}}\frac{b_{k+1}\eta_{k}(1-L^{2}\eta_{k}^{2})(1-\tau_{k})}{\tau_{k}}\langle\tilde{w}^{k},\hat{w}^{k+1}\rangle-\frac{b_{k+1}\eta_{k}(1-\tau_{k})^{2}}{\tau_{k}}\langle
w^{k},\tilde{w}^{k}\rangle\vspace{1ex}\\\
&=&\frac{b_{k+1}\eta_{k}(1-L^{2}\eta_{k}^{2})}{2\tau_{k}}\|\hat{w}^{k+1}-(1-\tau_{k})\tilde{w}^{k}\|^{2}+\frac{b_{k+1}\eta_{k}(1-\tau_{k})^{2}}{2\tau_{k}}\|w^{k}-\tilde{w}^{k}\|^{2}\vspace{1ex}\\\
&&+{~{}}\left(\frac{b_{k+1}\eta_{k}}{2\tau_{k}}-a_{k+1}\right)\|w^{k+1}\|^{2}+\left(a_{k}-\frac{b_{k}\eta_{k}(1-\tau_{k})}{2\tau_{k}}\right)\|w^{k}\|^{2}.\end{array}$

Let us choose $\eta_{k}:=\eta\in\left(0,\frac{1}{L}\right]$ as in (76),
$\tau_{k}:=\frac{1}{k+2}$ and
$a_{k}:=\frac{b_{k}\eta(1-\tau_{k})}{2\tau_{k}}=\frac{\eta b_{k}(k+1)}{2}$.
Then, we have $a_{k+1}=\frac{\eta b_{k+1}(k+2)}{2}=\frac{\eta
b_{k+1}}{2\tau_{k}}$. Moreover, since
$b_{k+1}=\frac{b_{k}}{1-\tau_{k}}=\frac{b_{k}(k+2)}{k+1}$. By induction, we
obtain $b_{k}=b_{0}(k+1)$ for some $b_{0}>0$. Using these parameters into the
last estimate, we obtain

$\begin{array}[]{lcl}\mathcal{V}_{k}-\mathcal{V}_{k+1}&\geq&\frac{b_{k+1}\eta(1-L^{2}\eta^{2})}{2\tau_{k}}\|\hat{w}^{k+1}-(1-\tau_{k})\tilde{w}^{k}\|^{2}+\frac{b_{k+1}\eta(1-\tau_{k})^{2}}{2\tau_{k}}\|w^{k}-\tilde{w}^{k}\|^{2}\geq
0.\end{array}$

Finally, using $\langle w^{k},x^{k}-x^{\star}\rangle\geq 0$ for