# Fast Nonlinear Two-Time-Scale Stochastic Approximation: Achieving ${\cal
O}(1/k)$ Finite-Sample Complexity

Thinh T. Doan Thinh T. Doan is with the Bradley Department of Electrical and
Computer Engineering, Virginia Tech, USA. Email<EMAIL_ADDRESS>
###### Abstract

This paper proposes to develop a new variant of the two-time-scale stochastic
approximation to find the roots of two coupled nonlinear operators, assuming
only noisy samples of these operators can be observed. Our key idea is to
leverage the classic Ruppert–Polyak averaging technique to dynamically
estimate the operators through their samples. The estimated values of these
averaging steps will then be used in the two-time-scale stochastic
approximation updates to find the desired solution. Our main theoretical
result is to show that under the strongly monotone condition of the underlying
nonlinear operators the mean-squared errors of the iterates generated by the
proposed method converge to zero at an optimal rate ${\cal O}(1/k)$, where $k$
is the number of iterations. Our result significantly improves the existing
result of two-time-scale stochastic approximation, where the best known
finite-time convergence rate is ${\cal O}(1/k^{2/3})$.

##  1 Introduction

In this paper, we study the problem of finding the root of two coupled
nonlinear operators. In particular, given two operators
$F:\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ and
$G:\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$, we seek to
find $x^{\star}$ and $y^{\star}$ such that

$\displaystyle\left\\{\begin{aligned} &F(x^{\star},y^{\star})=0\\\
&G(x^{\star},y^{\star})=0.\end{aligned}\right.$ (1.1)

We will be interested in the setting where $F,G$ are unknown. Instead, we
assume that there is a stochastic oracle that outputs noisy values of $F(x,y)$
and $G(x,y)$ for a given pair $(x,y)$, i.e., for any $x$ and $y$ we have
access to $F(x,y)\,+\,\xi$ and $G(x,y)\,+\,\psi$, where $\xi$ and $\psi$ are
two random variables. For solving problem (1.1) in this stochastic setting,
two-time-scale stochastic approximation (SA) [1], a generalized variant of the
classic SA [2], is among the most popular methods and has been studied
extensively in the literature. This method iteratively updates the variables
$x_{k}$ and $y_{k}$, the estimates of $x^{\star}$ and $y^{\star}$,
respectively, for all $k\geq 0$ as

$\displaystyle\begin{aligned}
x_{k+1}&=x_{k}-\alpha_{k}\left(F(x_{k},y_{k})+\xi_{k}\right)\\\
y_{k+1}&=y_{k}-\beta_{k}\left(G(x_{k},y_{k})+\psi_{k}\right),\end{aligned}$
(1.2)

where $x_{0}$ and $y_{0}$ are arbitrarily initialized in $\mathbb{R}^{d}$. In
(1.2), $\alpha_{k}$ and $\beta_{k}$ are two nonnegative step sizes chosen such
that $\beta_{k}\ll\alpha_{k}$, i.e., the second iterate is updated using step
sizes that are very small as compared to the ones used to update the first
iterate. Thus, the update of $x_{k}$ is referred to as the “fast-time scale”
while the update of $y_{k}$ is called the “slow-time scale”, explaining the
name of two-time-scale SA. The time-scale difference here is loosely defined
as the ratio between the two step sizes, i.e., $\beta_{k}/\alpha_{k}$. In
addition, the update of the fast iterate depends on the slow iterate and vice
versa, that is, they are coupled to each other. To handle this coupling, the
two step sizes have to be properly chosen to guarantee the convergence of the
method. Indeed, an important problem in this area is to properly select the
two step sizes so that the two iterates converge as fast as possible. This
task, however, is challenging due to the coupling between the iterates and the
impact of sampling noise on their updates.

At present, the asymptotic convergence of two-time-scale SA, i.e., whether
$(x_{k},y_{k})$ converges to $(x^{\star},y^{\star})$, is well-studied by using
the method of ordinary differential equation (ODE). This result can be found
in several textbooks, see for example [1, 3]. The asymptotic convergence rates
of this method have been also studied via the central limit theorem in both
linear and nonlinear settings [4, 5, 6]. In particular, these works show that
two-time-scale SA asymptotically converges in distribution at a rate ${\cal
O}(1/k)$. On the other hand, the finite-sample complexity (or finite-time
convergence rate) of this method in the mean-squared regime is rather
incomplete. In the linear setting, i.e., when $F,G$ are linear on $x,y$, the
mean-square error of the iterates generated by (1.2) converges to zero at a
rate ${\cal O}(1/k)$ [7, 8]. When $F,G$ are nonlinear, only a suboptimal
convergence rate, ${\cal O}(1/k^{2/3})$, is known [9]. Motivated by the work
in [9], the authors in [10] show that the two-time-scale SA in (1.2) converges
at a rate ${\cal O}(1/k)$, however, this result requires a restrictive
condition on the high order of the smoothness of the operators (see the
discussion in Remark 4 for more detail). This leaves the question of whether
the two-time-scale SA in (1.2) can find the solution $(x^{\star},y^{\star})$
of (1.1) with an ${\cal O}(1/k)$ finite-time sample complexity. The
overarching goal of this paper is to (partially) answer this question, where
our main contribution is:

#### Main Contributions.

This paper proposes to develop a new variant of the two-time-scale SA in
(1.2), leveraging the classic Ruppert-Polyak averaging technique. The mean
square errors of the iterates generated by the proposed methods are shown to
converge to zero at a rate ${\cal O}(1/k)$, which significantly improves the
best known result ${\cal O}(1/k^{2/3})$ in the setting of strongly monotone,
nonlinear operators.

###  1.1 Prior Works

#### Stochastic Approximation (SA).

The two-time-scale SA in (1.2) is a generalized variant of the classic SA
introduced by [2] with broad applications in various areas including
statistics, stochastic optimization, machine learning, and reinforcement
learning [11, 1, 12, 13, 14]. In general, SA is a simulation-based approach
for finding the root (or fixed point) of some unknown operator $F$ represented
by the form of an expectation, i.e., $F(x)=\mathbb{E}_{\pi}[F(x,\xi)]$, where
$\xi$ is some random variable with a distribution $\pi$. Specifically, this
method seeks a point $x^{\star}$ such that $F(x^{\star})=0$ based on the noisy
observations $F(x;\xi)$. The estimate $x$ of $x^{\star}$ is iteratively
updated by moving along the direction of $F(x;\xi)$ scaled by some step size.
Through a careful choice of this step size, the “noise” induced by the random
samples $\xi$ can be averaged out across iterations, and the algorithm
converges to $x^{*}$. Both asymptotic convergence and finite-sample complexity
of SA are very well-studied under general conditions of the operators and
noise models; see for example [1, 3, 14, 15, 16, 17, 18, 19, 20, 21]. One can
apply SA to solve problem (1.1), i.e., by considering $W=[FG]$ as a single
operator. This approach requires that $W$ has to satisfy a similar condition
as used in the existing literature of SA. On the other hand, in this paper we
consider a more general condition of $F,G$ that are not applicable to the
existing results of SA (see the discussion in Remark 2 for more details).

#### Two-Time-Scale Stochastic Approximation

is proposed to find the solution $(x^{\star},y^{\star})$ of (1.1) under weaker
conditions of the operators $F,G$ when SA is not applicable. Indeed, in many
applications presented in Section 4, two-time-scale SA is observed to perform
better than its single-time-scale counterpart, either more stable or converge
faster. Unfortunately, this advantage comes at a price: the convergence
analysis of the two-time-scale SA in (1.2), especially its finite-sample
complexity, is more challenging to study than SA. To date, the convergence of
this method, both asymptotic and finite-time convergence rate, is only well-
studied in the linear setting; see for example in [4, 22, 23, 24, 25, 7, 26,
8]. In particular, the asymptotic rate of (1.2) is probably first studied in
[4], where the authors show an ${\cal O}(1/k)$ convergence in distribution via
central limit theorem. Following this work, several results on the convergence
properties of (1.2) have been studied, e.g., convergence with high probability
in [26] and convergence in mean-square errors with constant step sizes in
[24]. The author in [25] extends the work in [24] to time-varying step sizes
and proposes a restarting scheme for step sizes to improve the transient time
of (1.2). Both works in [24, 25] only achieve a rate ${\cal O}(1/k^{2/3})$,
which is considered as suboptimal. The convergence rate of (1.2) in the linear
setting is later shown to be ${\cal O}(1/k)$ in [7], utilizing the setting in
[4]. Recently, this result is further refined with tighter constants in the
upper bounds of the convergence rates in [8].

Unlike the linear setting, theoretical results on the convergence rates of
(1.2) in the general nonlinear setting are less understood. In [5], the
authors provide an asymptotic convergence rate of (1.2) in distribution via
central limit theorem when the noise is i.i.d. This results is recently
generalized to the Markov setting, i.e., noise is dependent and samples are
biased, in [6]. On the other hand, the first finite-time complexity of (1.2)
in the nonlinear setting is studied in [9, 27], where the author shows a rate
of ${\cal O}(1/k^{2/3})$ in the mean-squared regime. This result is shown by
leveraging the singular perturbation technique in control theory [28] that is
first proposed to study the finite-time complexity of linear two-time-scale SA
in [24]. This result is recently refined in [10], where the authors obtain an
${\cal O}(1/k)$ convergence rate. However, the analysis in [10] requires a
rather restrictive condition on the high order of the smoothness of the
operators, which may not hold for any nonlinear operator (see the discussion
in Remark 4). The focus of this work is to develop a new variant of two-time-
scale SA that gives an ${\cal O}(1/k)$ convergence rate without using an
assumption on the high order of the smoothness of the operators.

##  2 Fast Two-Time-Scale Stochastic Approximation

The main difficulty in analyzing the convergence of the classic two-time-scale
SA in (1.2) is to handle the coupling between the two iterates, i.e., the
update of $x$ depends on $y$ and vice versa. A common approach to address this
problem is to choose proper step sizes $\beta_{k}\ll\alpha_{k}$, e.g.,
$\beta_{k}\sim 1/k\leq\alpha_{k}\sim 1/k^{2/3}$. Then one can show that the
update of $x$ is asymptotically independent on $y$ and vice versa. However,
under the impact of stochastic noise due to the sampling of operators, such a
choice of step sizes yields a suboptimal finite-time convergence rate of two-
time-scale SA, e.g., ${\cal O}(1/k^{2/3})$. In this paper, we propose an
alternative approach that will help to decouple the impact of noise on the
update of the iterates. In particular, we consider the following variant of
two-time-scale SA in (1.2)

$\displaystyle\left\\{\begin{array}[]{l}f_{k+1}=(1-\lambda_{k})f_{k}+\lambda_{k}(F(x_{k},y_{k})+\xi_{k})\\\
x_{k+1}=x_{k}-\alpha_{k}f_{k},\end{array}\right.\;\text{and}\;\left\\{\begin{array}[]{l}g_{k+1}=(1-\gamma_{k})g_{k}+\gamma_{k}(G(x_{k},y_{k})+\psi_{k})\\\
y_{k+1}=y_{k}-\beta_{k}g_{k}.\end{array}\right.$ (2.7)

Here, $f_{k}$ and $g_{k}$ are used to estimate the time-weighted average of
the operators by using the two nonnegative step sizes $\gamma_{k}$ and
$\lambda_{k}$. These estimates are then used to update the iterates $x,y$.
When $\lambda_{k}=\gamma_{k}=1$, the updates in (2.7) reduce to the ones in
(1.2). On the other hand, when $\lambda_{k},\gamma_{k}<1$ one can view the
updates of $f_{k}$ and $g_{k}$ as a variant of the classic Ruppert-Polyak
averaging technique [29, 30]. Ruppert-Polyak averaging technique is often
applied to estimate the time-weighted average values of the iterates in
optimization and SA literature. The average variable often achieves a better
performance than the one without using averaging; see for example [14, 8, 31]
and the references therein. In the proposed method in (2.7), we utilize the
Ruppert-Polyak averaging technique on the operators, which will help to
decompose the coupling between sampling noise and the iterate updates. Indeed,
we will show that by choosing $\lambda_{k},\gamma_{k}<1$, one can separate the
impact of the noise on the updates of the iterates, therefore, achieving an
optimal rate ${\cal O}(1/k)$ under a proper choice of step sizes
$\alpha_{k},\beta{k}$.

######  Remark 1.

In our approach, we use the Ruppert-Polyak averaging steps on the samples of
the operators, which is different to the existing works where these steps are
often applied to the main iterates. The benefit of using this averaging
technique to the iterates has been observed in the literature. On the other
hand, we are not aware of any prior work that uses this technique to improve
the estimate of the operators in parallel with the iterate updates. An
interesting question is to study the performance of these two types of
averaging steps when they are implemented in parallel.

###  2.1 Notation and Main Assumptions

In this section, we present the main assumptions which we will be used to
study the finite-time convergence rate of (2.7). These assumptions are fairly
standard and have been used in prior works under different variants [4, 22,
23, 24, 25, 7, 5]. Our first assumption is on the smoothness of the operators.

######  Assumption 1.

Given $y\in\mathbb{R}^{d}$ there exists an operator
$H:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ such that $x=H(y)$ satisfies

$\displaystyle F(H(y),y)=0.$

In addition, we assume that the operators $F,G,H$ are Lipschitz continuous
with a positive constant $L$, i.e., the following conditions satisfy for all
$x_{1},x_{2},y_{1},y_{2}\in\mathbb{R}^{d}$

$\displaystyle\begin{aligned} \|H(y_{1})-H(y_{2})\|&\leq L\|y_{1}-y_{2}\|,\\\
\|F(x_{1},y_{1})-F(x_{2},y_{2})\|&\leq L(\|x_{1}-x_{2}\|+\|y_{1}-y_{2}\|),\\\
\|G(x_{1},y_{1})-G(x_{2},y_{2})\|&\leq
L\left(\|x_{1}-x_{2}\|+\|y_{1}-y_{2}\|\right).\end{aligned}$ (2.8)

Second, we will study the setting when $F,G$ are strongly monotone, which
implies that $(x^{\star},y^{\star})$ is unique.

######  Assumption 2.

Given $y$, $F$ is strongly monotone w.r.t $x$, i.e., there exists a constant
$\mu_{F}>0$ such that

$\displaystyle\left\langle
x-z,F(x,y)-F(z,y)\right\rangle\geq\mu_{F}\|x-z\|^{2},\quad\forall
x,z\in\mathbb{R}^{d}.$ (2.9)

Moreover, $G$ is $1$-point strongly monotone w.r.t $y^{\star}$, i.e., there
exists a constant $\mu_{G}\in(0,\mu_{F}]$ such that

$\displaystyle\left\langle
y-y^{\star},G(H(y),y)\right\rangle\geq\mu_{G}\|y-y^{\star}\|^{2},\quad\forall
y\in\mathbb{R}^{d}.$ (2.10)

######  Remark 2.

First, Assumption 1 is often referred to as the global Lipschitz continuity.
Our analysis in this paper can be easily extended to the local Lipschitz
condition as studied in [6] (with extra constant terms in our theoretical
bounds), i.e., $H$ (similarly, $F,G$) satisfies

$\displaystyle\|H(y)\|\leq L(1+\|y\|)\quad\forall y\in\mathbb{R}^{d}.$

Second, it is also worth to point out that in general Assumption 2 does not
imply that the operator $W=[F,G]$ is strongly monotone. Indeed, if $W$ is
strongly monotone, one can apply the single-time-scale SA to solve (1.1),
e.g., by considering SA to update $z=[x,y]$. This statement can be verified in
the linear setting, where $F$ and $G$ are linear operators

$\displaystyle F(x,y)$ $\displaystyle={\bf A}_{11}x+{\bf A}_{12}y-b_{1}$
$\displaystyle G(x,y)$ $\displaystyle={\bf A}_{12}x+{\bf A}_{22}y-b_{2}.$

In this case, Assumption 2 implies that ${\bf A}_{11}$ and $\Delta={\bf
A}_{22}-{\bf A}_{21}{\bf A}_{11}^{-1}{\bf A}_{12}$ are negative definite (but
not necessarily symmetric), and $H(y)={\bf A}_{11}^{-1}(b_{1}-{\bf A}_{12}y)$.
It is straightforward to see that $(x^{*},y^{*})$ is unique and satisfies

$\displaystyle x^{*}$ $\displaystyle={\bf A}_{11}^{-1}(b_{1}-{\bf
A}_{12}y^{*})$ (2.11) $\displaystyle y^{*}$
$\displaystyle=\Delta^{-1}\left(b_{2}-{\bf A}_{21}{\bf
A}_{11}^{-1}b_{1}\right).$ (2.12)

However, the matrix ${\bf A}$ defined below is not negative definite

$\displaystyle{\bf A}=\left[\begin{array}[]{cc}{\bf A}_{11}&{\bf A}_{12}\\\
{\bf A}_{12}&{\bf A}_{22}\end{array}\right].$

To see this, we consider the following example

$\displaystyle{\bf A}_{11}$
$\displaystyle=\left[\begin{array}[]{cc}-3.5&-8.5\\\
1&0\end{array}\right],\quad{\bf A}_{12}=\left[\begin{array}[]{cc}2&0\\\
0&2\end{array}\right],\quad{\bf A}_{21}$
$\displaystyle=\left[\begin{array}[]{cc}5&0\\\ 0&5\end{array}\right],\quad{\bf
A}_{22}=\left[\begin{array}[]{cc}-10&-1.6\\\ 20&-1\end{array}\right],$

Then the eigenvalues of ${\bf A}_{11}$, $\Delta$, and ${\bf A}$ are

$\displaystyle\text{eig}({\bf A}_{11})$ $\displaystyle=-1.75\pm
2.33i,\quad\text{eig}(\Delta)=-3.44\pm 14.24i$ $\displaystyle\text{eig}({\bf
A})$ $\displaystyle=\\{-8.53\pm 4.51i,\;1.28\pm 4.24i\\},$

which obviously shows that ${\bf A}_{11}$ and $\Delta$ are negative but ${\bf
A}$ is not.

Finally, we consider i.i.d noise models, that is, $\xi_{k}$ and $\psi_{k}$ are
Martingale difference. We denote by ${\cal Q}_{k}$ the filtration containing
all the history generated by (1.2) up to time $k$, i.e.,

${\cal
Q}_{k}=\\{x_{0},y_{0},\xi_{0},\psi_{0},\xi_{1},\psi_{1},\ldots,\xi_{k},\psi_{k}\\}.$

######  Assumption 3.

The random variables $\xi_{k}$ and $\psi_{k}$, for all $k\geq 0$, are
independent of each other and across time, with zero mean and common variances
given as follows

$\mathbb{E}[\xi_{k}^{T}\xi_{k}\,|\,{\cal
Q}_{k-1}]=\Gamma_{11},\quad\mathbb{E}[\psi_{k}^{T}\psi_{k}\,|\,{\cal
Q}_{k-1}]=\Gamma_{22}.$ (2.13)

######  Remark 3.

We note that our analysis in this paper can be extended to cover the Markov
setting where $\xi_{k}$ and $\psi_{k}$ are generated by Markov processes. In
this case, $\xi_{k}$ and $\psi_{k}$ are dependent and the samples are biased.
In this setting, if the Markov processes are ergodic, i.e., they have
geometric mixing time, one can apply the technique proposed in [18] to handle
the dependent and biased samples. This will result in an extra factor ${\cal
O}(\log(k))$ in our theoretic bounds, where $k$ is the number of iterations.

##  3 Main Results

In this section, we present in details the main results of this paper, that
is, we provide a finite-time analysis for the convergence rates of the two-
time-scale stochastic approximation in (2.7). Under assumptions presented in
Section 2.1, we show that the mean square errors of these iterates,
$\mathbb{E}[\|x_{k}-x^{*}\|^{2}]$ and $\mathbb{E}[\|y_{k}-x^{*}\|^{2}]$,
converge to zero at a rate ${\cal O}(1/k)$. To show this result, we introduce
the following four residual variables,

$\displaystyle\begin{aligned} \Delta
f_{k}=f_{k}-F(x_{k},y_{k})\quad\text{and}\quad{\hat{x}}_{k}&=x_{k}-H(y_{k}),\\\
\Delta
g_{k}=g_{k}-G(x_{k},y_{k})\quad\text{and}\quad{\hat{y}}_{k}&=y_{k}-y^{\star}.\end{aligned}$
(3.14)

One can rewrite the updates of $x,y$ in (2.7) as

$\displaystyle x_{k+1}=x_{k}-\alpha_{k}F(x_{k},y_{k})-\alpha_{k}\Delta
f_{k}\quad\text{and}\quad y_{k+1}-\beta_{k}G(x_{k},y_{k})-\beta_{k}\Delta
g_{k},$

Thus, if $(\Delta f_{k},\Delta g_{k})$ goes to zero then the impact of
sampling noise on $x,y$ will be negligible. In addition, if
$({\hat{x}}_{k},{\hat{y}}_{k})$ go to zero then
$(x_{k},y_{k})\rightarrow(x^{\star},y^{\star})$. Thus, to establish the
convergence of $(x_{k},y_{k})$ to $(x^{\star},y^{\star})$ one can instead
study the convergence of $(\Delta f_{k},\Delta g_{k})$ and
$({\hat{x}}_{k},{\hat{y}}_{k})$ to zero. Indeed, our main result presented in
Theorem 1 below is to study an upper bound for the convergence rate of these
quantities. We will do this by considering the following choice of step sizes

$\displaystyle\gamma_{k}=\lambda_{k}=\frac{1}{4(k+1)},\quad\alpha_{k}=\frac{C_{\alpha}}{k+1},\quad\beta_{k}=\frac{C_{\beta}}{k+h+1},$
(3.15)

where $C_{\beta},C_{\alpha},h$ are some positive constants, which will be
specified in Theorem 1. Here, we choose $\lambda_{k}=\gamma_{k}$ for an ease
of exposition. In addition, we will consider the setting when
$\beta_{k}\leq\alpha_{k}$, i.e., the iterate $x_{k}$ is implemented at a
faster time scale than $y_{k}$.

###  3.1 Preliminaries

We first present some preliminary results, which will help to facilitate the
analysis of Theorem 1. In particular, the following lemmas are to provide some
important bounds for the size of the variables in our updates.

######  Lemma 1.

Let $\\{f_{k},g_{k},x_{k},y_{k}\\}$ be generated by (2.7). Then we have

$\displaystyle\begin{aligned} \|g_{k}\|&\leq\|\Delta
g_{k}\|+L\|{\hat{x}}_{k}\|+(1+L)^{2}\|{\hat{y}}_{k}\|,\\\
\|f_{k}\|&\leq\|\Delta f_{k}\|+L\|{\hat{x}}_{k}\|,\\\
\|x_{k+1}-x_{k}\|&\leq\alpha_{k}(\|\Delta f_{k}\|+L\|{\hat{x}}_{k}\|),\\\
\|y_{k+1}-y_{k}\|&\leq\beta_{k}(\|\Delta
g_{k}\|+L\|{\hat{x}}_{k}\|+(1+L)^{2}\|{\hat{y}}_{k}\|).\end{aligned}$ (3.16)

###### Proof.

Using the Lipschitz continuity of $G,H$ and (3.14) we first have

$\displaystyle\|g_{k}\|$ $\displaystyle=\|\Delta
g_{k}+G(x_{k},y_{k})-G(H(y_{k}),y_{k})+G(H(y_{k}),y_{k})-G(H(y^{\star}),y^{\star})\|$
$\displaystyle\leq\|\Delta
g_{k}\|+L\|{\hat{x}}_{k}\|+L(1+L)\|{\hat{y}}_{k}\|\leq\|\Delta
g_{k}\|+L\|{\hat{x}}_{k}\|+(1+L)^{2}\|{\hat{y}}_{k}\|.$

Second, since $F(H(y_{k}),y_{k})=0$ we have

$\displaystyle\|f_{k}\|=\|\Delta
f_{k}+F(x_{k},y_{k})-F(H(y_{k}),y_{k})\|\leq\|\Delta
f_{k}\|+L\|{\hat{x}}_{k}\|.$

Finally, by using Eqs. (2.7) the last two inequalities in (3.16) can be
obtained from the relations above. ∎

######  Lemma 2.

Let $\\{f_{k},g_{k},x_{k},y_{k}\\}$ be generated by (2.7), and
$\\{\lambda_{k}\\}$ be chosen as in (3.15). Then we have

$\displaystyle\|\Delta f_{k+1}\|^{2}$
$\displaystyle\leq(1-\lambda_{k})\|\Delta
f_{k}\|^{2}-\Big{(}\frac{1}{4}-\frac{8L^{2}\alpha_{k}^{2}}{\lambda_{k}^{2}}\Big{)}\lambda_{k}\|\Delta
f_{k}\|^{2}+2(1-\lambda_{k})\lambda_{k}\epsilon_{k}^{T}\Delta f_{k}$
$\displaystyle\quad+2\lambda_{k}^{2}\|\epsilon_{k}\|^{2}+\frac{8L^{2}}{\lambda_{k}}\Big{(}\alpha_{k}^{2}\|\Delta
g_{k}\|^{2}+2L^{2}\alpha_{k}^{2}\|{\hat{x}}_{k}\|^{2}+(1+L)^{4}\beta_{k}^{2}\|{\hat{y}}_{k}\|^{2}\Big{)}.$
(3.17)

###### Proof.

Recall that $\Delta f_{k}=f_{k}-F(x_{k},y_{k})$. By (2.7) we consider

$\displaystyle\Delta f_{k+1}$
$\displaystyle=(1-\lambda_{k})f_{k}+\lambda_{k}F(x_{k},y_{k})-F(x_{k+1},y_{k+1})+\lambda_{k}\epsilon_{k}$
$\displaystyle=(1-\lambda_{k})\Delta
f_{k}+F(x_{k},y_{k})-F(x_{k+1},y_{k+1})+\lambda_{k}\epsilon_{k},$

which by taking the quadratic norm on both sides gives

$\displaystyle\|\Delta f_{k+1}\|^{2}$
$\displaystyle=(1-\lambda_{k})^{2}\|\Delta
f_{k}\|^{2}+\|F(x_{k},y_{k})-F(x_{k+1},y_{k+1})+\lambda_{k}\epsilon_{k}\|^{2}$
$\displaystyle\quad+2(1-\lambda_{k})\lambda_{k}\epsilon_{k}^{T}\Delta
f_{k}+2(1-\lambda_{k})\Delta f_{k}^{T}(F(x_{k},y_{k})-F(x_{k+1},y_{k+1}))$
$\displaystyle\leq(1-\lambda_{k})^{2}\|\Delta
f_{k}\|^{2}+2\|F(x_{k},y_{k})-F(x_{k+1},y_{k+1})\|^{2}+2\lambda_{k}^{2}\|\epsilon_{k}\|^{2}$
$\displaystyle\quad+2(1-\lambda_{k})\lambda_{k}\epsilon_{k}^{T}\Delta
f_{k}+\frac{\lambda_{k}}{2}\|\Delta
f_{k}\|^{2}+\frac{2}{\lambda_{k}}\|F(x_{k},y_{k})-F(x_{k+1},y_{k+1})\|^{2}$
$\displaystyle=(1-\lambda_{k})\|\Delta
f_{k}\|^{2}-\Big{(}\frac{1}{2}-\lambda_{k}\Big{)}\lambda_{k}\|\Delta
f_{k}\|^{2}+2\lambda_{k}^{2}\|\epsilon_{k}\|^{2}+2(1-\lambda_{k})\lambda_{k}\epsilon_{k}^{T}\Delta
f_{k}$
$\displaystyle\quad+\frac{2}{\lambda_{k}}\|F(x_{k},y_{k})-F(x_{k+1},y_{k+1})\|^{2},$
(3.18)

where the inequality is due the Cauchy-Schwarz inequality and $\lambda_{k}<1$.
By (3.16) and using $\beta_{k}\leq\alpha_{k}$ we have

$\displaystyle\|F(x_{k},y_{k})-F(x_{k+1},y_{k+1})\|^{2}\leq
4L^{2}\Big{(}\alpha_{k}^{2}\|\Delta f_{k}\|^{2}+\alpha_{k}^{2}\|\Delta
g_{k}\|^{2}+2L^{2}\alpha_{k}^{2}\|{\hat{x}}_{k}\|^{2}+(1+L)^{4}\beta_{k}^{2}\|{\hat{y}}_{k}\|^{2}\Big{)},$

which when substituting into (3.18) and using $\lambda_{k}\leq 1/4$ we obtain
(3.17), i.e.,

$\displaystyle\|\Delta f_{k+1}\|^{2}$
$\displaystyle\leq(1-\lambda_{k})\|\Delta
f_{k}\|^{2}-\Big{(}\frac{1}{2}-\lambda_{k}\Big{)}\lambda_{k}\|\Delta
f_{k}\|^{2}+2\lambda_{k}^{2}\|\epsilon_{k}\|^{2}+2(1-\lambda_{k})\lambda_{k}\epsilon_{k}^{T}\Delta
f_{k}$
$\displaystyle\quad+\frac{8L^{2}}{\lambda_{k}}\Big{(}\alpha_{k}^{2}\|\Delta
f_{k}\|^{2}+\alpha_{k}^{2}\|\Delta
g_{k}\|^{2}+2L^{2}\alpha_{k}^{2}\|{\hat{x}}_{k}\|^{2}+(1+L)^{4}\beta_{k}^{2}\|{\hat{y}}_{k}\|^{2}\Big{)}$
$\displaystyle\leq(1-\lambda_{k})\|\Delta
f_{k}\|^{2}-\Big{(}\frac{1}{4}-\frac{8L^{2}\alpha_{k}^{2}}{\lambda_{k}^{2}}\Big{)}\lambda_{k}\|\Delta
f_{k}\|^{2}+2(1-\lambda_{k})\lambda_{k}\epsilon_{k}^{T}\Delta f_{k}$
$\displaystyle\quad+2\lambda_{k}^{2}\|\epsilon_{k}\|^{2}+\frac{8L^{2}}{\lambda_{k}}\Big{(}\alpha_{k}^{2}\|\Delta
g_{k}\|^{2}+2L^{2}\alpha_{k}^{2}\|{\hat{x}}_{k}\|^{2}+(1+L)^{4}\beta_{k}^{2}\|{\hat{y}}_{k}\|^{2}\Big{)}.$

∎

######  Lemma 3.

Let $\\{f_{k},g_{k},x_{k},y_{k}\\}$ be generated by (2.7), and
$\\{\gamma_{k}\\}$ be chosen as in (3.15). Then we have

$\displaystyle\|\Delta g_{k+1}\|^{2}$ $\displaystyle\leq(1-\gamma_{k})\|\Delta
g_{k}\|^{2}-\Big{(}\frac{1}{4}-\frac{8L^{2}\alpha_{k}^{2}}{\gamma_{k}^{2}}\Big{)}\gamma_{k}\|\Delta
g_{k}\|^{2}+2(1-\gamma_{k})\gamma_{k}\psi_{k}^{T}\Delta g_{k}$
$\displaystyle\quad+2\gamma_{k}^{2}\|\psi_{k}\|^{2}+\frac{8L^{2}}{\gamma_{k}}\Big{(}\alpha_{k}^{2}\|\Delta
f_{k}\|^{2}+2L^{2}\alpha_{k}^{2}\|{\hat{x}}_{k}\|^{2}+(1+L)^{4}\beta_{k}^{2}\|{\hat{y}}_{k}\|^{2}\Big{)}.$
(3.19)

###### Proof.

Recall that $\Delta g_{k}=g_{k}-G(x_{k,y_{k}})$. By (2.7) we have

$\displaystyle\Delta g_{k+1}=(1-\gamma_{k})\Delta
g_{k}+G(x_{k},y_{k})-G(x_{k+1},y_{k+1})+\gamma_{k}\psi_{k},$

which by taking the quadratic norm on both sides gives

$\displaystyle\|\Delta g_{k+1}\|^{2}$
$\displaystyle=(1-\gamma_{k})^{2}\|\Delta
g_{k}\|^{2}+\|G(x_{k},y_{k})-G(x_{k+1},y_{k+1})+\gamma_{k}\psi_{k}\|^{2}$
$\displaystyle\quad+2(1-\gamma_{k})\gamma_{k}\psi_{k}^{T}\Delta
g_{k}+2(1-\gamma_{k})\Delta g_{k}^{T}(G(x_{k},y_{k})-G(x_{k+1},y_{k+1}))$
$\displaystyle\leq(1-\gamma_{k})^{2}\|\Delta
G_{k}\|^{2}+2\|G(x_{k},y_{k})-G(x_{k+1},y_{k+1})\|^{2}+2\gamma_{k}^{2}\|\psi_{k}\|^{2}$
$\displaystyle\quad+2(1-\gamma_{k})\gamma_{k}\psi_{k}^{T}\Delta
g_{k}+\frac{\gamma_{k}}{2}\|\Delta
g_{k}\|^{2}+\frac{2}{\gamma_{k}}\|G(x_{k},y_{k})-G(x_{k+1},y_{k+1})\|^{2}$
$\displaystyle=(1-\gamma_{k})\|\Delta
g_{k}\|^{2}-\Big{(}\frac{1}{2}-\gamma_{k}\Big{)}\gamma_{k}\|\Delta
g_{k}\|^{2}+2\gamma_{k}^{2}\|\psi_{k}\|^{2}+2(1-\gamma_{k})\gamma_{k}\psi_{k}^{T}\Delta
g_{k}$
$\displaystyle\quad+\frac{2}{\gamma_{k}}\|G(x_{k},y_{k})-G(x_{k+1},y_{k+1})\|^{2},$
(3.20)

where in the inequality is due to the Cauchy-Schwarz inequality and
$\gamma_{k}<1$. By (3.16) and using $\beta_{k}\leq\alpha_{k}$ we have

$\displaystyle\|G(x_{k},y_{k})-G(x_{k+1},y_{k+1})\|^{2}\leq
4L^{2}\Big{(}\alpha_{k}^{2}\|\Delta f_{k}\|^{2}+\alpha_{k}^{2}\|\Delta
g_{k}\|^{2}+2L^{2}\alpha_{k}^{2}\|{\hat{x}}_{k}\|^{2}+(1+L)^{4}\beta_{k}^{2}\|{\hat{y}}_{k}\|^{2}\Big{)},$

which when substituting into (3.20) and using $\beta_{k}\leq 1/4$ we obtain
(3.19), i.e.,

$\displaystyle\|\Delta g_{k+1}\|^{2}$ $\displaystyle\leq(1-\gamma_{k})\|\Delta
g_{k}\|^{2}-\Big{(}\frac{1}{2}-\gamma_{k}\Big{)}\gamma_{k}\|\Delta
g_{k}\|^{2}+2\gamma_{k}^{2}\|\psi_{k}\|^{2}+2(1-\gamma_{k})\gamma_{k}\psi_{k}^{T}\Delta
g_{k}$
$\displaystyle\quad+\frac{8L^{2}}{\gamma_{k}}\Big{(}\alpha_{k}^{2}\|\Delta
f_{k}\|^{2}+\alpha_{k}^{2}\|\Delta
g_{k}\|^{2}+2L^{2}\alpha_{k}^{2}\|{\hat{x}}_{k}\|^{2}+(1+L)^{4}\beta_{k}^{2}\|{\hat{y}}_{k}\|^{2}\Big{)}$
$\displaystyle\leq(1-\gamma_{k})\|\Delta
g_{k}\|^{2}-\Big{(}\frac{1}{4}-\frac{8L^{2}\alpha_{k}^{2}}{\gamma_{k}^{2}}\Big{)}\gamma_{k}\|\Delta
g_{k}\|^{2}+2(1-\gamma_{k})\gamma_{k}\psi_{k}^{T}\Delta g_{k}$
$\displaystyle\quad+2\gamma_{k}^{2}\|\psi_{k}\|^{2}+\frac{8L^{2}}{\gamma_{k}}\Big{(}\alpha_{k}^{2}\|\Delta
f_{k}\|^{2}+2L^{2}\alpha_{k}^{2}\|{\hat{x}}_{k}\|^{2}+(1+L)^{4}\beta_{k}^{2}\|{\hat{y}}_{k}\|^{2}\Big{)}.$

∎

######  Lemma 4.

Let $\\{f_{k},g_{k},x_{k},y_{k}\\}$ be generated by (2.7), and $\alpha_{k}$ be
chosen as in (3.15) with $\alpha_{0}$ satisfies

$\displaystyle\alpha_{0}\leq\frac{\mu_{F}}{L^{2}}\cdot$ (3.21)

Then we have

$\displaystyle\|{\hat{x}}_{k+1}\|^{2}$
$\displaystyle\leq\left(1-\mu_{F}\alpha_{k}\right)\|{\hat{x}}_{k}\|^{2}-\left[\frac{\mu_{F}}{2}\alpha_{k}-\Big{(}3L+\frac{(1+L)^{4}}{\mu_{G}}\Big{)}\beta_{k}-\Big{(}3L^{4}+2L^{2}\Big{)}\alpha_{k}^{2}\right]\|{\hat{x}}_{k}\|^{2}$
$\displaystyle\quad+\frac{\mu_{G}}{4}\beta_{k}\|{\hat{y}}_{k}\|^{2}+(1+L)^{4}\big{(}3L^{2}\beta_{k}+\alpha_{k}\big{)}\beta_{k}\|{\hat{y}}_{k}\|^{2}$
$\displaystyle\quad+\Big{[}\frac{2}{\mu_{F}}\alpha_{k}+4\alpha_{k}^{2}\Big{]}\|\Delta
f_{k}\|^{2}+\Big{[}\frac{\alpha_{k}}{L}+(3L^{2}+1)\alpha_{k}^{2}\Big{]}\|\Delta
g_{k}\|^{2}.$ (3.22)

###### Proof.

We consider

$\displaystyle{\hat{x}}_{k+1}$
$\displaystyle=x_{k+1}-H(y_{k+1})=x_{k}-\alpha_{k}f_{k}-H(y_{k+1})$
$\displaystyle={\hat{x}}_{k}-\alpha_{k}F(x_{k},y_{k})-\alpha_{k}\Delta
f_{k}+H(y_{k})-H\left(y_{k}-\beta_{k}g_{k}\right),$

which implies

$\displaystyle\|{\hat{x}}_{k+1}\|^{2}$
$\displaystyle=\left\|{\hat{x}}_{k}-\alpha_{k}F(x_{k},y_{k})-\alpha_{k}\Delta
f_{k}+H(y_{k})-H\left(y_{k}-\beta_{k}g_{k}\right)\right\|^{2}$
$\displaystyle=\|{\hat{x}}_{k}-\alpha_{k}F(x_{k},y_{k})\|^{2}+\left\|H(y_{k})-H(y_{k}-\beta_{k}g_{k})\right\|^{2}+\alpha_{k}^{2}\|\Delta
f_{k}\|^{2}$
$\displaystyle\qquad+2\left({\hat{x}}_{k}-\alpha_{k}F(x_{k},y_{k})\right)^{T}\left(H(y_{k})-H(y_{k}-\beta_{k}g_{k})\right)-2\alpha_{k}\left({\hat{x}}_{k}-\alpha_{k}F(x_{k},y_{k})\right)^{T}\Delta
f_{k}$
$\displaystyle\qquad-2\alpha_{k}\left(H(y_{k})-H(y_{k}-\beta_{k}g_{k})\right)^{T}\Delta
f_{k}$ $\displaystyle=A_{1}+A_{2}+\alpha_{k}^{2}\|\Delta
f_{k}\|^{2}+A_{3}+A_{4}+A_{5},$ (3.23)

where $A_{1},\ldots,A_{5}$ are the terms on the right-hand side of the second
equation, in that order. We next analyze these terms. First, using
$F(H(y_{k}),y_{k})=0$ and ${\hat{x}}_{k}=x_{k}-H(y_{k})$ we have

$\displaystyle
A_{1}=\|{\hat{x}}_{k}-\alpha_{k}F(x_{k},y_{k})\|^{2}=\|{\hat{x}}_{k}\|^{2}-2\alpha_{k}{\hat{x}}_{k}^{T}F(x_{k},y_{k})+\|\alpha_{k}F(x_{k},y_{k})\|^{2}$
$\displaystyle=\|{\hat{x}}_{k}\|^{2}-2\alpha_{k}(x_{k}-H(y_{k}))^{T}\left(F(x_{k},y_{k})-F(H(y_{k}),y_{k})\right)+\alpha_{k}^{2}\|F(x_{k},y_{k})-F(H(y_{k}),y_{k})\|^{2}$
$\displaystyle\leq\|{\hat{x}}_{k}\|^{2}-2\mu_{F}\alpha_{k}\|x_{k}-H(y_{k})\|^{2}+L^{2}\alpha_{k}^{2}\|x_{k}-H(y_{k})\|^{2}$
$\displaystyle=\left(1-2\mu_{F}\alpha_{k}+L^{2}\alpha_{k}^{2}\right)\|{\hat{x}}_{k}\|^{2},$
(3.24)

where in the first inequality we use the strong monotonicity and Lispchitz
continuity of $F$. Second, using the Lipschitz continuity of $H$ and (3.16) we
obtain

$\displaystyle A_{2}$
$\displaystyle=\left\|H(y_{k})-H(y_{k}-\beta_{k}g_{k})\right\|^{2}\leq
L^{2}\beta_{k}^{2}\|g_{k}\|^{2}\leq L^{2}\beta_{k}^{2}\big{(}\|\Delta
g_{k}\|+L\|{\hat{x}}_{k}\|+L(1+L)\|{\hat{y}}_{k}\|\big{)}^{2}$
$\displaystyle\leq 3L^{2}\beta_{k}^{2}\big{(}\|\Delta
g_{k}\|^{2}+L^{2}\|{\hat{x}}_{k}\|^{2}+L^{2}(1+L)^{2}\|{\hat{y}}_{k}\|^{2}\big{)}.$
(3.25)

Third, using $\alpha_{k}\leq\alpha_{0}\leq\mu_{F}/L^{2}$ into (3.24) gives
$\|A_{1}\|\leq\|{\hat{x}}_{k}\|^{2}$, which by (3.16) and
$F(H(y_{k}),y_{k})=0$ gives

$\displaystyle A_{3}$
$\displaystyle=2\left({\hat{x}}_{k}-\alpha_{k}F(x_{k},y_{k})\right)^{T}\left(H(y_{k})-H(y_{k}-\beta_{k}g_{k})\right)\leq
2\beta_{k}\|{\hat{x}}_{k}-\alpha_{k}F(x_{k},y_{k})\|\|g_{k}\|$
$\displaystyle\leq 2\beta_{k}\|{\hat{x}}_{k}\|\|g_{k}\|\leq
2\beta_{k}\|{\hat{x}}_{k}\|\big{(}\|\Delta
g_{k}\|+L\|{\hat{x}}_{k}\|+(1+L)^{2}\|{\hat{y}}_{k}\|\big{)}$
$\displaystyle=2L\beta_{k}\|{\hat{x}}_{k}\|^{2}+2\beta_{k}\|{\hat{x}}_{k}\|\|\Delta
g_{k}\|+2(1+L)^{2}\beta_{k}\|{\hat{x}}_{k}\|\|{\hat{y}}_{k}\|.$ (3.26)

Similarly, we consider $A_{4}$

$\displaystyle
A_{4}=-2\alpha_{k}\left({\hat{x}}_{k}-\alpha_{k}F(x_{k},y_{k})\right)^{T}\Delta
f_{k}\leq 2\alpha_{k}\|{\hat{x}}_{k}\|\|\Delta f_{k}\|.$ (3.27)

Finally, using (3.16) one more time we have

$\displaystyle A_{5}$
$\displaystyle=-2\alpha_{k}\left(H(y_{k})-H(y_{k}-\beta_{k}g_{k})\right)^{T}\Delta
f_{k}\leq 2\alpha_{k}\beta_{k}\|g_{k}\|\|\Delta f_{k}\|$ $\displaystyle\leq
2\alpha_{k}\beta_{k}\big{(}\|\Delta
g_{k}\|+L\|{\hat{x}}_{k}\|+(1+L)^{2}\|{\hat{y}}_{k}\|\big{)}\|\Delta f_{k}\|.$
(3.28)

Substituting Eqs. (3.24)–(3.28) into (3.23) we obtain (3.22), i.e.,

$\displaystyle\|{\hat{x}}_{k+1}\|^{2}$
$\displaystyle\leq\left(1-2\mu_{F}\alpha_{k}+L^{2}\alpha_{k}^{2}\right)\|{\hat{x}}_{k}\|^{2}+3L^{2}\beta_{k}^{2}\big{(}\|\Delta
g_{k}\|^{2}+L^{2}\|{\hat{x}}_{k}\|^{2}+L^{2}(1+L)^{2}\|{\hat{y}}_{k}\|^{2}\big{)}+\alpha_{k}^{2}\|\Delta
f_{k}\|^{2}$
$\displaystyle\quad+2L\beta_{k}\|{\hat{x}}_{k}\|^{2}+2\beta_{k}\|{\hat{x}}_{k}\|\|\Delta
g_{k}\|+2(1+L)^{2}\beta_{k}\|{\hat{x}}_{k}\|\|{\hat{y}}_{k}\|$
$\displaystyle\quad+2\alpha_{k}\|{\hat{x}}_{k}\|\|\Delta
f_{k}\|+2\alpha_{k}\beta_{k}\big{(}\|\Delta
g_{k}\|+L\|{\hat{x}}_{k}\|+(1+L)^{2}\|{\hat{y}}_{k}\|\big{)}\|\Delta f_{k}\|$
$\displaystyle\leq\left(1-2\mu_{F}\alpha_{k}+L^{2}\alpha_{k}^{2}\right)\|{\hat{x}}_{k}\|^{2}+3L^{2}\beta_{k}^{2}\big{(}\|\Delta
g_{k}\|^{2}+L^{2}\|{\hat{x}}_{k}\|^{2}+L^{2}(1+L)^{2}\|{\hat{y}}_{k}\|^{2}\big{)}+\alpha_{k}^{2}\|\Delta
f_{k}\|^{2}$
$\displaystyle\quad+2L\beta_{k}\|{\hat{x}}_{k}\|^{2}+L\beta_{k}\|{\hat{x}}_{k}\|^{2}+\frac{\beta_{k}}{L}\|\Delta
g_{k}\|^{2}+\frac{(1+L)^{4}}{\mu_{G}}\beta_{k}\|{\hat{x}}_{k}\|^{2}+\frac{\mu_{G}}{4}\beta_{k}\|{\hat{y}}_{k}\|^{2}$
$\displaystyle\quad+\frac{\mu_{F}}{2}\alpha_{k}\|{\hat{x}}_{k}\|^{2}+\frac{2}{\mu_{F}}\alpha_{k}\|\Delta
f_{k}\|^{2}+L^{2}\alpha_{k}\beta_{k}\|{\hat{x}}_{k}\|^{2}+3\alpha_{k}\beta_{k}\|\Delta
f_{k}\|^{2}$ $\displaystyle\quad+\alpha_{k}\beta_{k}\big{(}\|\Delta
g_{k}\|^{2}+(1+L)^{4}\|\|{\hat{y}}_{k}\|^{2}\big{)}$
$\displaystyle\leq\left(1-\mu_{F}\alpha_{k}\right)\|{\hat{x}}_{k}\|^{2}-\left[\frac{\mu_{F}}{2}\alpha_{k}-\Big{(}3L+\frac{(1+L)^{4}}{\mu_{G}}\Big{)}\beta_{k}-\Big{(}3L^{4}+2L^{2}\Big{)}\alpha_{k}^{2}\right]\|{\hat{x}}_{k}\|^{2}$
$\displaystyle\quad+\frac{\mu_{G}}{4}\beta_{k}\|{\hat{y}}_{k}\|^{2}+(1+L)^{4}\big{(}3L^{2}\beta_{k}+\alpha_{k}\big{)}\beta_{k}\|{\hat{y}}_{k}\|^{2}$
$\displaystyle\quad+\Big{[}\frac{2}{\mu_{F}}\alpha_{k}+4\alpha_{k}^{2}\Big{]}\|\Delta
f_{k}\|^{2}+\Big{[}\frac{\alpha_{k}}{L}+(3L^{2}+1)\alpha_{k}^{2}\Big{]}\|\Delta
g_{k}\|^{2},$

where the last inequality we use $\beta_{k}\leq\alpha_{k}$. ∎

######  Lemma 5.

Let $\\{f_{k},g_{k},x_{k},y_{k}\\}$ be generated by (2.7), and $\alpha_{k}$ be
chosen as in (3.15) with $\alpha_{0}$ satisfies

$\displaystyle\beta_{0}\leq\frac{\mu_{G}}{4L^{2}(1+L)^{2}}\cdot$ (3.29)

Then we have

$\displaystyle\|{\hat{y}}_{k+1}\|^{2}$
$\displaystyle\leq(1-\mu_{G}\beta_{k})\|{\hat{y}}_{k}\|^{2}-\frac{\mu_{G}}{4}\beta_{k}\|{\hat{y}}_{k}\|^{2}+2\Big{(}\frac{2}{\mu_{G}}+\beta_{k}\Big{)}\beta_{k}\big{(}L^{2}\|{\hat{x}}_{k}\|^{2}+\|\Delta
g_{k}\|^{2}\big{)}.$ (3.30)

###### Proof.

We consider

$\displaystyle{\hat{y}}_{k+1}$
$\displaystyle=y_{k+1}-y^{\star}={\hat{y}}_{k}-\beta_{k}g_{k}={\hat{y}}_{k}-\beta_{k}G(H(y_{k}),y_{k})+\beta_{k}(G(H(y_{k}),y_{k})-G(x_{k},y_{k}))-\beta_{k}\Delta
g_{k},$

which gives

$\displaystyle\|{\hat{y}}_{k+1}\|^{2}$
$\displaystyle=\|{\hat{y}}_{k}-\beta_{k}G(H(y_{k}),y_{k})+\beta_{k}(G(H(y_{k}),y_{k})-G(x_{k},y_{k}))-\beta_{k}\Delta
g_{k}\|^{2}$
$\displaystyle=\|{\hat{y}}_{k}-\beta_{k}G(H(y_{k}),y_{k})\|^{2}+2\beta_{k}^{2}\|\Delta
g_{k}\|^{2}+2\beta_{k}^{2}\|G(H(y_{k}),y_{k})-G(x_{k},y_{k})\|^{2}$
$\displaystyle\quad+2\beta_{k}\big{(}{\hat{y}}_{k}-\beta_{k}G(H(y_{k}),y_{k})\big{)}^{T}\big{(}G(H(y_{k}),y_{k})-G(x_{k},y_{k})-\Delta
g_{k}\big{)}.$ (3.31)

First, using the strong monotonicity of $G$ and $G(H(y^{\star}),y^{\star})=0$
we have

$\displaystyle\|{\hat{y}}_{k}-\beta_{k}G(H(y_{k}),y_{k})\|^{2}$
$\displaystyle=\|{\hat{y}}_{k}\|^{2}+\beta_{k}^{2}\|G(H(y_{k}),y_{k})-G(H(y^{\star}),y^{\star})\|^{2}-2\beta_{k}(y_{k}-y^{\star})^{T}G(x_{k},y_{k})$
$\displaystyle\leq(1-2\mu_{G}\beta_{k}+L^{2}(L+1)^{2}\beta_{k}^{2})\|{\hat{y}}_{k}\|^{2}.$
(3.32)

Next, using (3.32) and since $\beta_{k}\leq\mu_{G}/(4L^{2}(1+L)^{2})$ we
obtain

$\displaystyle
2\beta_{k}\big{(}{\hat{y}}_{k}-\beta_{k}G(H(y_{k}),y_{k})\big{)}^{T}\big{(}G(H(y_{k}),y_{k})-G(x_{k},y_{k})-\Delta
g_{k}\big{)}$ $\displaystyle\leq
2\beta_{k}\|{\hat{y}}_{k}-\beta_{k}G(H(y_{k}),y_{k})\|\big{(}L\|{\hat{x}}_{k}\|+\|\Delta
g_{k}\|\big{)}\leq
2\beta_{k}\|{\hat{y}}_{k}\|\big{(}L\|{\hat{x}}_{k}\|+\|\Delta g_{k}\|\big{)}$
$\displaystyle\leq\frac{\mu_{G}}{4}\beta_{k}\|{\hat{y}}_{k}\|^{2}+\frac{4}{\mu_{G}}\beta_{k}\big{(}L\|{\hat{x}}_{k}\|+\|\Delta
g_{k}\|\big{)}^{2}$ (3.33)

Substituting (3.32) and (3.33) into (3.31) and using (3.29) we obtain (3.30),
i.e.,

$\displaystyle\|{\hat{y}}_{k+1}\|^{2}$
$\displaystyle\leq(1-2\mu_{G}\beta_{k}+L^{2}(L+1)^{2}\beta_{k}^{2})\|{\hat{y}}_{k}\|^{2}+2\beta_{k}^{2}\|\Delta
g_{k}\|^{2}+2L^{2}\beta_{k}^{2}\|{\hat{x}}_{k}\|^{2}$
$\displaystyle\quad+\frac{\mu_{G}}{4}\beta_{k}\|{\hat{y}}_{k}\|^{2}+\frac{4}{\mu_{G}}\beta_{k}\big{(}L\|{\hat{x}}_{k}\|+\|\Delta
g_{k}\|\big{)}^{2}$
$\displaystyle\leq(1-\mu_{G}\beta_{k})\|{\hat{y}}_{k}\|^{2}-\big{(}\frac{3\mu_{G}}{4}-L^{2}(L+1)^{2}\beta_{k}\big{)}\beta_{k}\|{\hat{y}}_{k}\|^{2}+2\Big{(}\frac{4}{\mu_{G}}+\beta_{k}\Big{)}\beta_{k}\big{(}L^{2}\|{\hat{x}}_{k}\|^{2}+\|\Delta
g_{k}\|^{2}\big{)}$
$\displaystyle\leq(1-\mu_{G}\beta_{k})\|{\hat{y}}_{k}\|^{2}-\frac{\mu_{G}}{2}\beta_{k}\|{\hat{y}}_{k}\|^{2}+2\Big{(}\frac{2}{\mu_{G}}+\beta_{k}\Big{)}\beta_{k}\big{(}L^{2}\|{\hat{x}}_{k}\|^{2}+\|\Delta
g_{k}\|^{2}\big{)},$

where the last inequality we use $\beta_{k}\leq\mu_{G}/(4L^{2}(1+L)^{2})$. ∎

###  3.2 Finite-Time Convergence Rates

Leveraging the preliminary results in the previous section, we are ready to
present the main result of this paper. In Theorem 1 below, we will show that
the iterates generated by the proposed method in (2.7) will converge to the
desired solution $(x^{\star},y^{\star})$ in expectation at a rate ${\cal
O}(1/k)$. Our analysis will be based on the following candidate Lyapunov
function

$\displaystyle V_{k}=\|\Delta f_{k}\|^{2}+\|\Delta
g_{k}\|^{2}+\|{\hat{x}}_{k}\|^{2}+\|{\hat{y}}_{k}\|^{2}.$ (3.34)

We will show that each term on the right-hand side of (3.34) will converge to
zero at the same rate. Therefore, we do not need to use the coupling Lyapunov
function as studied in [24, 9].

######  Theorem 1.

Let $\\{f_{k},g_{k},x_{k},y_{k}\\}$ be generated by (2.7), and
$\alpha_{k},\beta_{k}$ satisfy (3.15) with $C_{\alpha},C_{\beta},h$ satisfy

$\displaystyle
C_{\alpha}\leq\min\left\\{\frac{\mu_{G}}{288L^{2}(1+L)^{2}},\;\frac{1}{64L^{2}(1+L)^{4}}\right\\},\quad
C_{\beta}=\frac{2}{\mu_{G}},\quad
h\geq\max\left\\{\frac{8(1+L)^{4}}{\mu_{G}^{3}\alpha_{0}},\;\frac{8L^{2}(1+L)^{2}}{\mu_{G}^{2}}\right\\}.$
(3.35)

Then we obtain

$\displaystyle\mathbb{E}[V_{k+1}]\leq\frac{h^{2}\mathbb{E}[V_{0}]}{(k+h+1)^{2}}+\frac{h^{2}(\Gamma_{11}+\Gamma_{22})(k+1)}{8(k+1+h)^{2}}\cdot$
(3.36)

###### Proof.

It is straightforward to verify that under the condition in (3.35)
$\alpha_{0}$ and $\beta_{0}$ satisfy (3.21) and (3.29), respectively. Thus,
taking the conditional expectation w.r.t ${\cal Q}_{k}$ on both sides of Eqs.
(3.17), (3.19), (3.22), and (3.30), and using the conditions of step sizes in
(3.15), i.e., $\lambda_{k}=\gamma_{k}$, we obtain

$\displaystyle\mathbb{E}[V_{k+1}\mid{\cal Q}_{k}]=\mathbb{E}[\|\Delta
f_{k+1}\|^{2}+\|\Delta
g_{k+1}\|^{2}+\|{\hat{x}}_{k+1}\|^{2}+\|{\hat{y}}_{k+1}\|^{2}]$
$\displaystyle\leq(1-\lambda_{k})\|\Delta
f_{k}\|^{2}-\Big{(}\frac{1}{4}-\frac{8L^{2}\alpha_{k}^{2}}{\lambda_{k}^{2}}\Big{)}\lambda_{k}\|\Delta
f_{k}\|^{2}$
$\displaystyle\quad+2\lambda_{k}^{2}\mathbb{E}\big{[}\|\epsilon_{k}\|^{2}\mid{\cal
Q}_{k}\big{]}+\frac{8L^{2}}{\lambda_{k}}\Big{(}\alpha_{k}^{2}\|\Delta
g_{k}\|^{2}+2L^{2}\alpha_{k}^{2}\|{\hat{x}}_{k}\|^{2}+(1+L)^{4}\beta_{k}^{2}\|{\hat{y}}_{k}\|^{2}\Big{)}$
$\displaystyle\quad+(1-\gamma_{k})\|\Delta
g_{k}\|^{2}-\Big{(}\frac{1}{4}-\frac{8L^{2}\alpha_{k}^{2}}{\gamma_{k}^{2}}\Big{)}\gamma_{k}\|\Delta
g_{k}\|^{2}$
$\displaystyle\quad+2\mathbb{E}\big{[}\gamma_{k}^{2}\|\psi_{k}\|^{2}\mid{\cal
Q}_{k}\big{]}+\frac{8L^{2}}{\gamma_{k}}\Big{(}\alpha_{k}^{2}\|\Delta
f_{k}\|^{2}+2L^{2}\alpha_{k}^{2}\|{\hat{x}}_{k}\|^{2}+(1+L)^{4}\beta_{k}^{2}\|{\hat{y}}_{k}\|^{2}\Big{)}$
$\displaystyle\quad+\left(1-\mu_{F}\alpha_{k}\right)\|{\hat{x}}_{k}\|^{2}-\left[\frac{\mu_{F}}{2}\alpha_{k}-\Big{(}3L+\frac{(1+L)^{4}}{\mu_{G}}\Big{)}\beta_{k}-\Big{(}3L^{4}+2L^{2}\Big{)}\alpha_{k}^{2}\right]\|{\hat{x}}_{k}\|^{2}$
$\displaystyle\quad+\frac{\mu_{G}}{4}\beta_{k}\|{\hat{y}}_{k}\|^{2}+(1+L)^{4}\big{(}3L^{2}\beta_{k}+\alpha_{k}\big{)}\beta_{k}\|{\hat{y}}_{k}\|^{2}$
$\displaystyle\quad+\Big{[}\frac{2}{\mu_{F}}\alpha_{k}+4\alpha_{k}^{2}\Big{]}\|\Delta
f_{k}\|^{2}+\Big{[}\frac{\alpha_{k}}{L}+(3L^{2}+1)\alpha_{k}^{2}\Big{]}\|\Delta
g_{k}\|^{2}$
$\displaystyle\quad+(1-\mu_{G}\beta_{k})\|{\hat{y}}_{k}\|^{2}-\frac{\mu_{G}}{2}\beta_{k}\|{\hat{y}}_{k}\|^{2}+2\Big{(}\frac{4}{\mu_{G}}+\beta_{k}\Big{)}\beta_{k}\big{(}L^{2}\|{\hat{x}}_{k}\|^{2}+\|\Delta
g_{k}\|^{2}\big{)}$ $\displaystyle=(1-\lambda_{k})\|\Delta
f_{k}\|^{2}+(1-\lambda_{k})\|\Delta
g_{k}\|^{2}+\left(1-\mu_{F}\alpha_{k}\right)\|{\hat{x}}_{k}\|^{2}+(1-\mu_{G}\beta_{k})\|{\hat{y}}_{k}\|^{2}$
$\displaystyle\quad+2\lambda_{k}^{2}\mathbb{E}\big{[}\|\epsilon_{k}\|^{2}+\|\psi_{k}\|^{2}\mid{\cal
Q}_{k}\big{]}-\Big{[}\frac{1}{4}-\frac{16L^{2}\alpha_{k}^{2}}{\lambda_{k}^{2}}-\frac{2\alpha_{k}}{\mu_{F}\lambda_{k}}-\frac{4\alpha_{k}^{2}}{\lambda_{k}}\Big{]}\lambda_{k}\|\Delta
f_{k}\|^{2}$
$\displaystyle\quad-\Big{[}\frac{1}{4}-\frac{16L^{2}\alpha_{k}^{2}}{\lambda_{k}^{2}}-\frac{\alpha_{k}}{L\lambda_{k}}-\frac{3(L^{2}+1)\alpha_{k}^{2}}{\lambda_{k}}-\frac{8\beta_{k}}{\mu_{G}\lambda_{k}}-\frac{2\beta_{k}^{2}}{\lambda_{k}}\Big{]}\lambda_{k}\|\Delta
g_{k}\|^{2}$
$\displaystyle\quad-\left[\frac{\mu_{F}}{2}\alpha_{k}-\Big{(}3L+\frac{(1+L)^{4}}{\mu_{G}}\Big{)}\beta_{k}-\Big{(}3L^{4}+2L^{2}\Big{)}\alpha_{k}^{2}+\frac{16L^{4}\alpha_{k}^{2}}{\lambda_{k}}+\frac{4L^{2}\beta_{k}}{\mu_{G}}+2L^{2}\beta_{k}^{2}\right]\|{\hat{x}}_{k}\|^{2}$
$\displaystyle\quad-\Big{[}\frac{\mu_{G}}{4}-\frac{16L^{2}(1+L)^{4}\beta_{k}}{\lambda_{k}}-(1+L)^{4}\big{(}3L^{2}\beta_{k}+\alpha_{k}\big{)}\Big{]}\beta_{k}\|{\hat{y}}_{k}\|^{2}$
$\displaystyle\leq(1-\lambda_{k})\|\Delta f_{k}\|^{2}+(1-\lambda_{k})\|\Delta
g_{k}\|^{2}+\left(1-\mu_{F}\alpha_{k}\right)\|{\hat{x}}_{k}\|^{2}+(1-\mu_{G}\beta_{k})\|{\hat{y}}_{k}\|^{2}$
$\displaystyle\quad+2\lambda_{k}^{2}\mathbb{E}\big{[}\|\epsilon_{k}\|^{2}+\|\psi_{k}\|^{2}\mid{\cal
Q}_{k}\big{]}-\Big{[}\frac{1}{4}-\frac{16L^{2}\alpha_{k}^{2}}{\lambda_{k}^{2}}-\frac{2\alpha_{k}}{\mu_{F}\lambda_{k}}-\frac{4\alpha_{k}^{2}}{\lambda_{k}}\Big{]}\lambda_{k}\|\Delta
f_{k}\|^{2}$
$\displaystyle\quad-\Big{[}\frac{1}{4}-\frac{16L^{2}\alpha_{k}^{2}}{\lambda_{k}^{2}}-\frac{9\alpha_{k}}{\mu_{G}\lambda_{k}}-\frac{3(L+1)^{2}\alpha_{k}^{2}}{\lambda_{k}}\Big{]}\lambda_{k}\|\Delta
g_{k}\|^{2}$
$\displaystyle\quad-\Big{[}\frac{\mu_{F}}{2}-\frac{2(1+L)^{4}\beta_{k}}{\mu_{G}\alpha_{k}}-\frac{16L^{4}\alpha_{k}}{\lambda_{k}}-4L^{2}(L+1)^{2}\alpha_{k}\Big{]}\alpha_{k}\|{\hat{x}}_{k}\|^{2}$
$\displaystyle\quad-\Big{[}\frac{\mu_{G}}{4}-3(L+1)^{6}\alpha_{k}-\frac{16L^{2}(1+L)^{4}\alpha_{k}}{\lambda_{k}}\Big{]}\beta_{k}\|{\hat{y}}_{k}\|^{2},$
(3.37)

where the last inequality we use $\beta_{k}\leq\alpha_{k}$ and
$\mu_{G}\leq\mu_{F}\leq L$. Using the conditions in (3.35) we can show that
each term in the squared brackets on the right-hand side of (3.37) is
nonnegative. Therefore, since $\beta_{k}\leq\alpha_{k}\leq\lambda_{k}$ and by
taking the expectation on both sides of (3.37) we have

$\displaystyle\mathbb{E}[V_{k+1}]$
$\displaystyle\leq(1-\mu_{G}\beta_{k})\mathbb{E}[V_{k}]+2\lambda_{k}^{2}\big{(}\mathbb{E}[\|\epsilon_{k}\|]+\mathbb{E}[\|\psi_{k}\|^{2}]\big{)}$
$\displaystyle=\left(1-\frac{2}{k+h+1}\right)\mathbb{E}[V_{k}]+\frac{\Gamma_{11}+\Gamma_{22}}{8(k+1)^{2}}=\frac{k+h-1}{k+h+1}\mathbb{E}[V_{k}]+\frac{\Gamma_{11}+\Gamma_{22}}{8(k+1)^{2}},$

where the second inequality is due to Assumption 3. Multiplying both sides of
the preceding relation by $(k+h+1)^{2}$ and since $h>1$ we obtain

$\displaystyle(k+h+1)^{2}\mathbb{E}[V_{k+1}]$
$\displaystyle\leq(k+h+1)(k+h-1)\mathbb{E}[V_{k}]+\frac{(\Gamma_{11}+\Gamma_{22})(k+h+1)^{2}}{8(k+1)^{2}}$
$\displaystyle\leq(k+h)^{2}\mathbb{E}[V_{0}]+\frac{h^{2}(\Gamma_{11}+\Gamma_{22})}{8}$
$\displaystyle\leq
h^{2}\mathbb{E}[V_{0}]+\frac{h^{2}(\Gamma_{11}+\Gamma_{22})}{8}(k+1),$

which when diving both sides by $(k+1+h)^{2}$ yields (3.36).

######  Remark 4.

In [10], the authors show that the updates in (1.2) can achieve an ${\cal
O}(1/k)$ convergence rate under a strong assumption on $H$, i.e., $H$ is
differentiable and satisfies

$\displaystyle\|H(y_{1})-H(y_{2})-\nabla H(y_{2})^{T}(y_{1}-y_{2})\|\leq
L\|y_{1}-y_{2}\|^{\zeta}\quad\text{for }\zeta\in[1.5,2].$

This assumption basically implies that the impact of the noise in $y$ to $x$
update is negligible, i.e., on the high order of step sizes. Indeed, under
this assumption the analysis in [9] (with some simple modification) also
implies that the two-time-scale SA in (1.2) has a convergence rate ${\cal
O}(1/k)$. While this assumption satisfies in some settings, e.g., see examples
presented in [10], it may not hold in general (e.g., $H(y)=|y|$ where
$y\in\mathbb{R}$). On the other hand, in this paper we propose a new approach
to achieve an ${\cal O}(1/k)$ convergence under a more general condition on
$H$, i.e., $H$ is Lipschitz continuous. Finally, we note that when $H$ is
linear as presented in Remark 2, the high-order smooth condition above is
satisfied. Thus, one can utilize the techniques in [10, 9] to show an ${\cal
O}(1/k)$ convergence rate of (1.2). This can provide an alternative approach
to achieve the same convergence rate as studied in [7, 8] for the linear
setting.

∎

##  4 Motivating applications

In this section, we present a few existing applications of two-time-scale SA
in (2.7), which serve as the main motivation of this paper. Our improved
convergence rate in Theorem 1 suggests that one can apply the idea of the
proposed method in (2.7) to develop new algorithms for solving problems in
these applications.

#### Reinforcement Learning.

Two-time-scale SA has been used extensively to model reinforcement learning
methods, for example, gradient temporal difference (TD) learning and actor-
critic methods [32, 33, 34, 35, 36, 37, 38, 39, 40]. For example, the gradient
TD learning for solving the policy evaluation problem under nonlinear function
approximations studied in [34] can be viewed as a variant of (1.2). In this
problem, we want to estimate the cumulative rewards $V$ of a stationary policy
using function approximations $V_{y}$, that is, our goal is to find $y$ so
that $V_{y}$ is as close as possible to the true value $V$. Here, $V_{y}$ can
be represented by a neural network where $y$ is the weight vector of the
network. Let $\zeta$ be the environmental sate, $\gamma$ be the discount
factor, $\phi(\zeta)=\nabla V_{y}(\zeta)$ be the feature vector of state
$\zeta$, and $r$ be the reward return by the environment. Given a sequence of
samples $\\{\zeta_{k},r_{k}\\}$, one version of GTD are given as

$\displaystyle x_{k+1}$
$\displaystyle=x_{k}+\alpha_{k}(\delta_{k}-\phi(\zeta_{k})^{T}x_{k})\phi(\zeta_{k})$
$\displaystyle y_{k+1}$
$\displaystyle=y_{k}+\beta_{k}\Big{[}\left(\phi(\zeta_{k})-\gamma\phi(\zeta_{k+1})\right)\phi(\zeta_{k})^{T}x_{k}-h_{k}\Big{]},$

where $\delta_{k}$ and $h_{k}$ are defined as

$\displaystyle\delta_{k}$ $\displaystyle=r_{k}+\gamma
V_{y_{k}}(\zeta_{k+1})-V_{y_{k}}(\zeta_{k})$ $\displaystyle h_{k}$
$\displaystyle=(\delta_{k}-\phi(\zeta_{k})^{T}x_{k})\nabla^{2}V_{y_{k}}(\zeta_{k})x_{k},$

which is clearly a variant of (1.2) under some proper choice of $F$ and $G$.
It has been observed that the GTD method is more stable and performs better
when compared to the single-time-scale counterpart (TD learning) under off-
policy learning and nonlinear function approximations [13]. The proposed
method in this paper suggests that it may be better to use the Ruppert-Polyak
averaging steps on the samples to reduce the impacts of noise before updating
the main iterates.

#### Game Theory.

Two-time-scale SA can also be used to represent different learning methods in
game theory, in particular, in zero-sum settings. For example, consider a two-
player Markov game characterized by ${\cal M}=({\cal S},{\cal A},{\cal
B},{\cal P},\gamma,r)$. Here, ${\cal S}$ is the finite state space, ${\cal A}$
and ${\cal B}$ are the finite action spaces of the two players,
$\gamma\in(0,1)$ is the discount factor, and $r:{\cal S}\times{\cal
A}\times{\cal B}\rightarrow[0,1]$ is the reward function. Let $\Delta_{{\cal
F}}$ denote the probability simplex over a set ${\cal F}$, and ${\cal P}:{\cal
S}\times{\cal A}\times{\cal B}\rightarrow\Delta_{{\cal S}}$ be the transition
probability kernel, with ${\cal P}(s^{\prime}\mid s,a,b)$ specifying the
probability of the game transitioning from state $s$ to $s^{\prime}$ when the
first player selects action $a\in{\cal A}$ and the second player selects
$b\in{\cal B}$. The policies of the two players are denoted by
$\pi\in\Delta_{{\cal A}}^{{\cal S}}$ and $\phi\in\Delta_{{\cal B}}^{{\cal
S}}$, with $\pi(a\mid s)$, $\phi(b\mid s)$ denoting the probability of
selecting action $a$, $b$ in state $s$ according to $\pi$, $\phi$. Given a
policy pair $(\pi,\phi)$, we measure its performance in state $s\in{\cal S}$
by the value function

$\displaystyle V^{\pi,\phi}(s)=\mathbb{E}_{a_{k}\sim\pi(\cdot\mid
s_{k}),b_{k}\sim\phi(\cdot\mid s_{k}),s_{k+1}\sim{\cal P}(\cdot\mid
s_{k},a_{k},b_{k})}\Big{[}\sum\nolimits_{k=0}^{\infty}\gamma^{k}r\left(s_{k},a_{k},b_{k}\right)\mid
s_{0}=s\Big{]}.$

Under a fixed initial distribution $\rho\in\Delta_{{\cal S}}$, we define by
$J(\pi,\phi)\triangleq\mathbb{E}_{s_{0}\sim\rho}[V^{\pi,\phi}(s_{0})]$ the
discounted cumulative reward under $(\pi,\phi)$. The goal is to find a Nash
equilibrium $(\pi^{\star},\phi^{\star})$, which can be shown to be the
solution of the following minimax optimization problem

$\displaystyle\max_{\pi\in\Delta_{{\cal A}}^{{\cal
S}}}\min_{\phi\in\Delta_{{\cal B}}^{{\cal
S}}}J(\pi,\phi)=\min_{\phi\in\Delta_{{\cal B}}^{{\cal
S}}}\max_{\pi\in\Delta_{{\cal A}}^{{\cal
S}}}J(\pi,\phi)=J(\pi^{\star},\phi^{\star}).$

When the transition probability kernel ${\cal P}$ is known, one can apply
dynamic programming approach to solve this problem. On the other hand, when
${\cal P}$ is unknown one can use reinforcement learning methods to find the
Nash equilibrium. For example, we can apply the decentralized Q-learning
method recently studied in [41], which is a variant of two-time-scale SA. This
method maintains two variables, namely, local value functions and local
q-function for each player, that are updated at different time scales by using
different step sizes. In particular, let $v_{1}^{\pi}(s),q_{1}^{\pi}(s,a)$ be
the local value function and state-action value function at player $1$,
respectively. Similar notation is defined for player $2$. Player $1$ will
update its local functions as (similar to player $2$)

$\displaystyle q_{i}^{k+1}(s_{k},a_{k})$
$\displaystyle=q_{i}^{k+1}(s_{k},a_{k})+\alpha_{k}\left(r(s_{k},a_{k},b_{k})+\gamma
v_{i}^{k}(s_{k}^{\prime})-q_{i}^{k+1}(s_{k},a_{k})\right)$ $\displaystyle
v_{i}^{k+1}(s_{k})$
$\displaystyle=v_{i}^{k}(s_{k})+\beta_{k}\big{(}\big{[}\bar{\pi}_{i}^{k}\big{]}^{T}q_{i}^{k+1}(s_{k},\cdot)-v_{i}^{k},(s_{k})\big{)}$

where $\beta_{k}\leq\alpha_{k}$ and $\bar{\pi}_{i}^{k}$ is the best response
policy w.r.t $q_{i}^{k+1}$. While asymptotic convergence of this method has
been studied in [41] based on the ODE method, its finite-time complexity is
completely unknown. Thus, an interesting question is to study whether one can
apply the technique studied in this paper to derive the convergence rate of
the decentralized Q-learning method. Another question is to see whether one
can improve the existing convergence rates of policy gradient descent-ascent
algorithm for zero-sum Markov games studied in [42], where the authors show
that this algorithm converges at a rate ${\cal O}(1/k^{1/3})$. Although this
result is not comparable to the one studied in Theorem 1, the objective
function $J(\pi,\phi)$ is nonconvex and nonconcave, therefore, one cannot
immediately apply our analysis to the zero-sum Markov game setting.

#### Minimax Optimization.

We consider the following minimax optimization problems

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

where $f:\mathbb{R}^{m}\times\mathbb{R}^{n}\rightarrow\mathbb{R}$ is a
(non)convex function w.r.t $x$ for a fixed $y$ and (non)concave w.r.t $y$ for
a fixed $x$. The minimax problem has broad applications in different areas
including game theory [43], training generative adversarial networks (GANs)
[44, 45], adversarial and robust machine learning [46, 47], resource
allocation over networks [48], and distributed optimization [49, 50]; to name
just a few. For solving this problem, (alternating) gradient descent-ascent is
the most popular method used in the existing literature; see for example [51,
52, 53, 54, 55, 56, 57] and the references therein. This method iteratively
updates the estimates $x_{k},y_{k}$ of the desired solution as

$\displaystyle
x_{k+1}=x_{k}-\alpha_{k}\nabla_{x}f(x_{k},y_{k};\xi_{k})\quad\text{and}\quad
y_{k+1}=y_{k}-\beta_{k}\nabla_{y}f(x_{k},y_{k};\psi_{k}),$ (4.38)

where $\alpha_{k},\beta_{k}$ are chosen differently to guarantee the
convergence of this method. It is obvious that gradient descent ascent is a
variant of the two-time-scale SA.

#### Other Applications.

Finally, two-time-scale SA has been used in distributed optimization to
address the issues of communication constraints [58, 59, 60, 61] and in
distributed control to handle clustered networks [62, 63, 64, 65, 66, 67, 68].
For example, the best convergence rate of distributed consensus-based policy
gradient methods under random quantization is ${\cal O}(1/\sqrt{k})$ when the
objective function is strongly convex [61], a variant of strong monotonicity
in Assumption 2. Thus, it will be interesting to see whether we can apply the
proposed technique in this paper to improve the convergence rate of
distributed optimization method under random quantization.

##  5 Concluding Remarks

In this paper, by leveraging the Ruppert-Polyak averaging techniques we
propose to develop a new variant of the classic two-time-scale SA that
achieves a finite-time convergence rate ${\cal O}(1/k)$ in mean-squared
regime. Our result significantly improves the existing convergence rate of
two-time-scale SA, which is ${\cal O}(1/k^{2/3})$. The key idea of our
approach is to apply the averaging technique to estimate the operators from
their samples before updating the main iterates. This can help to reduce the
impact of the sampling noise on the updates of the iterates, therefore,
improving its convergence towards the desired solutions. An interesting future
direction of this paper is to apply the proposed method to improve the
performance of existing control, optimization, and learning algorithms that
can be represented as variants of two-time-scale SA.

## References

  * [1] V. S. Borkar, _Stochastic Approximation: A Dynamical Systems Viewpoint_. Cambridge University Press, 2008.
  * [2] H. Robbins and S. Monro, “A stochastic approximation method,” _The Annals of Mathematical Statistics_ , vol. 22, no. 3, pp. 400–407, 1951.
  * [3] A. Benveniste, M. Métivier, and P. Priouret, _Adaptive algorithms and stochastic approximations_. Springer Science & Business Media, 2012, vol. 22.
  * [4] V. R. Konda and J. N. Tsitsiklis, “Convergence rate of linear two-time-scale stochastic approximation,” _The Annals of Applied Probability_ , vol. 14, no. 2, pp. 796–819, 2004.
  * [5] A. Mokkadem and M. Pelletier, “Convergence rate and averaging of nonlinear two-time-scale stochastic approximation algorithms,” _The Annals of Applied Probability_ , vol. 16, no. 3, pp. 1671–1702, 2006.
  * [6] J. Hu, V. Doshi, and D. Y. Eun, “Central limit theorem for two-timescale stochastic approximation with markovian noise: Theory and applications,” _To appear in AISTAT Artificial Intelligence and Statistics (AISTATS). arXiv preprint arXiv:2401.09339_ , 2024.
  * [7] M. Kaledin, E. Moulines, A. Naumov, V. Tadic, and H.-T. Wai, “Finite time analysis of linear two-timescale stochastic approximation with Markovian noise,” in _Proceedings of Thirty Third Conference on Learning Theory_ , vol. 125, 2020, pp. 2144–2203.
  * [8] S. U. Haque, S. Khodadadian, and S. T. Maguluri, “Tight finite time bounds of two-time-scale linear stochastic approximation with markovian noise,” _arXiv preprint arXiv:2401.00364_ , 2023.
  * [9] T. T. Doan, “Nonlinear two-time-scale stochastic approximation: Convergence and finite-time performance,” _IEEE Transactions on Automatic Control_ , vol. 68, no. 8, pp. 4695–4705, 2023.
  * [10] Y. Han, X. Li, and Z. Zhang, “Finite-time decoupled convergence in nonlinear two-time-scale stochastic approximation,” _arXiv preprint arXiv:2401.03893_ , 2024.
  * [11] D. Bertsekas and J. Tsitsiklis, _Neuro-Dynamic Programming_ , 2nd ed. Athena Scientific, Belmont, MA, 1999.
  * [12] T. Hastie, R. Tibshirani, and J. Friedman, _The Elements of Statistical Learning : Data Mining, Inference, and Prediction_. Springer, 2009.
  * [13] R. S. Sutton and A. G. Barto, _Reinforcement Learning: An Introduction_ , 2nd ed. MIT Press, Cambridge, MA, 2018\.
  * [14] G. Lan, _Lectures on Optimization Methods for Machine Learning_. Springer-Nature, 2020.
  * [15] L. Bottou, F. Curtis, and J. Nocedal, “Optimization methods for large-scale machine learning,” _SIAM Review_ , vol. 60, no. 2, pp. 223–311, 2018.
  * [16] J. Bhandari, D. Russo, and R. Singal, “A finite time analysis of temporal difference learning with linear function approximation,” in _COLT_ , 2018\.
  * [17] B. Karimi, B. Miasojedow, E. Moulines, and H.-T. Wai, “Non-asymptotic analysis of biased stochastic approximation scheme,” in _Conference on Learning Theory, COLT 2019, 25-28 June 2019, Phoenix, AZ, USA_ , 2019, pp. 1944–1974.
  * [18] R. Srikant and L. Ying, “Finite-time error bounds for linear stochastic approximation and TD learning,” in _COLT_ , 2019.
  * [19] B. Hu and U. Syed, “Characterizing the exact behaviors of temporal difference learning algorithms using markov jump linear system theory,” in _Advances in Neural Information Processing Systems 32_ , 2019.
  * [20] S. Chen, A. Devraj, A. Busic, and S. Meyn, “Explicit mean-square error bounds for monte-carlo and linear stochastic approximation,” ser. Proceedings of Machine Learning Research, vol. 108, 26–28 Aug 2020, pp. 4173–4183.
  * [21] Z. Chen, S. Zhang, T. T. Doan, S. T. Maguluri, and J.-P. Clarke, “Performance of Q-learning with Linear Function Approximation: Stability and Finite-Time Analysis,” Available at: https://arxiv.org/abs/1905.11425, 2019.
  * [22] G. Dalal, G. Thoppe, B. Szörényi, and S. Mannor, “Finite sample analysis of two-timescale stochastic approximation with applications to reinforcement learning,” in _COLT_ , 2018.
  * [23] T. T. Doan and J. Romberg, “Linear two-time-scale stochastic approximation a finite-time analysis,” in _2019 57th Annual Allerton Conference on Communication, Control, and Computing (Allerton)_ , 2019, pp. 399–406.
  * [24] H. Gupta, R. Srikant, and L. Ying, “Finite-time performance bounds and adaptive learning rate selection for two time-scale reinforcement learning,” in _Advances in Neural Information Processing Systems_ , 2019.
  * [25] T. T. Doan, “Finite-time analysis and restarting scheme for linear two-time-scale stochastic approximation,” _SIAM Journal on Control and Optimization_ , vol. 59, no. 4, pp. 2798–2819, 2021.
  * [26] G. Dalal, B. Szorenyi, and G. Thoppe, “A tale of two-timescale reinforcement learning with the tightest finite-time bound,” _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 34, no. 04, pp. 3701–3708, Apr. 2020\.
  * [27] S. Zeng, T. T. Doan, and J. Romberg, “A two-time-scale stochastic optimization framework with applications in control and reinforcement learning,” _arXiv preprint arXiv:2109.14756_ , 2021.
  * [28] P. Kokotović, H. K. Khalil, and J. O’Reilly, _Singular Perturbation Methods in Control: Analysis and Design_. Society for Industrial and Applied Mathematics, 1999.
  * [29] D. Ruppert, “Efficient estimations from a slowly convergent robbins-monro process,” Cornell University Operations Research and Industrial Engineering, Tech. Rep., 1988.
  * [30] B. T. Polyak and A. B. Juditsky, “Acceleration of stochastic approximation by averaging,” _SIAM journal on control and optimization_ , vol. 30, no. 4, pp. 838–855, 1992.
  * [31] S. Chen, A. Devraj, A. Bernstein, and S. Meyn, “Accelerating optimization and reinforcement learning with quasi stochastic approximation,” in _2021 American Control Conference (ACC)_. IEEE, 2021, pp. 1965–1972.
  * [32] R. Sutton, H. R. Maei, and C. Szepesvári, “A convergent o(n) temporal-difference algorithm for off-policy learning with linear function approximation,” in _Advances in Neural Information Processing Systems 21_ , 2009.
  * [33] R. Sutton, H. R. Maei, D. Precup, S. Bhatnagar, D. Silver, C. Szepesvári, and E. Wiewiora, “Fast gradient-descent methods for temporal-difference learning with linear function approximation,” in _Proceedings of the 26th International Conference On Machine Learning, ICML_ , vol. 382, 01 2009.
  * [34] H. R. Maei, C. Szepesvári, S. Bhatnagar, D. Precup, D. Silver, and R. S. Sutton, “Convergent temporal-difference learning with arbitrary smooth function approximation,” in _Proceedings of the 22nd International Conference on Neural Information Processing Systems_ , 2009, p. 1204–1212.
  * [35] V. R. Konda and J. N. Tsitsiklis, “On actor-critic algorithms,” _SIAM J. Control Optim._ , vol. 42, no. 4, 2003.
  * [36] T. Xu, Z. Wang, and Y. Liang, “Improving sample complexity bounds for actor-critic algorithms,” Available at: https://arxiv.org/abs/2010.00355, 2020.
  * [37] Y. Wu, W. Zhang, P. Xu, and Q. Gu, “A finite time analysis of two time-scale actor critic methods,” Available at: https://arxiv.org/abs/2005.01350, 2020.
  * [38] M. Hong, H.-T. Wai, Z. Wang, and Z. Yang, “A two-timescale framework for bilevel optimization: Complexity analysis and application to actor-critic,” Available at: https://arxiv.org/abs/2007.05170, 2020.
  * [39] S. Khodadadian, T. T. Doan, S. T. Maguluri, and J. Romberg, “Finite Sample Analysis of Two-Time-Scale Natural Actor-Critic Algorithm,” Available at: https://arxiv.org/abs/2101.10506, 2021.
  * [40] S. Zeng, T. T. Doan, and J. Romberg, “Finite-time complexity of online primal-dual natural actor-critic algorithm for constrained markov decision processes,” in _2022 IEEE 61st Conference on Decision and Control (CDC)_ , 2022, pp. 4028–4033.
  * [41] M. Sayin, K. Zhang, D. Leslie, T. Basar, and A. Ozdaglar, “Decentralized q-learning in zero-sum markov games,” _Advances in Neural Information Processing Systems_ , vol. 34, pp. 18 320–18 334, 2021.
  * [42] S. Zeng, T. Doan, and J. Romberg, “Regularized gradient descent ascent for two-player zero-sum markov games,” _Advances in Neural Information Processing Systems_ , vol. 35, pp. 34 546–34 558, 2022.
  * [43] T. Basar and G. J. Olsder, _Dynamic Noncooperative Game Theory, 2nd Edition_. Society for Industrial and Applied Mathematics, 1998.
  * [44] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio, “Generative adversarial networks,” _Commun. ACM_ , vol. 63, no. 11, p. 139–144, 2020.
  * [45] L. Mescheder, S. Nowozin, and A. Geiger, “The numerics of gans,” in _Proceedings of the 31st International Conference on Neural Information Processing Systems_ , ser. NIPS’17, 2017, p. 1823–1833.
  * [46] A. Kurakin, I. J. Goodfellow, and S. Bengio, “Adversarial machine learning at scale,” in _5th International Conference on Learning Representations, ICLR_ , 2017.
  * [47] Q. Qian, S. Zhu, J. Tang, R. Jin, B. Sun, and H. Li, “Robust optimization over multiple domains,” _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 33, pp. 4739–4746, Jul. 2019.
  * [48] Y.-F. Liu, Y.-H. Dai, and Z.-Q. Luo, “Max-min fairness linear transceiver design for a multi-user mimo interference channel,” _IEEE Transactions on Signal Processing_ , vol. 61, no. 9, pp. 2413–2423, 2013.
  * [49] G. Lan, S. Lee, and Y. Zhou, “Communication-efficient algorithms for decentralized and stochastic optimization,” _Mathematical Programming_ , vol. 180, pp. 237–284, 2020.
  * [50] T.-H. Chang, M. Hong, H.-T. Wai, X. Zhang, and S. Lu, “Distributed learning in the nonconvex world: From batch data to streaming and beyond,” _IEEE Signal Processing Magazine_ , vol. 37, no. 3, pp. 26–38, 2020.
  * [51] T. Lin, C. Jin, and M. Jordan, “On gradient descent ascent for nonconvex-concave minimax problems,” in _Proceedings of the 37th International Conference on Machine Learning_ , vol. 119. PMLR, 13–18 Jul 2020, pp. 6083–6093.
  * [52] S. Lu, I. Tsaknakis, M. Hong, and Y. Chen, “Hybrid block successive approximation for one-sided non-convex min-max problems: Algorithms and applications,” _IEEE Transactions on Signal Processing_ , vol. 68, pp. 3676–3691, 2020.
  * [53] J. Yang, N. Kiyavash, and N. He, “Global convergence and variance reduction for a class of nonconvex-nonconcave minimax problems,” in _Advances in Neural Information Processing Systems_ , vol. 33, 2020, pp. 1153–1165.
  * [54] Z. Xu, H.-L. Zhang, Y. Xu, and G. Lan, “A unified single-loop alternating gradient projection algorithm for nonconvex-concave and convex-nonconcave minimax problems,” _ArXiv_ , vol. abs/2006.02032, 2020.
  * [55] J. Zhang, P. Xiao, R. Sun, and Z. Luo, “A single-loop smoothed gradient descent-ascent algorithm for nonconvex-concave min-max problems,” in _Advances in Neural Information Processing Systems_ , vol. 33, 2020, pp. 7377–7389.
  * [56] T. Doan, “Convergence rates of two-time-scale gradient descent-ascent dynamics for solving nonconvex min-max problems,” in _Learning for Dynamics and Control Conference_. PMLR, 2022, pp. 192–206.
  * [57] A. Cherukuri, B. Gharesifard, and J. Cortés, “Saddle-point dynamics: Conditions for asymptotic stability of saddle points,” _SIAM Journal on Control and Optimization_ , vol. 55, no. 1, pp. 486–511, 2017.
  * [58] T. T. Doan, C. L. Beck, and R. Srikant, “On the convergence rate of distributed gradient methods for finite-sum optimization under communication delays,” _Proceedings ACM Meas. Anal. Comput. Syst._ , vol. 1, no. 2, pp. 37:1–37:27, 2017.
  * [59] T. T. Doan, S. T. Maguluri, and J. Romberg, “Convergence rates of distributed gradient methods under random quantization: A stochastic approximation approach,” _IEEE on Transactions on Automatic Control_ , 2020.
  * [60] A. Reisizadeh, A. Mokhtari, H. Hassani, and R. Pedarsani, “Quantized decentralized consensus optimization,” in _2018 IEEE Conference on Decision and Control (CDC)_ , 2018, pp. 5838–5843.
  * [61] M. M. Vasconcelos, T. T. Doan, and U. Mitra, “Improved convergence rate for a distributed two-time-scale gradient method under random quantization,” in _2021 60th IEEE Conference on Decision and Control (CDC)_ , 2021, pp. 3117–3122.
  * [62] D. Romeres, F. Dörfler, and F. Bullo, “Novel results on slow coherency in consensus and power networks,” in _Proc. of 2013 European Control Conference_ , 2013, pp. 742–747.
  * [63] J. Chow and P. Kokotovic, “Time scale modeling of sparse dynamic networks,” _IEEE Transactions on Automatic Control_ , vol. 30, no. 8, pp. 714–722, 1985\.
  * [64] E. Biyik and M. Arcak, “Area aggregation and time-scale modeling for sparse nonlinear networks,” _Systems and Control Letters_ , vol. 57, no. 2, pp. 142–149, 2008.
  * [65] A. M. Boker, C. Yuan, F. Wu, and A. Chakrabortty, “Aggregate control of clustered networks with inter-cluster time delays,” in _Proc. of 2016 American Control Conference_ , 2016, pp. 5340–5345.
  * [66] T. V. Pham, T. T. Doan, and D. H. Nguyen, “Distributed two-time-scale methods over clustered networks,” in _2021 American Control Conference (ACC)_. IEEE, 2021, pp. 4625–4630.
  * [67] A. Dutta, A. M. Boker, and T. T. Doan, “Convergence rates of distributed consensus over cluster networks: A two-time-scale approach,” in _2022 IEEE 61st Conference on Decision and Control (CDC)_. IEEE, 2022, pp. 7035–7040.
  * [68] A. Dutta, N. Masrourisaadat, and T. T. Doan, “Convergence rates of decentralized gradient dynamics over cluster networks: Multiple-time-scale lyapunov approach,” in _2022 IEEE 61st Conference on Decision and Control (CDC)_ , 2022, pp. 6497–6502.