#  MDPGT: Momentum-based Decentralized Policy Gradient Tracking

Zhanhong Jiang1, Xian Yeow Lee2, Sin Yong Tan2, Kai Liang Tan2,
Aditya Balu2, Young M. Lee1, Chinmay Hegde3, Soumik Sarkar2

###### Abstract

We propose a novel policy gradient method for multi-agent reinforcement
learning, which leverages two different variance-reduction techniques and does
not require large batches over iterations. Specifically, we propose a
momentum-based decentralized policy gradient tracking (MDPGT) where a new
momentum-based variance reduction technique is used to approximate the local
policy gradient surrogate with importance sampling, and an intermediate
parameter is adopted to track two consecutive policy gradient surrogates.
Moreover, MDPGT provably achieves best available sample complexity of
$\mathcal{O}(N^{-1}\epsilon^{-3})$ for converging to an $\epsilon$-stationary
point of the global average of $N$ local performance functions (possibly
nonconcave). This outperforms the state-of-the-art sample complexity in
decentralized model-free reinforcement learning and when initialized with a
single trajectory, the sample complexity matches those obtained by the
existing decentralized policy gradient methods. We further validate the
theoretical claim for the Gaussian policy function. When the required error
tolerance $\epsilon$ is small enough, MDPGT leads to a linear speed up, which
has been previously established in decentralized stochastic optimization, but
not for reinforcement leaning. Lastly, we provide empirical results on a
multi-agent reinforcement learning benchmark environment to support our
theoretical findings.

## Introduction

Multi-agent reinforcement learning (MARL) is an emerging topic which has been
explored both in theoretical (Nguyen et al. 2014; Zhang et al. 2018; Qu et al.
2019; Zhang et al. 2021c) and empirical settings (Helou, Kalathil, and Xie
2020; Mukherjee, Bai, and Chakrabortty 2020; Zhou et al. 2020). Several
appealing applications of MARL can be seen in (Zhang, Yang, and Başar 2019;
Nguyen, Nguyen, and Nahavandi 2020) and relevant references therein.

While MARL can primarily be cast into two different categories, i.e.,
cooperative (Li, Chen, and Chen 2020; Wang et al. 2020; Li et al. 2020) and
competitive (Chen et al. 2020), our focus is in the cooperative setting; see
(Wei et al. 2021) for details on the competitive setting. Cooperative MARL is
typically modeled as a networked multi-agent Markov decision process (MDP)
(Zhang, Yang, and Basar 2018; Chu, Chinchali, and Katti 2020; Zhang et al.
2018) in which the agents share a centralized reward function (Simões, Lau,
and Reis 2020; Ackermann et al. 2019). However, in practice, this is not
necessarily the case, and instead a more general yet challenging scenario is
that agents have heterogeneous reward functions. Inherently, the ultimate goal
in such a cooperative MARL setting is for agents to maximize the global
average of local long-term returns. To address this problem, various
algorithms have been proposed, including distributed-learning (Arslan and
Yüksel 2016; Nguyen and Mukhopadhyay 2017) and distributed actor-critic (Li et
al. 2020; Ryu, Shin, and Park 2020). More recent works have successfully
showed finite-sample analysis for decentralized batch MARL (Zhang et al.
2021c) and leveraged advances in analysis of descent-ascent algorithms (Lu et
al. 2021).

These preliminary attempts have facilitated the theoretical understanding of
cooperative MARL by showing explicit sample complexity bounds, which match
that of standard (vanilla) stochastic gradient descent (SGD). Additionally,
recent works (Huang et al. 2020; Xu, Gao, and Gu 2019) in centralized RL have
revealed that with simple _variance reduction_ techniques, this sample
complexity can be reduced to $\mathcal{O}(\epsilon^{-3})$ to reach an
$\epsilon$-stationary point (i.e., $\mathbb{E}[\|\nabla
J(\mathbf{x})\|]\leq\epsilon$, where $J$ is the return function and
$\mathbf{x}\in\mathbb{R}^{d}$ is the decision variable to be optimized), which
has been admitted as the best complexity in decentralized optimization (Das et
al. 2020; Karimireddy et al. 2020). However, no similar matching bounds have
yet been reported in the decentralized (cooperative MARL) setting. Hence, this
motivates the question:

> _Can we achieve a sample complexity of $\mathcal{O}(\epsilon^{-3})$ in
> decentralized MARL via variance reduction?_

In this paper, we answer this question affirmatively by proposing a variance-
reduced policy gradient tracking approach, MDPGT (Algorithm 1), and analyzing
it in Theorem 1. Additionally, we propose a variation (based on a different
initialization) that enables state-of-the-art (SOTA) sample complexity for
decentralized MARL (Zhang et al. 2021c; Lu et al. 2021). See Table 1 for SOTA
comparisons. Specifically:

  1. 1.

We propose MDPGT, in which we use a stochastic policy gradient surrogate, a
convex combination of the vanilla stochastic policy gradient and an importance
sampling-based stochastic recursive algorithm (SARAH) (Nguyen et al. 2017) for
the local gradient update. Instead of directly applying the stochastic policy
gradient surrogate in the parameter update, an intermediate parameter is
adopted to track the difference between two consecutive stochastic policy
gradient surrogates. For smooth nonconcave performance functions, we show that
MDPGT with the mini-batch initialization can converge to an
$\epsilon$-stationary point in $\mathcal{O}(N^{-1}\epsilon^{-3})$ gradient-
based updates which matches the best available known upper bounds (Huang et
al. 2020).

  2. 2.

We modify the initialization of the proposed algorithm MDPGT by using a single
trajectory instead of a mini-batch of trajectories. Surprisingly, we find that
only one trajectory results in a larger sampling complexity
$\mathcal{O}(N^{-1}\epsilon^{-4})$, which, however, is the same as obtained by
the SOTA (Zhang et al. 2021c; Lu et al. 2021) with a linear speed up when
$\epsilon$ is sufficiently small. Additionally, our algorithm shows that when
updating the policy parameter in MDPGT, the mini-batch size is
$\mathcal{O}(1)$ instead of being $\epsilon$-related (Xu, Gao, and Gu 2019; Qu
et al. 2019), which can significantly improve practical efficiency.

  3. 3.

To facilitate the theoretical understanding of MDPGT, we leverage a benchmark
gridworld environment for numerical simulation and compare our proposed
algorithm to a baseline decentralized policy gradient (DPG) and the momentum-
based decentralized policy gradient (MDPG, described in Algorithm 2 in the
supplementary materials), which is a new variant created in this work for the
purpose of empirical comparison. We show that our theoretical claims are valid
based on the experiments.

Table 1: Comparisons between existing and proposed approaches.

Method | Complexity | Decentralized | Variance Reduction | Linear Speed Up | I.S.
---|---|---|---|---|---
MBPG (Huang et al. 2020) | $\mathcal{O}(\epsilon^{-3})$ | ✗ | ✓ | ✗ | ✓
Value Prop (Qu et al. 2019) | $\mathcal{O}(\epsilon^{-4})$ | ✓ | ✗ | ✗ | ✗
DCPG (Zeng et al. 2020) | $\mathcal{O}(\epsilon^{-4})$ | ✓ | ✗ | ✗ | ✗
Safe-Dec-PG (Lu et al. 2021) | $\mathcal{O}(\epsilon^{-4})$ | ✓ | ✓ | ✗ | ✗
DFQI (Zhang et al. 2021c) | $\mathcal{O}(\epsilon^{-4})$ | ✓ | ✓ | ✗ | ✗
Dec-TD(0)+GT (Lin and Ling 2021) | N/A | ✓ | ✓ | ✗ | ✗
MDPGT (ours) | $\mathcal{O}(\epsilon^{-4})$ | ✓ | ✓ | ✓ | ✓
MDPGT-MI (ours) | $\mathcal{O}(\epsilon^{-3})$ | ✓ | ✓ | ✓ | ✓

  * 1

Complexity: Sampling complexity for achieving $\mathbb{E}[\|\nabla
J(\mathbf{x})\|]\leq\epsilon$.

  * 2

Linear Speed Up: If an algorithm has $\mathcal{O}(1/\sqrt{K})$ convergence,
then its sampling complexity of attaining an $\mathcal{O}(\epsilon)$ accurate
solution is $\epsilon^{-2}$. Similarly, $\mathcal{O}(1/\sqrt{NK})$ corresponds
to $N^{-1}\epsilon^{-2}$, which is $N$ times faster than the former.
Typically, $K$ has to satisfy a certain condition.

  * 3

MDPGT-MI is MDPGT with mini-batch initialization. We use this notation for
conveniently classifying two different initialization approaches. In the rest
of paper, we still adopt MDPGT to unify these two approaches.

  * 4

I.S. denotes the utilization of importance sampling.

Related Works. Most previous decentralized MARL papers (Zhang et al. 2018;
Suttle et al. 2020; Li et al. 2020; Chen et al. 2020; Bono et al. 2018) tend
to focus on convergence to the optimal return. Exceptions include (Qu et al.
2019), where they proved non-asymptotic convergence rates with nonlinear
function approximation using value propagation. This enables us to
approximately derive the number of stochastic gradient evaluations. However,
the algorithm involves the complex inner-outer structure and requires the size
of the mini-batch to be $\sqrt{K}$, with $K$ being the number of iterations,
which may not be practically implementable. Zhang et al. (2021c) obtain
$\mathcal{O}(\epsilon^{-4})$ for the cooperative setting by using gradient
tracking (GT), which is a bias correction technique dedicated to decentralized
optimization, but with several specifically imposed assumptions, such as
stationary sample paths, which may not be realistic. Lu et al. (2021) also
utilize GT but require dual parameter updates to achieve
$\mathcal{O}(\epsilon^{-4})$; our approach is different and simpler. In this
context, we mention that centralized counterparts of MARL (Huang et al. 2020;
Xu, Gao, and Gu 2019; Papini et al. 2018) have achieved sample complexity of
$\mathcal{O}(\epsilon^{-3})$. However, in both Xu, Gao, and Gu (2019) and
Papini et al. (2018), the size of mini-batch is $\epsilon$-related, which is
more computationally sophisticated than those in both (Huang et al. 2020) and
our proposed method. We provide additional discussion of related work in the
supplementary materials.

## Preliminaries

We first formulate MARL, followed by an overview of variance reduction
techniques and decentralized policy gradients.

### MARL Formulation

In this context, we consider a networked system involving multiple agents (say
$N$) that collaboratively solve dynamic optimization problems. Specifically,
the system can be quantified as a graph, i.e.,
$\mathcal{G}=(\mathcal{V},\mathcal{E})$, where $\mathcal{V}=\\{1,2,...,N\\}$
is the vertex set and $\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$ is
the edge set. Throughout the paper, we assume that $\mathcal{G}$ is static and
undirected, though in a few previous works (Lin et al. 2019; Suttle et al.
2020; Zhang et al. 2018) $\mathcal{G}$ could be directed. The goal of this
work is to provide the rigorous theoretical analysis for our decentralized
MARL algorithm, with the property of the graph not being the main focus. When
a pair of agents $i$ and $j$ can communicate with each other, we have
$(i,j)\in\mathcal{E}$. We also define the neighborhood of a specific agent
$i$, $Nb(i)$, such that
$Nb(i)\triangleq\\{j|j\in\mathcal{V},(i,j)\in\mathcal{E}\textnormal{or}\;j=i\\}$.
Only agents in $Nb(i)$ are able to communicate with the agent $i$. We next
present the definition of networked MARL on top of $\mathcal{G}$.

With multiple agents, the networked Markov decision process is thus
characterized by a tuple
$(\mathcal{S},\\{\mathcal{A}^{i}\\}_{i\in\mathcal{V}},\mathcal{P},\\{r^{i}\\}_{i\in\mathcal{V}},\mathcal{G},\gamma)$,
where $\mathcal{S}$ indicates a global state space shared by all agents in
$\mathcal{V}$ with $|\mathcal{S}|<\infty$, $\mathcal{A}^{i}$ signifies the
action space specified for agent $i$, and $\gamma\in(0,1]$ is the discount
factor. Moreover, in the cooperative MARL setting, the environment is driven
by the joint action space instead of individual action spaces. Thus,
$\mathcal{A}\triangleq\prod_{i\in\mathcal{V}}\mathcal{A}^{i}$ is defined as
the joint action space over all agents in $\mathcal{V}$.
$\mathcal{P}:\mathcal{S}\times\mathcal{A}\times\mathcal{S}\to[0,1]$ represents
the probability function to transition the current state to the next state.
$\\{r^{i}\\}_{i\in\mathcal{V}}:\mathcal{S}\times\mathcal{A}\to\mathbb{R}$ is
the local reward function of agent $i$ and $r^{i}\in[-R,R](R>0)$.
Additionally, states and actions are assumed to be globally observable, while
the rewards are only locally observable. Such an assumption corresponds to the
definition of the cooperative MARL setting and has been generic in previous
works (Zhang et al. 2018; Zhang, Yang, and Basar 2018; Zhang et al. 2021c).

We next describe how agents behave in such an environment. Suppose that the
current state of the environment is $s_{k}\in\mathcal{S}$, where $k$ is the
time step. Each agent $i$ chooses its own action
$a^{i}_{k}\in\mathcal{A}^{i}$, based on the local policy,
$\pi^{i}:\mathcal{S}\times\mathcal{A}^{i}\to[0,1]$. For a parameterized
policy, we denote by $\pi^{i}_{\mathbf{x}^{i}}(s,a^{i})$, which indicates the
probability of agent $i$ choosing action $a^{i}$ given the current state $s$
and $\mathbf{x}^{i}\in\mathbb{R}^{d_{i}}$ here is the policy parameter.
Stacking all local policy parameter together yields
$\mathbf{x}=[(\mathbf{x}^{1})^{\top},(\mathbf{x}^{2})^{\top},...,(\mathbf{x}^{N})^{\top}]^{\top}\in\mathbb{R}^{\sum_{i\in\mathcal{V}}d_{i}}$.
Hence, the joint policy can be denoted as
$\pi_{\mathbf{x}}:\mathcal{S}\times\mathcal{A}\to[0,1]$, where
$\pi_{\mathbf{x}}(s,a)\triangleq\prod_{i\in\mathcal{V}}\pi^{i}_{\mathbf{x}^{i}}(s,a^{i})$
and $a\in\mathcal{A}$. In this context, the decisions are decentralized due to
locally observable rewards, locally evaluated policies and locally executed
actions. To simplify the notations, we drop the $\mathbf{x}^{i}$ for
$\pi^{i}_{\mathbf{x}^{i}}$ and $\mathbf{x}$ for $\pi_{\mathbf{x}}$
respectively for local and joint policies throughout the rest of the paper.
With the joint policy $\pi$ and the state transition function $\mathcal{P}$,
the environment evolves from $s$ to $s^{\prime}$ with the probability
$\mathcal{P}(s^{\prime}|s,a)$. Another assumption imposed in this paper for
the policy function is that for all
$i\in\mathcal{V},s\in\mathcal{S},a^{i}\in\mathcal{A}^{i}$, $\pi^{i}(s,a^{i})$
is continuously differentiable w.r.t. all
$\mathbf{x}^{i}\in\mathbb{R}^{d_{i}}$. Such an assumption will assist in
characterizing the smoothness of the objective function.

The goal for each agent is to learn a local policy $\pi^{i}_{*}$ such that the
joint policy $\pi_{*}$ is able to maximize the global average of expected
cumulative discounted rewards, i.e.,

$\pi_{*}=\textnormal{argmax}_{\mathbf{x}\in\mathbb{R}^{d}}J(\mathbf{x})\triangleq\frac{1}{N}\sum_{i\in\mathcal{V}}\mathbb{E}\bigg{[}\sum^{H}_{h=0}\gamma^{h}r^{i}_{h}\bigg{]},$
(1)

where $H$ is the horizon and $d=\sum_{i\in\mathcal{V}}d_{i}$. Several works
(Zhang et al. 2018; Qu et al. 2019; Zhang et al. 2021c; Lin et al. 2019;
Suttle et al. 2020) have made their attempts to resolve this optimization
problem, leading to different algorithms. Since each agent only has access to
local information, a communication protocol needs to be introduced in the
system, as done in decentralized optimization. With that, well-known
centralized policy-based algorithms can be extended as MARL algorithms.
Nevertheless, one issue that has not been sufficiently explored is the
inherent policy gradient variance, which could even be more significant in the
MARL algorithms. Consequently, this work propose novel MARL algorithms to
investigate how to reduce the policy gradient variance during the
optimization.

### Variance Reduction and Bias Correction

In stochastic optimization, variance reduction techniques have been well
studied and applied to either centralized or decentralized gradient descent
type of algorithms, such as SVRG (Johnson and Zhang 2013), SARAH (Nguyen et
al. 2017), SPIDER (Fang et al. 2018), Hybrid-SARAH (Tran-Dinh et al. 2019) and
STORM (Cutkosky and Orabona 2019). In another line of work, the GT technique
(Pu and Nedić 2020; Sun, Daneshmand, and Scutari 2019) was proposed
specifically for consensus-based decentralized optimization techniques to
improve the convergence rate by tracking and correcting each agent’s locally
aggregated gradients. In our work, we leverage both Hybrid-SARAH and GT to
reduce the policy gradient variance and correct the policy gradient bias
respectively in the MARL and achieve the best convergence rate. Hybrid-SARAH
performs with a trade-off parameter to balance the effect between vanilla
stochastic gradient and SARAH. More detail on these techniques are elaborated
in the supplementary materials.

So far, we are not aware of existing results that have successfully shown
SARAH or Hybrid-SARAH type of variance reduction techniques well suited for
decentralized non-oblivious learning problems, e.g., MARL. Consequently, the
regular Hybrid-SARAH technique cannot be directly applied to MARL; we address
this challenge in sections below.

### Decentralized Policy Gradient

Given a time horizon $H$, we define a trajectory specifically for agent $i$,
as $\tau^{i}\triangleq\\{s_{0},a^{i}_{0},...,s_{H-1},a^{i}_{H-1}\\}$ under any
stationary policy. By following the trajectory $\tau^{i}$, a cumulative
discounted reward is given as
$\mathcal{R}_{i}(\tau^{i})\triangleq\sum_{h=0}^{H}\gamma^{h}r^{i}_{h}$ such
that an individual return can be obtained as:

$J_{i}(\mathbf{x}^{i})\triangleq\mathbb{E}_{\tau^{i}\sim
p_{i}(\tau^{i}|\mathbf{x}^{i})}[\mathcal{R}_{i}(\tau^{i})]=\int\mathcal{R}_{i}(\tau^{i})p_{i}(\tau^{i}|\mathbf{x}^{i})d\tau^{i},$
(2)

where $p_{i}(\tau^{i}|\mathbf{x}^{i})$ is the probability distribution over
$\tau^{i}$ that is equivalent to the following expression given the initial
distribution $\rho_{0}^{i}=\rho^{i}(s_{0})$. Without loss of generalization,
we can assume that the initial distribution is identical for all agents,
namely $\rho(s_{0})$. Then, we have,

$p_{i}(\tau^{i}|\mathbf{x}^{i})=\rho_{0}(s_{0})\prod_{h=0}^{H-1}\mathcal{P}(s_{h+1}|s_{h},a^{i}_{h})\pi^{i}(a^{i}_{h}|s_{h}).$
(3)

For each agent $i$, the goal is to find an optimal policy $\pi^{i}_{*}$ to
maximize the return $J_{i}(\mathbf{x}^{i})$. As discussed above, the
underlying dynamic distribution results in a non-oblivious learning problem,
which is more significant in MARL. To resolve this issue, the decentralized
policy gradient is a decent choice. As background knowledge of MARL, we next
present how to arrive at the local stochastic policy gradient, which will help
characterize the analysis for the proposed algorithms.

Computing the gradient of $J_{i}(\mathbf{x}^{i})$ w.r.t $\mathbf{x}^{i}$
yields the following formula:

$\begin{split}\nabla
J_{i}(\mathbf{x}^{i})&=\int\mathcal{R}_{i}(\tau^{i})\frac{\nabla
p_{i}(\tau^{i}|\mathbf{x}^{i})}{p_{i}(\tau^{i}|\mathbf{x}^{i})}p_{i}(\tau^{i}|\mathbf{x}^{i})d\tau^{i}\\\
&=\mathbb{E}_{\tau^{i}\sim
p_{i}(\tau^{i}|\mathbf{x}^{i})}[\nabla\textnormal{log}p_{i}(\tau^{i}|\mathbf{x}^{i})\mathcal{R}_{i}(\tau^{i})]\end{split}$
(4)

In practice, $p_{i}(\tau^{i}|\mathbf{x}^{i})$ is typically unknown such that
the accurate full policy gradient for agent $i$ is difficult to obtain. Thus,
similar to decentralized stochastic gradient descent (Jiang et al. 2017), we
calculate the policy gradient by sampling a mini-batch of trajectories
$\mathcal{B}=\\{\tau^{i}_{m}\\}^{|\mathcal{B}|}_{m=1}$ from the distribution
$p_{i}(\tau^{i}|\mathbf{x}^{i})$ such that

$\hat{\nabla}J_{i}(\mathbf{x}^{i})=\frac{1}{|\mathcal{B}|}\sum_{m\in\mathcal{B}}\nabla\textnormal{log}p_{i}(\tau^{i}_{m}|\mathbf{x}^{i})\mathcal{R}_{i}(\tau^{i}_{m}).$
(5)

In addition, combining Eq. 3, we can observe that
$\nabla\textnormal{log}p_{i}(\tau^{i}_{m}|\mathbf{x}^{i})$ is independent of
the probability transition $\mathcal{P}$. Hence, Eq. 5 is written as

$\begin{split}\hat{\nabla}J_{i}(\mathbf{x}^{i})&=\frac{1}{|\mathcal{B}|}\sum_{m\in\mathcal{B}}\mathbf{g}_{i}(\tau^{i}_{m}|\mathbf{x}^{i})\\\
&=\frac{1}{|\mathcal{B}|}\sum_{m\in\mathcal{B}}\bigg{(}\sum_{h=0}^{H-1}\nabla_{\mathbf{x}^{i}}\textnormal{log}\pi^{i}(a^{i,m}_{h},s^{m}_{h})\bigg{)}\cdot\\\
&\bigg{(}\sum_{h=0}^{H-1}\gamma^{h}r^{i}_{h}(a^{i,m}_{h},s^{m}_{h})\bigg{)}\end{split}$
(6)

In the above equation, $\mathbf{g}_{i}(\tau^{i}|\mathbf{x}^{i})$ is the
unbiased estimate of $\nabla J_{i}(\mathbf{x}^{i})$, i.e.,
$\mathbb{E}[\mathbf{g}^{i}(\tau^{i}|\mathbf{x}^{i})]=\nabla
J_{i}(\mathbf{x}^{i})$. Some well-known policy gradient estimators can be
obtained through Eq. 6, such as decentralized REINFORCE, which is the direct
extension of its centralized version. We refer interested readers to (Huang et
al. 2020) for more details.

## Our Proposed Approach: MDPGT

### Hybrid Importance Sampling SARAH

In this subsection, we propose a hybrid importance sampling version of SARAH,
termed HIS-SARAH, for decentralized policy gradient updates. First, we define
the importance sampling weight (Metelli et al. 2020) as follows:

$\upsilon(\tau|\mathbf{x}^{\prime},\mathbf{x})=\frac{p(\tau|\mathbf{x}^{\prime})}{p(\tau|\mathbf{x})}=\prod_{h=0}^{H-1}\frac{\pi_{\mathbf{x}^{\prime}}(a_{h}|s_{h})}{\pi_{\mathbf{x}}(a_{h}|s_{h})}.$
(7)

As mentioned in the last section, due to the non-oblivious learning problem,
$\mathbb{E}_{\tau\sim
p(\tau|\mathbf{x})}[\mathbf{g}(\tau|\mathbf{x})-\mathbf{g}(\tau|\mathbf{x}^{\prime})]\neq\nabla
J(\mathbf{x})-\nabla J(\mathbf{x}^{\prime})$. With Eq. 7 we have
$\mathbb{E}_{\tau\sim
p(\tau|\mathbf{x})}[\mathbf{g}(\tau|\mathbf{x})-\upsilon(\tau|\mathbf{x}^{\prime},\mathbf{x})\mathbf{g}(\tau|\mathbf{x}^{\prime})]=\nabla
J(\mathbf{x})-\nabla J(\mathbf{x}^{\prime})$, which has been analyzed in
(Huang et al. 2020) for centralized policy optimization methods and will be a
key relationship in our proof. We denote by $\mathbf{u}^{i}$ the stochastic
policy gradient surrogate for agent $i$. Thus, applying Eq. 7 in a
decentralized manner for Hybrid-SARAH (See Supplementary materials for
definition) gives the following update law at a time step $k$:

$\begin{split}\mathbf{u}^{i}_{k}&=\beta\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})+(1-\beta)[\mathbf{u}_{k-1}^{i}+\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})\\\
&-\upsilon_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{x}^{i}_{k})\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1})].\end{split}$
(8)

The second term on the right hand side of Eq. 8 differ in the extra importance
sampling weight compared to Eq. 13 in the supplementary materials.
Intuitively,
$\upsilon_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{x}^{i}_{k})$ resolves
the non-stationarity in the MARL and retains the regular variance reduction
property of HIS-SARAH as applied in supervised learning problems. Clearly,
each $\mathbf{u}_{k}^{i}$ is a conditionally biased estimator $\nabla
J_{i}(\mathbf{x}^{i}_{k})$, i.e., $\mathbb{E}[\mathbf{u}_{k}^{i}]\neq\nabla
J_{i}(\mathbf{x}^{i}_{k})$ typically. Nevertheless, it can be shown that
$\mathbb{E}[\mathbf{u}_{k}^{i}]=\mathbb{E}[\nabla J_{i}(\mathbf{x}^{i}_{k})]$,
which implies that $\mathbf{u}_{k}^{i}$ acts as a surrogate for the underlying
exact full policy gradient. Therefore, $\mathbf{u}_{k}^{i}$ will be called
directly the stochastic policy gradient surrogate for the rest of analysis.
With Eq. 8 in hand, we now are ready to present the algorithmic framework in
the following.

### Algorithmic Framework

We first present MDPGT (in Algorithm 1), which only takes a trajectory to
initialize the policy gradient surrogate, leading to significant randomness
due to the conditionally biased estimator property at the starting point of
optimization, but still retaining the same sampling complexity as compared to
the SOTA of MARL. To have a better initialization, we also present another way
of initialization by sampling a mini-batch of trajectories from the
distribution (in blue in Algorithm 1). Surprisingly, we will see that with a
proper size of mini-batch initialization, the sampling complexity of our
proposed approach complies with the best result in centralized RL, which
improves the SOTA of MARL.

Result: $\tilde{\mathbf{x}}_{K}$ chosen uniformly random from
$\\{\mathbf{x}^{i}_{k},i\in\mathcal{V}\\}^{K}_{k=1}$

Input:
$\mathbf{x}^{i}_{0}=\bar{\mathbf{x}}_{0}\in\mathbb{R}^{d},\eta\in\mathbb{R}^{+},\beta\in(0,1),\mathbf{W}\in\mathbb{R}^{N\times
N},\mathbf{v}^{i}_{0}=\mathbf{0}_{d},\mathbf{u}_{-1}^{i}=\mathbf{0}_{d},K,\mathcal{B}\in\mathbb{Z}^{+},k=1$;

Initialize the local policy gradient surrogate by sampling a trajectory
$\tau^{i}_{0}$ from
$p_{i}(\tau^{i}|\mathbf{x}^{i}_{0}):\mathbf{u}^{i}_{0}=\mathbf{g}_{i}(\tau^{i}_{0}|\mathbf{x}^{i}_{0})$,
or by sampling a mini-batch of trajectories
$\\{\tau_{0}^{i,m}\\}_{m=1}^{|\mathcal{B}|}$ from
$p_{i}(\tau^{i}|\mathbf{x}^{i}_{0}):\mathbf{u}^{i}_{0}=\frac{1}{|\mathcal{B}|}\sum_{m=1}^{|\mathcal{B}|}\mathbf{g}_{i}(\tau^{i,m}_{0}|\mathbf{x}^{i}_{0})$;

Initialize the local policy gradient tracker: $\mathbf{v}^{i}_{1}=\sum_{j\in
Nb(i)}\omega_{ij}\mathbf{v}^{j}_{0}+\mathbf{u}^{i}_{0}-\mathbf{u}^{i}_{-1}$;

Initialize the local estimate of the policy network parameter:
$\mathbf{x}^{i}_{1}=\sum_{j\in
Nb(i)}\omega^{ij}(\mathbf{x}^{j}_{0}+\eta\mathbf{v}^{j}_{1})$;

while  _$k <K$_ do

for  _each agent_ do

Sample a trajectory $\tau^{i}_{k}$ from $p_{i}(\tau^{i}|\mathbf{x}^{i}_{k})$
and compute the local policy gradient surrogate using Eq. 8;

Update the local policy gradient tracker $\mathbf{v}^{i}_{k+1}=\sum_{j\in
Nb(i)}\omega_{ij}\mathbf{v}^{j}_{k}+\mathbf{u}^{i}_{k}-\mathbf{u}^{i}_{k-1}$;

Update the local estimate of the policy network parameters
$\mathbf{x}^{i}_{k+1}=\sum_{j\in
Nb(i)}\omega_{ij}(\mathbf{x}^{j}_{k}+\eta\mathbf{v}^{j}_{k+1})$;

end for

$k=k+1$;

end while

Algorithm 1 MDPGT

A brief outline of Algorithm 1 is as follows. The initialization of the policy
gradient surrogate $\mathbf{u}^{i}_{0}$ can either be based on only a
trajectory sampled from $p_{i}(\tau^{i}|\mathbf{x}^{i}_{0})$ or a mini-batch.
Subsequently, the policy gradient tracker and network parameters are
initialized based on $\mathbf{u}^{i}_{0}$. The core part of the algorithm
consists of each individual update for
$\mathbf{u}^{i}_{k},\mathbf{v}^{i}_{k}$, and $\mathbf{x}^{i}_{k}$. By
controlling the value of $\beta$ in Eq. 8, MDPGT can degenerate to either
vanilla decentralized policy gradient (with $\beta=1$) or decentralized
version of SRVR-PG (Xu, Gao, and Gu 2019) (with $\beta=0$), both with the
gradient tracking step. In our work, to emphasize the impact of the trade-off
on the policy gradient surrogate, we keep $\beta\in(0,1)$, which makes $\beta$
act more closely as the momentum coefficient in accelerated SGD algorithms
(Singh et al. 2020).

We emphasize that we are unaware of theoretical results for decentralized
SRVR-PG. Hence, the proof techniques presented in this paper can also apply to
this case. Another implication from Algorithm 1 is that at the beginning of
each time step $k$, only one trajectory is required for computing the policy
gradient, allowing for the batch size to be independent of $\epsilon$, i.e.,
$\mathcal{O}(1)$, where we omit the number of agents $N$ when considering the
whole networked system.

## Theoretical Convergence

In this section, we present an analysis of MDPGT. Most of the assumptions
below are mild, and standard in the decentralized optimization and RL
literature. Due to space limitations, we defer auxiliary lemmas and proofs to
the supplementary materials.

######  Assumption 1.

Gradient and Hessian matrix of function $\textnormal{log}\pi^{i}(a^{i}|s)$ are
bounded, i.e., there exist constants $C_{g},C_{h}>0$ such that
$\|\nabla\textnormal{log}\pi^{i}(a^{i}|s)\|\leq C_{g}$ and
$\|\nabla^{2}\textnormal{log}\pi^{i}(a^{i}|s)\|\leq C_{h}$, for all
$i\in\mathcal{V}$.

Note that we skip the subscript $\mathbf{x}^{i}$ at $\pi^{i}$ for the notation
simplicity. In this context, we did not impose the bounded policy gradient
assumption, though it can be derived based on the above assumption, which has
been adopted in centralized RL algorithms (Zhang et al. 2021b; Huang et al.
2020; Xu, Gao, and Gu 2019). Additionally, it also helps derive the smoothness
of $J_{i}(\mathbf{x}^{i})$ that has typically been exposed as an assumption in
decentralized learning/optimization literature.

######  Assumption 2.

The mixing matrix $\mathbf{W}\in\mathbb{R}^{N\times N}$ is doubly stochastic
such that $\lambda\triangleq\|\mathbf{W}-\mathbf{P}\|\in[0,1)$, where
$\lambda$ signifies the second largest eigenvalue to measure the algebraic
connectivity of the graph, and
$\mathbf{P}=\frac{1}{N}\mathbf{1}^{\top}\mathbf{1}$ and $\mathbf{1}$ is a
column vector with each entry being 1.

######  Assumption 3.

Variance of importance sampling weight
$\upsilon_{i}(\tau^{i}|\mathbf{x}_{1},\mathbf{x}_{2})$ is bounded, i.e., there
exists a constant $\mathcal{M}>0$ such that
$\mathbb{V}(\upsilon_{i}(\tau^{i}|\mathbf{x}_{1},\mathbf{x}_{2}))\leq\mathcal{M}$,
for any $\mathbf{x}_{1},\mathbf{x}_{2}\in\mathbb{R}^{d_{i}}$ and $\tau^{i}\sim
p_{i}(\tau^{i}|\mathbf{x}_{1})$, for all $i\in\mathcal{V}$.

Assumption 2 is generic in most previous works on decentralized optimization,
though such a property has been relaxed in some works (Nedić and Olshevsky
2014). However, we have not been aware of any existing works in MARL doing
such a relaxation and its investigation can be of independent interest.
Assumption 3 is specifically introduced for importance sampling-based methods.
Such an assumption is critical to construct the relationship between
$\mathbb{V}(\upsilon_{i}(\tau^{i}|\mathbf{x}_{1},\mathbf{x}_{2}))$ and
$\|\mathbf{x}_{1}-\mathbf{x}_{2}\|^{2}$, through which the impact of the
variance of importance sampling on the convergence can be explicitly
quantified. Another typical assumption is for the bounded variance of
stochastic gradient such that
$\mathbb{E}[\|\mathbf{g}_{i}(\tau^{i}|\mathbf{x}^{i})-\nabla
J_{i}(\mathbf{x}^{i})\|^{2}]\leq\sigma^{2}_{i}$. However, under MARL setting,
such a result can be derived from Assumption 1 and we present the formal
result in Lemma 1. In this context, we also have
$\bar{\sigma}^{2}=\frac{1}{N}\sum^{N}_{i=1}\sigma^{2}_{i}$, for all
$i\in\mathcal{V}$. The explicit expression of $\bar{\sigma}^{2}$ is given in
the supplementary materials.

### Main Results

We present the main results to show specifically the convergence rates for
MDPGT when it is initialized by a mini-batch of trajectories. We denote by
$L>0$ the smoothness constant and $G>0$ the upper bound of
$\|\mathbf{g}_{i}(\tau^{i}|\mathbf{x}^{i})\|$ for all $i\in\mathcal{V}$. We
further define a constant $C_{\upsilon}>0$ such that
$\mathbb{V}(\upsilon_{i}(\tau^{i}|\mathbf{x}_{1},\mathbf{x}_{2}))\leq
C^{2}_{\upsilon}\|\mathbf{x}_{1}-\mathbf{x}_{2}\|^{2}$. The explicit
expressions of these constants are derived in lemmas in the supplementary
materials. Note that in our work, $G$ is not directly assumed, but instead
derived based on Assumption 1.

######  Theorem 1.

Let Assumptions 1,2 and 3 hold. Let the momentum coefficient
$\beta=\frac{96L^{2}+96G^{2}C_{\upsilon}}{N}\eta^{2}$. If MDPGT is initialized
by a mini-batch of trajectories with the size being $\mathcal{B}$ and the step
size satisfies the following condition

$\begin{split}0<\eta\leq\textnormal{min}&\bigg{\\{}\frac{(1-\lambda^{2})^{2}}{\lambda\sqrt{12844L^{2}+9792G^{2}C_{\upsilon}}},\frac{\sqrt{N(1-\lambda^{2})}\lambda}{31\sqrt{L^{2}+G^{2}C_{\upsilon}}},\\\
&\frac{1}{6\sqrt{6(L^{2}+G^{2}C_{\upsilon})}}\bigg{\\}},\end{split}$

then the output $\tilde{\mathbf{x}}_{K}$ satisfies: for all $K\geq 2$:

$\begin{split}&\mathbb{E}[\|\nabla
J(\tilde{\mathbf{x}}_{K})\|^{2}]\leq\frac{4(J^{*}-J(\bar{\mathbf{x}}_{0}))}{\eta
K}+\frac{4\bar{\sigma}^{2}}{N|\mathcal{B}|\beta
K}+\frac{8\beta\bar{\sigma}^{2}}{N}\\\
&+\frac{34\lambda^{2}}{KN(1-\lambda^{2})^{3}}\|\nabla\mathbf{J}(\bar{\mathbf{x}}_{0})\|^{2}+\frac{68\lambda^{2}\bar{\sigma}^{2}}{(1-\lambda^{2})^{3}|\mathcal{B}|K}+\frac{204\lambda^{2}\beta^{2}\bar{\sigma}^{2}}{(1-\lambda^{2})^{3}},\end{split}$
(9)

where $J^{*}$ is the upper bound of $J(\mathbf{x})$ and
$\|\nabla\mathbf{J}(\bar{\mathbf{x}}_{0})\|^{2}\triangleq\sum^{N}_{i=1}\|\nabla
J_{i}(\bar{\mathbf{x}}_{0})\|^{2}$.

Theorem 1 depicts that when $K\to\infty$, MDPGT enables convergence to a
steady-state error in a sublinear rate $\mathcal{O}(1/K)$ if $\eta$ and
$\beta$ are selected properly, i.e.,

$\mathbb{E}[\|\nabla
J(\tilde{\mathbf{x}}_{K})\|^{2}]\leq\frac{8\beta\bar{\sigma}^{2}}{N}+\frac{204\lambda^{2}\beta^{2}\bar{\sigma}^{2}}{(1-\lambda^{2})^{3}}.$
(10)

By observing Eq. 10, the steady-state error is determined by the number of
agents, the variance of stochastic policy gradient, and the spectral gap of
the graph $1-\lambda$. Increasing the number of agents leads to a small error
bound. Though different network topologies imply different error bounds, the
higher order term of $\beta$ can reduce the impact of the spectral gap on the
error bound. Another suggestion from Eq. 10 is that $\eta$ and $\beta$ can be
reduced to make the steady-state error arbitrarily small, while in return this
can affect the speed of convergence. Surprisingly, even though we have to
adopt the bounded stochastic policy gradient derived from Assumption 1 for
analysis, the error bound in Eq. 10 only depends heavily on the variance,
which is inherently consistent with most conclusions from decentralized
optimization in literature without the bounded stochastic gradient assumption.
While $J^{*}$ is essentially correlated to the upper bound of reward $R$, in
this context, we still adopt the implicit $J^{*}$ for convenience. We next
provide the analysis for the non-asymptotic behavior, defining appropriately
$\eta,\beta$, and $|\mathcal{B}|$.

######  Corollary 1.

Let
$\eta=\frac{N^{2/3}}{8LK^{1/3}},\beta=\frac{DN^{1/3}}{64L^{2}K^{2/3}},|\mathcal{B}|=bigg\lceil\frac{K^{1/3}}{N^{2/3}}bigg\rceil$
in Theorem 1. We have,

$\begin{split}&\mathbb{E}[\|\nabla
J(\tilde{\mathbf{x}}_{K})\|^{2}]\leq\underbrace{\frac{256L^{3}D(J^{*}-J(\bar{\mathbf{x}}_{0}))+2048L^{4}\bar{\sigma}^{2}+D^{2}\bar{\sigma}^{2}}{8L^{2}D(NK)^{2/3}}}_{T_{1}}\\\
&+\underbrace{\frac{34\lambda^{2}}{KN(1-\lambda^{2})^{3}}\|\nabla\mathbf{J}(\bar{\mathbf{x}}_{0})\|^{2}+\frac{\lambda^{2}\bar{\sigma}^{2}(51D^{2}+69632N^{2/3}L^{4})}{1024(1-\lambda^{2})^{3}K^{4/3}L^{4}}}_{T_{2}},\end{split}$
(11)

for all

$\begin{split}&K\geq\textnormal{max}\bigg{\\{}\frac{N^{2}D^{1.5}}{512L^{3}},\frac{29791\sqrt{N}(L^{2}+G^{2}C_{\upsilon})^{1.5}}{512L^{3}\lambda^{3}(1-\lambda^{2})^{1.5}},\\\
&\frac{(12844L^{2}+9792G^{2}C_{\upsilon})^{1.5}N^{2}\lambda^{3}}{512L^{3}(1-\lambda^{2})^{6}}\bigg{\\}},\end{split}$

where $D=96L^{2}+96G^{2}C_{\upsilon}$.

######  Remark 1.

An implication from Corollary 1 is that at the early stage of optimization,
before $T_{1}$ in Eq. 11 dominates, the complexity is tightly related to the
algebraic connectivity of the network topology in $T_{2}$, which is measured
by the spectral gap $1-\lambda$. However, $T_{2}$ is in a large order of
$1/K$. As the optimization moves towards the latter stage where $T_{1}$
dominates, the overall complexity is independent of the network topology.

For the ease of exposition, with Corollary 1, when $K$ is sufficiently large,
it is an immediate consequence as $\mathbb{E}[\|\nabla
J(\tilde{\mathbf{x}}_{K})\|^{2}]\leq\mathcal{O}((NK)^{-2/3})$. Thus, for
achieving $\mathbb{E}[\|\nabla J(\tilde{\mathbf{x}}_{K})\|]\leq\epsilon$, the
following relationship is obtained:

$\displaystyle\mathbb{E}[\|\nabla J(\tilde{\mathbf{x}}_{K})\|]$
$\displaystyle=\sqrt{(\mathbb{E}[\|\nabla
J(\tilde{\mathbf{x}}_{K})\|])^{2}}\leq\sqrt{\mathbb{E}[\|\nabla
J(\tilde{\mathbf{x}}_{K})\|^{2}]}\leq\epsilon.$

Combining the last two inequalities results in the ultimate sampling
complexity, i.e., $\mathcal{O}(N^{-1}\epsilon^{-3})$, which exhibits linear
speed up. More importantly, this is $N$ times smaller than the sampling
complexity of the centralized momentum-based policy gradient methods (Huang et
al. 2020) that performs on a single node. However, we have known from
Corollary 1 that typically $K$ has to be large, which will in the linear speed
up regime reduce $\eta$ and $\beta$.

We also investigate a worse initialization with only a single trajectory
sampled from $p_{i}(\tau^{i}|\mathbf{x}^{i}_{0})$. However, without a mini-
batch initialization, the eventual sampling complexity is
$\mathcal{O}(N^{-1}\epsilon^{-4})$ (see Theorem 2 and Corollary 3). Though
variance reduction techniques have not reduced the order of $\epsilon^{-1}$,
compared to the SOTA approaches, the linear speed up still enables the
complexity to be $N$ times smaller than that in (Xu, Gao, and Gu 2019; Huang
et al. 2020). Additionally, different from traditional decentralized learning
problems, MARL has more significant variances in the optimization procedure
due to the non-oblivious characteristic. Using just a single trajectory for
each agent to initialize is can be a poor scheme, but the adopted variance
reduction techniques can successfully maintain the SOTA sampling complexity in
a decentralized setting. Please refer to the supplementary materials for
formal results and proof.

Implication for Gaussian Policy. We study the sample complexity when the
policy function $\pi^{i}(a^{i}|s)$ of each agent is explicitly a Gaussian
distribution. For a bounded action space $\mathcal{A}^{i}\subset\mathbb{R}$, a
Gaussian policy parameterized by $\mathbf{x}_{i}$ is defined as

$\pi^{i}(a^{i}|s)=\frac{1}{\sqrt{2\pi}}\textnormal{exp}\bigg{(}-\frac{((\mathbf{x}^{i})^{\top}\phi_{i}(s)-a^{i})^{2}}{2\xi^{2}}\bigg{)},$
(12)

where $\xi^{2}$ is a constant standard deviation parameter and
$\phi_{i}(s):\mathcal{S}\to\mathbb{R}^{d_{i}}$ is a mapping from the state
space to the feature space. Thus, the following formal result can be obtained.
A more formal analysis and proof can be seen in the supplementary materials.

######  Corollary 2.

Let $\pi^{i}(a^{i}|s)$ be defined as a Gaussian distribution in Eq. 12 with
$|a^{i}|\leq C_{a}$, where $C_{a},C_{f}>0$, and $\|\phi_{i}(s)\|\leq C_{f}$,
and $\eta,\beta,|\mathcal{B}|$ be defined as in Corollary 1. The sampling
complexity of attaining the accuracy $\mathbb{E}[\|\nabla
J(\tilde{\mathbf{x}}_{K})\|]\leq\epsilon$ is
$\mathcal{O}((1-\gamma)^{-4.5}N^{-1}\epsilon^{-3})$.

## Numerical Experiments and Results

To validate our proposed algorithm, we performed experiments on a multi-agent
benchmark environment with a cooperative navigation task that has been
commonly used as a benchmark in several previous works (Qu et al. 2019; Zhang
et al. 2018; Lu et al. 2021). Our platform for cooperative navigation is
derived off the particle environment introduced by (Lowe et al. 2017). In our
modification, all agents are initialized at random locations with a specific
goal in a 2-dimensional grid world. Each agent observes the combined position
and velocity of itself and all other agents. The agent is capable of moving
up, down, left or right with the objective of navigating to their respective
goals. The reward function of each agent is defined as the negative euclidean
distance of the agent to the goal. Additionally, a penalty of -1 is imposed
whenever the agent collides with any other agents. All agent’s policy is
represented by a 3-layer dense neural network with 64 hidden units with $tanh$
activation functions. The agents were trained for 50,000 episodes with a
horizon of 50 steps and discount factor of 0.99. For the sake of brevity, we
present numerical results in only one environment setting with five agents.
Additional results with different number of agents and a simplified
environment and computing infrastructure details are available in the
supplementary materials111Codes to reproduce these results are available at
the following repository: https://github.com/xylee95/MD-PGT.

### Efficacy of MDPGT

Figure 1: Average rewards obtain by MDPGT, MDPG and DPG in a cooperative
navigation task for five agents. For MDPGT, $\beta$ = 0.5 and results shown
are averaged across five random seeds. The line plots denote the mean value
and shaded regions denote the standard deviation of rewards.

Figure 1 illustrates the average training rewards obtained by the five agents
in the cooperative navigation gridworld environment. As observed, both MDPG
and MDPGT significantly outperforms the baseline, denoted as DPG. Comparing
MDPG with MDPGT, we observe that while both algorithms initially have similar
performance, MDPGT begins to outperform MDPG around the 15,000 iteration.
Additionally, when we compare the standard deviation of the rewards, shown as
shaded regions, we observe that standard deviation of MDPGT is also smaller
than the standard deviation of MDPG. In summary, these results validate our
theoretical findings that utilizing gradient tracking as bias correction
technique does improve the performance of the algorithm. Additionally, the
improvement in terms of sampling complexity over DPG is empirically evident
through the result.

### Effect of Momentum Coefficient

Figure 2: Ablation study illustrating the effect of various momentum
coefficients, $\beta$ on the performance of MDPGT for five agents in the
cooperative navigation environment. Figure 3: Experiment results for five
agents in the cooperative navigation environment to compare the effects of
different network topologies. $\beta$ = 0.5 for all experiments shown.

Next, we perform an additional ablation study to investigate the effect of the
momentum coefficient $\beta$ on the performance of MDPGT. As shown in Figure
2, we see that the choice of momentum coefficient does indeed have an effect
on the performance. A $\beta$ that is low can induce a faster convergence
rate, but at the cost of a higher fluctuations in rewards, as seen by $\beta$
= 0.2 and 0.3. Conversely, a $\beta$ value that is too high will cause the
surrogate to degenerate into vanilla policy gradients and reflects a similar
performance as DPG, which matches the implication by Eq. 10. Ultimately, we
see that for this environment, $\beta$ = 0.4 and 0.5 offers the perfect
balance between convergence rate and stability/variance of the training.
Hence, $\beta$ can be viewed as hyper-parameter which can be tuned to trade
off between optimizing for convergence versus training stability.

### Effect of Different Topologies

Finally, we provide evidence which confirms the fact that our proposed method
holds under various networks topologies. To test our hypothesis, we train five
agents in the same cooperative navigation environment using three different
network topologies: fully-connected, ring and bi-partite topology. As seen in
Figure 3, the five agents achieves similar rewards despite communicating via
different network topologies. This validates our claim in Remark 1.

## Conclusions

This paper proposes a novel MARL algorithm that involves variance reduction
techniques to reduce the sampling complexity of decentralized policy-based
methods. Specifically we have developed the algorithmic framework and analyzed
it in a principled manner. An importance sampling-based stochastic recursive
momentum is presented as the policy gradient surrogate, which is taken as
input to a policy gradient tracker. Through theoretical analysis, we have
found that the proposed method can improve the sampling efficiency in the
decentralized RL settings compared to the SOTA methods. To the best of our
knowledge, this is the first time to achieve the best available rate,
$\mathcal{O}(\epsilon^{-3})$, for generic (possibly non-concave) performance
functions. Empirical results have shown the superiority of the proposed MDPGT
over the baseline decentralized policy gradient methods. Future research
directions include: 1) incorporating more complex decentralized environments
in our experiments to reveal potentially novel and interesting results; 2)
extending the proposed method to model-based decentralized RL settings to
improve further the sampling efficiency; 3) theoretically analyzing the
robustness of the proposed method under adversarial attacks.

## References

  * Ackermann et al. (2019) Ackermann, J.; Gabler, V.; Osa, T.; and Sugiyama, M. 2019.  Reducing overestimation bias in multi-agent domains using double centralized critics.  _arXiv preprint arXiv:1910.01465_.
  * Arslan and Yüksel (2016) Arslan, G.; and Yüksel, S. 2016.  Decentralized Q-learning for stochastic teams and games.  _IEEE Transactions on Automatic Control_ , 62(4): 1545–1558.
  * Bono et al. (2018) Bono, G.; Dibangoye, J. S.; Matignon, L.; Pereyron, F.; and Simonin, O. 2018.  Cooperative multi-agent policy gradient.  In _Joint European Conference on Machine Learning and Knowledge Discovery in Databases_ , 459–476. Springer.
  * Chen et al. (2020) Chen, B.; Xu, M.; Liu, Z.; Li, L.; and Zhao, D. 2020.  Delay-Aware Multi-Agent Reinforcement Learning for Cooperative and Competitive Environments.  _arXiv e-prints_ , arXiv–2005.
  * Chu, Chinchali, and Katti (2020) Chu, T.; Chinchali, S.; and Katti, S. 2020.  Multi-agent reinforcement learning for networked system control.  _arXiv preprint arXiv:2004.01339_.
  * Cortes, Mansour, and Mohri (2010) Cortes, C.; Mansour, Y.; and Mohri, M. 2010.  Learning Bounds for Importance Weighting.  In _Nips_ , volume 10, 442–450. Citeseer.
  * Cutkosky and Orabona (2019) Cutkosky, A.; and Orabona, F. 2019.  Momentum-based variance reduction in non-convex sgd.  _arXiv preprint arXiv:1905.10018_.
  * Das et al. (2020) Das, R.; Acharya, A.; Hashemi, A.; Sanghavi, S.; Dhillon, I. S.; and Topcu, U. 2020\.  Faster Non-Convex Federated Learning via Global and Local Momentum.  _arXiv preprint arXiv:2012.04061_.
  * Fang et al. (2018) Fang, C.; Li, C. J.; Lin, Z.; and Zhang, T. 2018.  Spider: Near-optimal non-convex optimization via stochastic path integrated differential estimator.  _arXiv preprint arXiv:1807.01695_.
  * Helou, Kalathil, and Xie (2020) Helou, R. E.; Kalathil, D.; and Xie, L. 2020.  Fully Decentralized Reinforcement Learning-based Control of Photovoltaics in Distribution Grids for Joint Provision of Real and Reactive Power.  _arXiv preprint arXiv:2008.01231_.
  * Huang et al. (2020) Huang, F.; Gao, S.; Pei, J.; and Huang, H. 2020.  Momentum-based policy gradient methods.  In _International Conference on Machine Learning_ , 4422–4433. PMLR.
  * Jiang et al. (2017) Jiang, Z.; Balu, A.; Hegde, C.; and Sarkar, S. 2017.  Collaborative deep learning in fixed topology networks.  _arXiv preprint arXiv:1706.07880_.
  * Johnson and Zhang (2013) Johnson, R.; and Zhang, T. 2013.  Accelerating stochastic gradient descent using predictive variance reduction.  _Advances in neural information processing systems_ , 26: 315–323.
  * Karimireddy et al. (2020) Karimireddy, S. P.; Jaggi, M.; Kale, S.; Mohri, M.; Reddi, S. J.; Stich, S. U.; and Suresh, A. T. 2020.  Mime: Mimicking centralized stochastic algorithms in federated learning.  _arXiv preprint arXiv:2008.03606_.
  * Li, Chen, and Chen (2020) Li, M.-L.; Chen, S.; and Chen, J. 2020.  Adaptive learning: A new decentralized reinforcement learning approach for cooperative multiagent systems.  _IEEE Access_ , 8: 99404–99421.
  * Li et al. (2020) Li, W.; Jin, B.; Wang, X.; Yan, J.; and Zha, H. 2020.  F2A2: Flexible Fully-decentralized Approximate Actor-critic for Cooperative Multi-agent Reinforcement Learning.  _arXiv preprint arXiv:2004.11145_.
  * Lin and Ling (2021) Lin, Q.; and Ling, Q. 2021.  Decentralized TD (0) with Gradient Tracking.  _IEEE Signal Processing Letters_.
  * Lin et al. (2019) Lin, Y.; Zhang, K.; Yang, Z.; Wang, Z.; Başar, T.; Sandhu, R.; and Liu, J. 2019.  A communication-efficient multi-agent actor-critic algorithm for distributed reinforcement learning.  In _2019 IEEE 58th Conference on Decision and Control (CDC)_ , 5562–5567. IEEE.
  * Lowe et al. (2017) Lowe, R.; Wu, Y.; Tamar, A.; Harb, J.; Abbeel, P.; and Mordatch, I. 2017.  Multi-Agent Actor-Critic for Mixed Cooperative-Competitive Environments.  _Neural Information Processing Systems (NIPS)_.
  * Lu et al. (2021) Lu, S.; Zhang, K.; Chen, T.; Basar, T.; and Horesh, L. 2021.  Decentralized Policy Gradient Descent Ascent for Safe Multi-Agent Reinforcement Learning.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 35, 8767–8775.
  * Metelli et al. (2020) Metelli, A. M.; Papini, M.; Montali, N.; and Restelli, M. 2020.  Importance Sampling Techniques for Policy Optimization.  _Journal of Machine Learning Research_ , 21(141): 1–75.
  * Mukherjee, Bai, and Chakrabortty (2020) Mukherjee, S.; Bai, H.; and Chakrabortty, A. 2020.  Model-free decentralized reinforcement learning control of distributed energy resources.  In _2020 IEEE Power & Energy Society General Meeting (PESGM)_, 1–5. IEEE.
  * Nedić and Olshevsky (2014) Nedić, A.; and Olshevsky, A. 2014.  Distributed optimization over time-varying directed graphs.  _IEEE Transactions on Automatic Control_ , 60(3): 601–615.
  * Nguyen et al. (2014) Nguyen, D. T.; Yeoh, W.; Lau, H. C.; Zilberstein, S.; and Zhang, C. 2014.  Decentralized multi-agent reinforcement learning in average-reward dynamic DCOPs.  In _Proceedings of the AAAI conference on artificial intelligence_ , volume 28.
  * Nguyen et al. (2017) Nguyen, L. M.; Liu, J.; Scheinberg, K.; and Takáč, M. 2017.  SARAH: A novel method for machine learning problems using stochastic recursive gradient.  In _International Conference on Machine Learning_ , 2613–2621. PMLR.
  * Nguyen and Mukhopadhyay (2017) Nguyen, T.; and Mukhopadhyay, S. 2017.  Selectively decentralized Q-learning.  In _2017 IEEE International Conference on Systems, Man, and Cybernetics (SMC)_ , 328–333. IEEE.
  * Nguyen, Nguyen, and Nahavandi (2020) Nguyen, T. T.; Nguyen, N. D.; and Nahavandi, S. 2020.  Deep reinforcement learning for multiagent systems: A review of challenges, solutions, and applications.  _IEEE transactions on cybernetics_ , 50(9): 3826–3839.
  * Papini et al. (2018) Papini, M.; Binaghi, D.; Canonaco, G.; Pirotta, M.; and Restelli, M. 2018.  Stochastic variance-reduced policy gradient.  In _International Conference on Machine Learning_ , 4026–4035. PMLR.
  * Pirotta, Restelli, and Bascetta (2013) Pirotta, M.; Restelli, M.; and Bascetta, L. 2013.  Adaptive step-size for policy gradient methods.  _Advances in Neural Information Processing Systems_ , 26: 1394–1402.
  * Pu and Nedić (2020) Pu, S.; and Nedić, A. 2020.  Distributed stochastic gradient tracking methods.  _Mathematical Programming_ , 1–49.
  * Qu et al. (2019) Qu, C.; Mannor, S.; Xu, H.; Qi, Y.; Song, L.; and Xiong, J. 2019.  Value propagation for decentralized networked deep multi-agent reinforcement learning.  _arXiv preprint arXiv:1901.09326_.
  * Ryu, Shin, and Park (2020) Ryu, H.; Shin, H.; and Park, J. 2020.  Multi-agent actor-critic with hierarchical graph attention network.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 34, 7236–7243.
  * Simões, Lau, and Reis (2020) Simões, D.; Lau, N.; and Reis, L. P. 2020.  Multi-agent actor centralized-critic with communication.  _Neurocomputing_ , 390: 40–56.
  * Singh et al. (2020) Singh, N.; Data, D.; George, J.; and Diggavi, S. 2020.  SQuARM-SGD: Communication-Efficient Momentum SGD for Decentralized Optimization.  _arXiv preprint arXiv:2005.07041_.
  * Sun, Daneshmand, and Scutari (2019) Sun, Y.; Daneshmand, A.; and Scutari, G. 2019.  Convergence rate of distributed optimization algorithms based on gradient tracking.  _arXiv preprint arXiv:1905.02637_.
  * Suttle et al. (2020) Suttle, W.; Yang, Z.; Zhang, K.; Wang, Z.; Başar, T.; and Liu, J. 2020.  A multi-agent off-policy actor-critic algorithm for distributed reinforcement learning.  _IFAC-PapersOnLine_ , 53(2): 1549–1554.
  * Sutton et al. (1999) Sutton, R. S.; McAllester, D. A.; Singh, S. P.; Mansour, Y.; et al. 1999.  Policy gradient methods for reinforcement learning with function approximation.  In _NIPs_ , volume 99, 1057–1063. Citeseer.
  * Tran-Dinh et al. (2019) Tran-Dinh, Q.; Pham, N. H.; Phan, D. T.; and Nguyen, L. M. 2019.  Hybrid stochastic gradient descent algorithms for stochastic nonconvex optimization.  _arXiv preprint arXiv:1905.05920_.
  * Wang et al. (2020) Wang, X.; Wang, C.; Li, X.; Leung, V. C.; and Taleb, T. 2020.  Federated deep reinforcement learning for internet of things with decentralized cooperative edge caching.  _IEEE Internet of Things Journal_ , 7(10): 9441–9455.
  * Wei et al. (2021) Wei, C.-Y.; Lee, C.-W.; Zhang, M.; and Luo, H. 2021.  Last-iterate Convergence of Decentralized Optimistic Gradient Descent/Ascent in Infinite-horizon Competitive Markov Games.  _arXiv preprint arXiv:2102.04540_.
  * Xin, Khan, and Kar (2020) Xin, R.; Khan, U. A.; and Kar, S. 2020.  A near-optimal stochastic gradient method for decentralized non-convex finite-sum optimization.  _arXiv preprint arXiv:2008.07428_.
  * Xin, Khan, and Kar (2021) Xin, R.; Khan, U. A.; and Kar, S. 2021.  A hybrid variance-reduced method for decentralized stochastic non-convex optimization.  _arXiv preprint arXiv:2102.06752_.
  * Xu, Gao, and Gu (2019) Xu, P.; Gao, F.; and Gu, Q. 2019.  Sample efficient policy gradient methods with recursive variance reduction.  _arXiv preprint arXiv:1909.08610_.
  * Zeng et al. (2020) Zeng, S.; Anwar, A.; Doan, T.; Romberg, J.; and Raychowdhury, A. 2020.  A decentralized policy gradient approach to multi-task reinforcement learning.  _arXiv preprint arXiv:2006.04338_.
  * Zhang et al. (2021a) Zhang, J.; Bedi, A. S.; Wang, M.; and Koppel, A. 2021a.  MARL with General Utilities via Decentralized Shadow Reward Actor-Critic.  _arXiv preprint arXiv:2106.00543_.
  * Zhang et al. (2021b) Zhang, J.; Ni, C.; Yu, Z.; Szepesvari, C.; and Wang, M. 2021b.  On the Convergence and Sample Efficiency of Variance-Reduced Policy Gradient Method.  _arXiv preprint arXiv:2102.08607_.
  * Zhang, Yang, and Basar (2018) Zhang, K.; Yang, Z.; and Basar, T. 2018.  Networked multi-agent reinforcement learning in continuous spaces.  In _2018 IEEE Conference on Decision and Control (CDC)_ , 2771–2776. IEEE.
  * Zhang, Yang, and Başar (2019) Zhang, K.; Yang, Z.; and Başar, T. 2019.  Decentralized Multi-Agent Reinforcement Learning with Networked Agents: Recent Advances.  _arXiv preprint arXiv:1912.03821_.
  * Zhang et al. (2018) Zhang, K.; Yang, Z.; Liu, H.; Zhang, T.; and Basar, T. 2018.  Fully decentralized multi-agent reinforcement learning with networked agents.  In _International Conference on Machine Learning_ , 5872–5881. PMLR.
  * Zhang et al. (2021c) Zhang, K.; Yang, Z.; Liu, H.; Zhang, T.; and Basar, T. 2021c.  Finite-sample analysis for decentralized batch multi-agent reinforcement learning with networked agents.  _IEEE Transactions on Automatic Control_.
  * Zhou et al. (2020) Zhou, P.; Chen, X.; Liu, Z.; Braud, T.; Hui, P.; and Kangasharju, J. 2020.  DRLE: Decentralized Reinforcement Learning at the Edge for Traffic Light Control.  _arXiv preprint arXiv:2009.01502_.

## Supplementary Materials

We present the additional analysis and detailed proof for the auxiliary lemmas
and main theorems as well as additional empirical results.

### Additional Related Works

The most recent work (Zhang et al. 2021b) in centralized RL introduced SIVR-
PG, which employed the hidden convex nature of general utility function and
leveraged techniques from composition optimization to attain the sampling
complexity of $\tilde{\mathcal{O}}(\epsilon^{-2})$. Nevertheless, SIVR-PG
requires an additional assumption that the unnormalized state-action occupancy
measure in the utility function is convex and its corresponding inverse
mapping exists, which may not rigorously hold for some RL settings. A
concurrent work (Zhang et al. 2021a) recently that has been published on arXiv
proposes a Decentralized Shadow Reward Actor Critic (DSAC) for general
utilities, which can be treated as the extension to the decentralized setting
of SIVR-PG. In their work, the authors develop the shadow reward that
estimates the derivative of the local utility with respect to their occupancy
measure and show that DSAC converges to $\epsilon$-stationarity in
$\mathcal{O}(\epsilon^{-2.5})$ or even faster $\mathcal{O}(\epsilon^{-2})$
with high probability. They further establish the global optimality by
adopting the diminishing step size. Though DSAC reveals faster convergence
rate, additional occupancy measure needs to be incorporated and updated in the
algorithmic framework. Also, the shadow reward they have defined is
essentially a derivative, which may vary significantly among different agents.
This could be another source of variance during updates. In the analysis, we
have observed that more assumptions are required to arrive at the theoretical
claims, including strongly convex critic objective function. Different from
most previous works, the constructed error bounds in their work are not based
on expectation. The parameters, including step size, batch size, trajectory
length, and even final iteration need to be carefully defined. In the
empirical evaluation, the authors present a cooperative navigation with safety
constraint and do not show any comparison with the decentralized baseline
methods. We therefore are not going to compare our proposed method with theirs
due to different emphasis.

###  Hybrid-SARAH and GT

Hybrid-SARAH. In this context, we denote by $\mathbf{g}(\zeta,\mathbf{x})$ the
stochastic gradient, where $\zeta$ is the random seed. As discussed above,
Hybrid-SARAH is the combination between vanilla stochastic gradient and SARAH
such that the following expression can be obtained for a time step $k$:

$\mathbf{u}_{k}=\beta\mathbf{g}(\zeta_{k},\mathbf{x}_{k})+(1-\beta)(\mathbf{u}_{k-1}+\mathbf{g}(\zeta_{k},\mathbf{x}_{k})-\mathbf{g}(\zeta_{k},\mathbf{x}_{k-1})),$
(13)

where $\beta\in[0,1]$ is the momentum parameter. One can observe from Eq. 13
that in the part of SARAH, the term $\mathbf{g}(\zeta_{k},\mathbf{x}_{k-1})$
requires the computation of $\mathbf{g}(\zeta_{k},\mathbf{x}_{k-1})$ and the
copy of last iteration of parameters $\mathbf{x}_{k-1}$, which could increase
the computational power and memory significantly if the dimension of
$\mathbf{x}$ is quite high, though this may not be an issue for a modern
computer. Additionally, this form of update in SARAH enables to improve the
convergence rates in different problems such as decentralized learning (Xin,
Khan, and Kar 2020) and reinforcement learning (Huang et al. 2020; Xu, Gao,
and Gu 2019). Therefore we will still adapt Hybrid-SARAH in our MARL
algorithms.

GT. Another technique that has specifically been developed for decentralized
optimization is GT, which is close to SARAH, but focusing on tracking and
correcting each agent’s locally aggregated gradients. However, in terms of
formulation, GT typically appears in a form with a consensus or mixing step
that can expressed as follows:

$\mathbf{v}^{i}_{k+1}=\sum_{j\in
Nb(i)}\omega_{ij}\mathbf{v}^{j}_{k}+\mathbf{g}_{i}(\zeta_{k}^{i},\mathbf{x}^{i}_{k})-\mathbf{g}_{i}(\zeta_{k-1}^{i},\mathbf{x}^{i}_{k-1}),$
(14)

where $\omega_{ij}\in\mathbf{W}\in\mathbb{R}^{N\times N}$ is the probability
of an edge in the communication network $\mathcal{G}$. $\mathbf{W}$ is the
mixing matrix to define the topology of the communication network and one
assumption will be imposed for the matrix. An immediate observation from Eq.
14 is that the update requires the copy of
$\mathbf{g}_{i}(\zeta_{k-1}^{i},\mathbf{x}^{i}_{k-1})$ from last iteration
without the extra computation as in SARAH. The difference between GT and SARAH
is not the focus of this work and instead we will combine both Hybrid-SARAH
and GT together to reduce the stochastic gradient variance and correct the
stochastic gradient bias in the MARL.

### Additional Analysis for Algorithm 1

We discuss briefly how MDPGT enables to reduce the policy gradient variance.
Let $\Delta^{i}_{k}\triangleq\mathbf{u}^{i}_{k}-\nabla
J_{i}(\mathbf{x}^{i}_{k})$. With simple mathematical manipulations, we can
obtain

$\begin{split}&\mathbb{E}[\Delta^{i}_{k}]=\mathbb{E}[(1-\beta)\Delta^{i}_{k-1}+\beta(\underbrace{\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})-\nabla
J_{i}(\mathbf{x}^{i}_{k})}_{T_{3}})+\\\
&(1-\beta)(\underbrace{\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})-\upsilon_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{x}^{i}_{k})\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1})-\nabla
J_{i}(\mathbf{x}^{i}_{k})+\nabla
J_{i}(\mathbf{x}^{i}_{k-1})}_{T_{4}})].\end{split}$ (15)

We can observe from the last equation that
$\mathbb{E}[\Delta^{i}_{k}]=(1-\beta)\mathbb{E}[\Delta^{i}_{k-1}]$ due to
$\mathbb{E}_{\tau^{i}_{k}\sim p_{i}(\tau^{i}|\mathbf{x}^{i}_{k})}[T_{3}]=0$
and $\mathbb{E}_{\tau^{i}_{k}\sim
p_{i}(\tau^{i}|\mathbf{x}^{i}_{k})}[T_{4}]=0$. Applying Cauthy-Schwartz
inequality, we further attain the upper bound for
$\mathbb{E}[\|\Delta^{i}_{k}\|^{2}]$ as follows,

$\begin{split}\mathbb{E}[\|\Delta^{i}_{k}\|^{2}]&\leq(1-\beta)^{2}\mathbb{E}[\|\Delta^{i}_{k-1}\|^{2}]+2\beta^{2}\mathbb{E}[\|T_{3}\|^{2}]\\\
&+2(1-\beta)^{2}\mathbb{E}[\|T_{4}\|^{2}].\end{split}$ (16)

In the sequel, we will establish the relationship between the importance
sampling weight
$\upsilon_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{x}^{i}_{k})$ and
$\|\mathbf{x}^{i}_{k}-\mathbf{x}^{i}_{k-1}\|^{2}$ such that the following
equality can be obtained

$\begin{split}\mathcal{O}(\|T_{4}\|^{2})&=\mathcal{O}(\|\mathbf{x}^{i}_{k}-\mathbf{x}^{i}_{k-1}\|^{2})\\\
&=\mathcal{O}(\|\mathbf{x}^{i}_{k}+\bar{\mathbf{x}}_{k}-\bar{\mathbf{x}}_{k}+\bar{\mathbf{x}}_{k-1}-\bar{\mathbf{x}}_{k-1}-\mathbf{x}^{i}_{k-1}\|^{2}).\end{split}$
(17)

Eq. 17 is bounded above by the following

$\mathcal{O}(3\|\mathbf{x}^{i}_{k}-\bar{\mathbf{x}}_{k}\|^{2}+3\|\mathbf{x}^{i}_{k-1}-\bar{\mathbf{x}}_{k-1}\|^{2}+3\eta^{2}\|\bar{\mathbf{u}}_{k-1}\|^{2}),$

where $\bar{*}$ is the ensemble average of all agents $i\in\mathcal{V}$ and
the first two terms signify the consensus estimate, which can be controlled by
$\beta$ and $\eta$. The bound also follows from the fact that
$\bar{\mathbf{u}}_{k-1}=\bar{\mathbf{v}}_{k}$ (see Lemma 4). Therefore, in
light of Eq. 16, setting $\beta$ and $\eta$ properly is able to reduce the
variance. This also motivates the adoption of policy gradient surrogate in the
work.

### Auxiliary Lemmas and Their Proof

We start with the results for the smoothness constant $L$ and the upper bound
of $\|\mathbf{g}_{i}(\tau^{i}|\mathbf{x}^{i})\|$, $G$.

######  Lemma 1.

Let $\mathbf{g}_{i}(\tau^{i}|\mathbf{x}^{i})$ be defined in Eq. 6 for all
$i\in\mathcal{V}$. Assumption 1 implies the following conclusions:

  * •

$\|\mathbf{g}_{i}(\tau^{i}|\mathbf{x}_{1})-\mathbf{g}_{i}(\tau^{i}|\mathbf{x}_{2})\|\leq
L\|\mathbf{x}_{1}-\mathbf{x}_{2}\|,\forall\mathbf{x}_{1},\mathbf{x}_{2}\in\mathbb{R}^{d_{i}}$,
with $L=\frac{C_{h}R}{(1-\gamma)^{2}}$;

  * •

$J({\mathbf{x}})$ is $L$-smooth, i.e., $\|\nabla^{2}J(\mathbf{x})\|\leq L$;

  * •

$\|\mathbf{g}_{i}(\tau^{i}|\mathbf{x}^{i})\|\leq
G,\forall\mathbf{x}^{i}\in\mathbb{R}^{d_{i}}$, with
$G=\frac{C_{g}R}{(1-\gamma)^{2}}$;

  * •

Variance of stochastic policy gradient
$\mathbf{g}_{i}(\tau^{i}|\mathbf{x}^{i})$ is bounded, i.e.,
$\mathbb{E}[\|\mathbf{g}_{i}(\tau^{i}|\mathbf{x}^{i})-\nabla
J_{i}(\mathbf{x}^{i})\|^{2}]\leq\sigma^{2}_{i},\bar{\sigma}^{2}=\frac{1}{N}\sum^{N}_{i=1}\sigma^{2}_{i}$,
for all $i\in\mathcal{V}$, where
$\bar{\sigma}^{2}=\frac{C_{g}^{2}R^{2}}{(1-\gamma)^{4}}$.

###### Proof.

Recall the definition of $\mathbf{g}_{i}(\tau^{i}|\mathbf{x}^{i})$ in Eq. 6
such that the PGT estimator (Sutton et al. 1999) is as follows:

$\mathbf{g}_{i}(\tau^{i}|\mathbf{x}^{i})=\sum^{H-1}_{h=0}\sum^{h}_{q=0}(\gamma^{q}r^{i}_{q}(a^{i}_{q},s_{q})-b_{q})\nabla_{\mathbf{x}^{i}}\textnormal{log}\pi^{i}(a^{i}_{h},s_{h}),$

where $b_{q}$ is a constant baseline and we specify the subscript
$\mathbf{x}^{i}$ to indicate the gradient w.r.t it. Later this is omitted
unless specified appropriately. Further, we have

$\begin{split}\|\nabla\mathbf{g}_{i}(\tau^{i}|\mathbf{x}^{i})\|&=\bigg{\|}\sum^{H-1}_{h=0}\nabla^{2}\textnormal{log}\pi^{i}(a^{i}_{h},s_{h})\bigg{(}\sum^{h}_{q=0}\gamma^{q}r^{i}_{q}(a^{i}_{q},s_{q})\bigg{)}\bigg{\|}\\\
&\leq\bigg{(}\sum^{H-1}_{q=0}\|\nabla^{2}\textnormal{log}\pi^{i}(a^{i}_{q},s_{q})\|\frac{R}{1-\gamma}\bigg{)}\\\
&\leq\frac{C_{h}R}{(1-\gamma)^{2}}\end{split}$ (18)

where $b_{q}=0$. The above inequalities follow from that $\gamma\in(0,1)$.
When $b_{q}\neq 0$, we can easily scale it with $\gamma^{h}$ and the above
result still holds but with different constant coefficient. Thus, the first
part of Lemma 1 can be proved.

As the PGT estimator is an unbiased estimator of the policy gradient $\nabla
J_{i}(\mathbf{x}^{i})$, we then have

$\nabla
J_{i}(\mathbf{x}^{i})=\mathbb{E}_{\tau^{i}}[\mathbf{g}_{i}(\tau^{i}|\mathbf{x}^{i})],\nabla^{2}J_{i}(\mathbf{x}^{i})=\mathbb{E}_{\tau^{i}}[\nabla\mathbf{g}_{i}(\tau^{i}|\mathbf{x}^{i})].$

Hence, the smoothness of $J_{i}(\mathbf{x}^{i})$ can be directly implied from
the Lipchitzness of $\mathbf{g}_{i}(\tau^{i}|\mathbf{x}^{i})$,

$\|\nabla^{2}J_{i}(\mathbf{x}^{i})\|=\|\mathbb{E}_{\tau^{i}}[\nabla\mathbf{g}_{i}(\tau^{i}|\mathbf{x}^{i})]\|\leq\|\nabla\mathbf{g}_{i}(\tau^{i}|\mathbf{x}^{i})\|\leq\frac{C_{h}R}{(1-\gamma)^{2}}.$

As the above relationship holds for each agent $i\in\mathcal{V}$, the second
part of Lemma 1 can be obtained.

Similarly, taking the norm of $\mathbf{g}_{i}(\tau^{i}|\mathbf{x}^{i})$ leads
to

$\|\mathbf{g}_{i}(\tau^{i}|\mathbf{x}^{i})\|\leq\bigg{\|}\sum^{H-1}_{h=0}\nabla\textnormal{log}\pi^{i}(a^{i}_{h},s_{h})\frac{\gamma^{h}R(1-\gamma^{H-h})}{1-\gamma}\bigg{\|}\leq\frac{C_{g}R}{(1-\gamma)^{2}}.$

Based on the above inequality, it is immediately obtained that

$\|\mathbf{g}_{i}(\tau^{i}|\mathbf{x}^{i})\|^{2}\leq\frac{C_{g}^{2}R^{2}}{(1-\gamma)^{4}}.$

For a random variable $X$, we have that
$\mathbb{E}[\|X-\mathbb{E}[X]\|^{2}]\leq\mathbb{E}[\|X\|^{2}]$ such that

$\mathbb{E}[\|\mathbf{g}_{i}(\tau^{i}|\mathbf{x}^{i})-\nabla
J_{i}(\mathbf{x}^{i})\|^{2}]\leq\frac{C_{g}^{2}R^{2}}{(1-\gamma)^{4}}.$

Also, $\bar{\sigma}^{2}=\frac{C_{g}^{2}R^{2}}{(1-\gamma)^{4}}$. The proof is
complete. ∎

Lemma 1 completes the definitions of $L$ and $G$. One immediate observation
for both $L$ and $G$ is that they are independent of the horizon $H$, which
hence are tighter. In the following, we introduce a relation for the
importance sampling weight. We first introduce the Rényi divergence between
two distributions $Q$ and $Z$ as follows:

$\mathcal{D}_{\alpha}(Q||Z)=\frac{1}{\alpha-1}\textnormal{log}_{2}\int_{x}Q(x)\bigg{(}\frac{Q(x)}{Z(x)}\bigg{)}^{\alpha-1}dx,$

which is non-negative for all $\alpha>0$. The exponentiated Rényi divergence
is $\mathcal{D}_{\alpha}(Q||Z)=2^{D_{\alpha}(Q||Z)}$.

######  Lemma 2.

(Lemma 1 in (Cortes, Mansour, and Mohri 2010)) Let
$\upsilon_{i}(\tau^{i}|\mathbf{x}_{1},\mathbf{x}_{2})=\frac{p_{i}(\tau^{i}|\mathbf{x}_{1})}{p_{i}(\tau^{i}|\mathbf{x}_{2})}$
be the importance sampling weight for distributions
$p_{i}(\tau^{i}|\mathbf{x}_{1})$ and $p_{i}(\tau^{i}|\mathbf{x}_{2})$. Then
the expectation, second moment and variance of
$\upsilon_{i}(\tau^{i}|\mathbf{x}_{1},\mathbf{x}_{2})$ satisfy the following
results:

$\mathbb{E}[\upsilon_{i}(\tau^{i}|\mathbf{x}_{1},\mathbf{x}_{2})]=1,\mathbb{E}[\upsilon_{i}^{2}(\tau^{i}|\mathbf{x}_{1},\mathbf{x}_{2})]=\mathcal{D}_{2}(p_{i}(\tau^{i}|\mathbf{x}_{1})||p_{i}(\tau^{i}|\mathbf{x}_{2})),$

and

$\mathbb{V}(\upsilon_{i}(\tau^{i}|\mathbf{x}_{1},\mathbf{x}_{2}))=\mathcal{D}_{2}(p_{i}(\tau^{i}|\mathbf{x}_{1})||p_{i}(\tau^{i}|\mathbf{x}_{2}))-1.$

With Lemma 2 in hand, we now are ready to state the relation between
$\mathbb{V}(\upsilon_{i}(\tau^{i}|\mathbf{x}_{1},\mathbf{x}_{2}))$ and
$\|\mathbf{x}_{1}-\mathbf{x}_{2}\|^{2}$.

######  Lemma 3.

For any $\mathbf{x}_{1},\mathbf{x}_{2}\in\mathbb{R}^{d_{i}},\forall
i\in\mathcal{V}$, let
$\upsilon_{i}(\tau^{i}|\mathbf{x}_{1},\mathbf{x}_{2})=\frac{p_{i}(\tau^{i}|\mathbf{x}_{1})}{p_{i}(\tau^{i}|\mathbf{x}_{2})}$.
With Assumptions 1 and 3, We have

$\mathbb{V}[\upsilon_{i}(\tau^{i}|\mathbf{x}_{1},\mathbf{x}_{2})]\leq
C_{\upsilon}\|\mathbf{x}_{1}-\mathbf{x}_{2}\|^{2},$ (19)

where $C_{\upsilon}=H(2HC_{g}^{2}+C_{h})(\mathcal{M}+1)$.

###### Proof.

According to Lemma 2, the variance of the importance sampling weight is

$\mathbb{V}[\upsilon_{i}(\tau^{i}|\mathbf{x}_{1},\mathbf{x}_{2})]=\mathcal{D}_{2}(p_{i}(\tau^{i}|\mathbf{x}_{1})||p_{i}(\tau^{i}|\mathbf{x}_{2}))-1.$

By definition, the following relationship can be obtained

$\mathcal{D}_{2}(p_{i}(\tau^{i}|\mathbf{x}_{1})||p_{i}(\tau^{i}|\mathbf{x}_{2}))=\int_{\tau^{i}}p_{i}(\tau^{i}|\mathbf{x}_{1})\frac{p_{i}(\tau^{i}|\mathbf{x}_{1})}{p_{i}(\tau^{i}|\mathbf{x}_{2})}d\tau^{i}=\int_{\tau^{i}}p_{i}(\tau^{i}|\mathbf{x}_{1})^{2}p_{i}(\tau^{i}|\mathbf{x}_{2})^{-1}d\tau^{i}.$

Taking the gradient of
$\mathcal{D}_{2}(p_{i}(\tau^{i}|\mathbf{x}_{1})||p_{i}(\tau^{i}|\mathbf{x}_{2}))$
w.r.t. $\mathbf{x}_{1}$, we then have

$\nabla_{\mathbf{x}_{1}}\mathcal{D}_{2}(p_{i}(\tau^{i}|\mathbf{x}_{1})||p_{i}(\tau^{i}|\mathbf{x}_{2}))=2\int_{\tau^{i}}p_{i}(\tau^{i}|\mathbf{x}_{1})\nabla_{\mathbf{x}_{1}}p_{i}(\tau^{i}|\mathbf{x}_{1})p_{i}(\tau^{i}|\mathbf{x}_{2})^{-1}d\tau^{i}.$

Letting $\mathbf{x}_{1}=\mathbf{x}_{2}$ in the last equation yields

$\nabla_{\mathbf{x}_{1}}\mathcal{D}_{2}(p_{i}(\tau^{i}|\mathbf{x}_{1})||p_{i}(\tau^{i}|\mathbf{x}_{2}))\bigg{|}_{\mathbf{x}_{1}=\mathbf{x}_{2}}=2\int_{\tau^{i}}\nabla_{\mathbf{x}_{1}}p_{i}(\tau^{i}|\mathbf{x}_{1})d\tau^{i}\bigg{|}_{\mathbf{x}_{1}=\mathbf{x}_{2}}=0.$

Applying the mean value theorem w.r.t. $\mathbf{x}_{1}$ results in

$\mathcal{D}_{2}(p_{i}(\tau^{i}|\mathbf{x}_{1})||p_{i}(\tau^{i}|\mathbf{x}_{2}))=1+\frac{1}{2}(\mathbf{x}_{1}-\mathbf{x}_{2})^{\top}\nabla^{2}_{\mathbf{x}_{3}}\mathcal{D}_{2}(p_{i}(\tau^{i}|\mathbf{x}_{3})||p_{i}(\tau^{i}|\mathbf{x}_{2}))(\mathbf{x}_{1}-\mathbf{x}_{2}),$
(20)

where $\mathbf{x}_{3}=t\mathbf{x}_{1}+(1-t)\mathbf{x}_{2},t\in[0,1]$. The last
inequality follows from the fact that
$\mathcal{D}_{2}(p_{i}(\tau^{i}|\mathbf{x}_{2})||p_{i}(\tau^{i}|\mathbf{x}_{2}))=1$.
We have now obtained another expression for
$\mathcal{D}_{2}(p_{i}(\tau^{i}|\mathbf{x}_{1})||p_{i}(\tau^{i}|\mathbf{x}_{2}))$
and to bound it, we shall compute the Hessian matrix
$\nabla^{2}_{\mathbf{x}_{3}}\mathcal{D}_{2}(p_{i}(\tau^{i}|\mathbf{x}_{3})||p_{i}(\tau^{i}|\mathbf{x}_{2}))$.
Taking the gradient
$\nabla_{\mathbf{x}_{3}}\mathcal{D}_{2}(p_{i}(\tau^{i}|\mathbf{x}_{3})||p_{i}(\tau^{i}|\mathbf{x}_{2}))$
of w.r.t. $\mathbf{x}_{3}$ leads to

$\begin{split}\nabla^{2}_{\mathbf{x}_{3}}\mathcal{D}_{2}(p_{i}(\tau^{i}|\mathbf{x}_{3})||p_{i}(\tau^{i}|\mathbf{x}_{2}))&=2\int_{\tau^{i}}\nabla_{\mathbf{x}_{3}}\textnormal{log}p_{i}(\tau^{i}|\mathbf{x}_{3})\nabla_{\mathbf{x}_{3}}\textnormal{log}p_{i}(\tau^{i}|\mathbf{x}_{3})^{\top}\frac{p_{i}(\tau^{i}|\mathbf{x}_{3})^{2}}{p_{i}(\tau^{i}|\mathbf{x}_{2})}d\tau^{i}\\\
&+2\int_{\tau^{i}}\nabla_{\mathbf{x}_{3}}^{2}p_{i}(\tau^{i}|\mathbf{x}_{3})p_{i}(\tau^{i}|\mathbf{x}_{3})p_{i}(\tau^{i}|\mathbf{x}_{2})^{-1}d\tau^{i}.\end{split}$

The above equation implies that to obtain the Hessian matrix, we need to
attain $\nabla_{\mathbf{x}_{3}}^{2}p_{i}(\tau^{i}|\mathbf{x}_{3})$, which
signifies the Hessian matrix of the trajectory distribution function. We next
derive the the Hessian matrix of log-density function.

$\nabla_{\mathbf{x}_{3}}^{2}\textnormal{log}p_{i}(\tau^{i}|\mathbf{x}_{3})=-p_{i}(\tau^{i}|\mathbf{x}_{3})^{2}\nabla_{\mathbf{x}_{3}}p_{i}(\tau^{i}|\mathbf{x}_{3})\nabla_{\mathbf{x}_{3}}p_{i}(\tau^{i}|\mathbf{x}_{3})^{\top}+p_{i}(\tau^{i}|\mathbf{x}_{3})^{-1}\nabla_{\mathbf{x}_{3}}^{2}p_{i}(\tau^{i}|\mathbf{x}_{3}).$

Combining the last two equations yields

$\begin{split}\|\nabla^{2}_{\mathbf{x}_{3}}\mathcal{D}_{2}(p_{i}(\tau^{i}|\mathbf{x}_{3})||p_{i}(\tau^{i}|\mathbf{x}_{2}))\|&=\bigg{\|}4\int_{\tau^{i}}\nabla_{\mathbf{x}_{3}}\textnormal{log}p_{i}(\tau^{i}|\mathbf{x}_{3})\nabla_{\mathbf{x}_{3}}p_{i}(\tau^{i}|\mathbf{x}_{3})^{\top}\frac{p_{i}(\tau^{i}|\mathbf{x}_{3})^{2}}{p_{i}(\tau^{i}|\mathbf{x}_{2})}d\tau^{i}\\\
&+2\int_{\tau^{i}}\nabla^{2}_{\mathbf{x}_{3}}\textnormal{log}p_{i}(\tau^{i}|\mathbf{x}_{3})\frac{p_{i}(\tau^{i}|\mathbf{x}_{3})^{2}}{p_{i}(\tau^{i}|\mathbf{x}_{2})}d\tau^{i}\bigg{\|}\\\
&\leq\int_{\tau^{i}}\frac{p_{i}(\tau^{i}|\mathbf{x}_{3})^{2}}{p_{i}(\tau^{i}|\mathbf{x}_{2})}(4\|\nabla_{\mathbf{x}_{3}}\textnormal{log}p_{i}(\tau^{i}|\mathbf{x}_{3})\|^{2}+2\|\nabla^{2}_{\mathbf{x}_{3}}\textnormal{log}p_{i}(\tau^{i}|\mathbf{x}_{3})\|)d\tau^{i}\\\
&\leq(4H^{2}C^{2}_{g}+2HC_{h})\mathbb{E}[\upsilon_{i}(\tau^{i}|\mathbf{x}_{1},\mathbf{x}_{2})^{2}]\\\
&\leq 2H(2HC^{2}_{g}+C_{h})(\mathcal{M}+1),\end{split}$

where the second inequality follows from Assumption 1 and Lemma 2 and the last
inequality follows from Assumption 3. Therefore, substituting the above result
into Eq. 20 attain the following desirable result

$\mathbb{V}[\upsilon_{i}(\tau^{i}|\mathbf{x}_{1},\mathbf{x}_{2})]=\mathcal{D}_{2}(p_{i}(\tau^{i}|\mathbf{x}_{1})||p_{i}(\tau^{i}|\mathbf{x}_{2}))-1\leq
C_{\upsilon}\|\mathbf{x}_{1}-\mathbf{x}_{2}\|^{2}.$

∎

Lemma 3 is critical in the proof as it will translate HIS-SARAH into a squared
norm that can be bounded above. We next present the outline of how to prove
Theorem 1 and apply the same proof techniques to show Theorem 2. We first
define some notations for the convenience of proof. In a generalized case, the
dimension for each $\mathbf{x}^{i}$ is not necessarily the same, while in this
work we assume that $d_{1}=d_{2}=...=d_{N}=\frac{d}{N}$ for the ease of
exposition and in this context, $d/N$ is assumed to be an integer. Recalling
the update laws for $\mathbf{v}^{i}_{k}$ and $\mathbf{x}^{i}_{k}$ in Algorithm
1 in a matrix form, we have:

$\displaystyle\mathbf{v}_{k+1}$
$\displaystyle=\underline{\mathbf{W}}\mathbf{v}_{k}+\mathbf{u}_{k}-\mathbf{u}_{k-1}$
(21a) $\displaystyle\mathbf{x}_{k+1}$
$\displaystyle=\underline{\mathbf{W}}(\mathbf{x}_{k}+\eta\mathbf{v}_{k+1})$
(21b)

where $\underline{\mathbf{W}}\triangleq\mathbf{W}\otimes\mathbf{I}_{d/N}$ and
$\mathbf{x}_{k},\mathbf{v}_{k},\mathbf{u}_{k}$ are square-integrable random
vectors in $\mathbb{R}^{d}$ that concatenate the local estimates of the
solution $\\{\mathbf{x}^{i}_{k}\\}_{i=1}^{N}$, gradient trackers
$\\{\mathbf{v}^{i}_{k}\\}_{i=1}^{N}$, and the stochastic policy gradient
surrogates $\\{\mathbf{u}^{i}_{k}\\}^{N}_{i=1}$. Additionally, we denote
$\nabla\mathbf{J}(\mathbf{x}_{k})=[\nabla
J_{1}(\mathbf{x}_{k}^{1})^{\top},...,\nabla
J_{N}(\mathbf{x}_{k}^{N})^{\top}]^{\top}$ and define the ensemble averages as
follows.

$\bar{\mathbf{x}}_{k}\triangleq\frac{1}{N}(\mathbf{1}^{\top}_{N}\otimes\mathbf{I}_{d/N})\mathbf{x}_{k},\bar{\mathbf{v}}_{k}\triangleq\frac{1}{N}(\mathbf{1}^{\top}_{N}\otimes\mathbf{I}_{d/N})\mathbf{v}_{k},$
$\bar{\mathbf{u}}_{k}\triangleq\frac{1}{N}(\mathbf{1}^{\top}_{N}\otimes\mathbf{I}_{d/N})\mathbf{u}_{k},\overline{\nabla\mathbf{J}}(\mathbf{x}_{k})\triangleq\frac{1}{N}(\mathbf{1}^{\top}_{N}\otimes\mathbf{I}_{d/N})\nabla\mathbf{f}(\mathbf{x}_{k}).$

With these definitions, we present some known results for gradient tracking
type of algorithms in decentralized optimization and refer interested readers
to (Xin, Khan, and Kar 2021) for the detail of proof.

######  Lemma 4.

The following relationships hold for MDPGT.

  * •

Define
$\mathbf{\Lambda}=\frac{1}{N}(\mathbf{1}^{\top}_{N}\otimes\mathbf{I}_{d/N})$.
Thus,
$\|\underline{\mathbf{W}}\mathbf{x}-\mathbf{\Lambda}\mathbf{x}\|\leq\lambda\|\mathbf{x}-\mathbf{\Lambda}\mathbf{x}\|$.

  * •

$\bar{\mathbf{v}}_{k+1}=\bar{\mathbf{u}}_{k}$, for all $k\geq 0$.

  * •

$\|\overline{\nabla\mathbf{J}}(\mathbf{x}_{k})-\nabla
J(\bar{\mathbf{x}}_{k})\|^{2}\leq\frac{L^{2}}{N}\|\mathbf{x}_{k}-\mathbf{\Lambda}\mathbf{x}_{k}\|^{2}$,
for all $k\geq 0$.

The first relationship in Lemma 4 has been well-known due to a doubly
stochastic mixing matrix $\underline{\mathbf{W}}$ and the second one is an
immediate consequence of the update law for the gradient tracker. The third
relationship is due to the second conclusion from Lemma 1. Multiplying the
update for $\mathbf{x}^{i}_{k+1}$ by $\Lambda$ yields the following equality

$\bar{\mathbf{x}}_{k+1}=\bar{\mathbf{x}}_{k}+\eta\bar{\mathbf{v}}_{k+1}=\bar{\mathbf{x}}_{k}+\eta\bar{\mathbf{u}}_{k},\;\forall
k\geq 0.$ (22)

We next establish a key technical lemma that sheds light on the convergence in
terms of the second moment of $\nabla J(\bar{\mathbf{x}}_{k})$.

######  Lemma 5.

Let $\bar{\mathbf{x}}_{k}$ be generated by Eq. 22. If the step size
$\eta\in(0,\frac{1}{2L}]$, then for all $K\geq 0$, we have:

$\begin{split}&\sum_{k=0}^{K}\|\nabla
J(\bar{\mathbf{x}}_{k})\|^{2}\leq\frac{2(J^{*}-J(\bar{\mathbf{x}}_{0}))}{\eta}-\frac{1}{2}\sum_{k=0}^{K}\|\bar{\mathbf{u}}_{k}\|^{2}\\\
&+2\sum_{k=0}^{K}\|\bar{\mathbf{u}}_{k}-\overline{\nabla\mathbf{J}}(\mathbf{x}_{k})\|^{2}+\frac{2L^{2}}{N}\sum_{k=0}^{K}\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}.\end{split}$
(23)

###### Proof.

Based on the second conclusion of Lemma 1, we have the following relationship:

$\begin{split}J(\bar{\mathbf{x}}_{k+1})&\geq
J(\bar{\mathbf{x}}_{k})+\langle\nabla
J(\bar{\mathbf{x}}_{k}),(\bar{\mathbf{x}}_{k+1}-\bar{\mathbf{x}}_{k})\rangle-\frac{L}{2}\|\bar{\mathbf{x}}_{k+1}-\bar{\mathbf{x}}_{k}\|^{2}\\\
&\geq J(\bar{\mathbf{x}}_{k})+\eta\langle\nabla
J(\bar{\mathbf{x}}_{k}),\bar{\mathbf{u}}_{k}\rangle-\frac{\eta^{2}L}{2}\|\bar{\mathbf{u}}_{k}\|^{2}.\end{split}$
(24)

According to the basic inequality
$\langle\mathbf{a},\mathbf{b}\rangle=\frac{1}{2}(\|\mathbf{a}\|^{2}+\|\mathbf{b}\|^{2}-\|\mathbf{a}-\mathbf{b}\|^{2})$,
the last inequality, the following can be obtained

$\begin{split}J(\bar{\mathbf{x}}_{k+1})&\geq
J(\bar{\mathbf{x}}_{k})+\frac{\eta}{2}\|\nabla
J(\bar{\mathbf{x}}_{k})\|^{2}+\bigg{(}\frac{\eta}{2}-\frac{\eta^{2}L}{2}\bigg{)}\|\bar{\mathbf{u}}_{k}\|^{2}-\frac{\eta}{2}\|\bar{\mathbf{u}}_{k}-\nabla
J(\bar{\mathbf{x}}_{k})\|^{2}\\\ &\geq
J(\bar{\mathbf{x}}_{k})+\frac{\eta}{2}\|\nabla
J(\bar{\mathbf{x}}_{k})\|^{2}+\bigg{(}\frac{\eta}{2}-\frac{\eta^{2}L}{2}\bigg{)}\|\bar{\mathbf{u}}_{k}\|^{2}-\frac{\eta}{2}\|\bar{\mathbf{u}}_{k}-\nabla
J(\bar{\mathbf{x}}_{k})+\overline{\nabla\mathbf{J}}(\mathbf{x}_{k})-\overline{\nabla\mathbf{J}}(\mathbf{x}_{k})\|^{2}\\\
&\geq J(\bar{\mathbf{x}}_{k})+\frac{\eta}{2}\|\nabla
J(\bar{\mathbf{x}}_{k})\|^{2}+\frac{\eta}{4}\|\bar{\mathbf{u}}_{k}\|^{2}-\eta\|\bar{\mathbf{u}}_{k}-\overline{\nabla\mathbf{J}}(\mathbf{x}_{k})\|^{2}-\frac{\eta
L^{2}}{N}\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2},\end{split}$ (25)

where the third inequality is due to $\|\mathbf{a}+\mathbf{b}\|^{2}\leq
2(\|\mathbf{a}\|^{2}+\|\mathbf{b}\|^{2})$, the smoothness of $J$ and
$\eta\in(0,\frac{1}{2L}]$. Rewriting Eq. 25 yields the following relationship

$\|\nabla
J(\bar{\mathbf{x}}_{k})\|^{2}\leq\frac{2(J(\bar{\mathbf{x}}_{k+1})-J(\bar{\mathbf{x}}_{k}))}{\eta}-\frac{1}{2}\|\bar{\mathbf{u}}_{k}\|^{2}+2\|\bar{\mathbf{u}}_{k}-\overline{\nabla\mathbf{J}}(\mathbf{x}_{k})\|^{2}+\frac{2L^{2}}{N}\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}.$
(26)

Summing up the last equation from $0$ to $K$ completes the proof. ∎

According to Lemma 5, it is clearly observed that in order to arrive at the
main results presented in the previous section, the upper bound of the right
hand side in the last inequality needs to be explicitly quantified. One may
argue that in the main theorems it is $\tilde{\mathbf{x}}_{K}$ instead of
$\bar{\mathbf{x}}_{K}$. Dividing Eq. 23 by $\frac{1}{K+1}$ allows us to
approximate $\tilde{\mathbf{x}}_{K}$ using the ensemble average. Additionally,
a sufficiently large $K$ implies all agents have converged to the (near-)
optimal solution, suggesting that the gap between $\tilde{\mathbf{x}}_{K}$ and
$\bar{\mathbf{x}}_{K}$ is quite close to 0. Though in the final result, the
convergence rate is specifically for $\tilde{\mathbf{x}}_{K}$, Lemma 5 helps
facilitate the understanding of convergence for MDPGT and will be critical to
derive the final error bounds in the main theorems. Thus, in the sequel, the
analysis is dedicated to finding the upper bounds for
$\sum_{k=0}^{K}\|\bar{\mathbf{u}}_{k}-\overline{\nabla\mathbf{J}}(\mathbf{x}_{k})\|^{2}$
and $\sum_{k=0}^{K}\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}$. To start the
analysis, we present two relationships that will help significantly
characterize the recursions of gradient variances, consensus, and gradient
tracking errors.

######  Lemma 6.

Let $\\{e_{k}\\}_{k\geq 0},\\{b_{k}\\}_{k\geq 0}$, and $\\{c_{k}\\}_{k\geq 0}$
be nonnegative sequences such that $e_{k}\leq
qe_{k-1}+qb_{k-1}+c_{k}+C,\forall k\geq 1$, where $q\in(0,1)$ and $C\geq 0$.
Then, we have for all $K\geq 1$,

$\sum_{k=1}^{K}e_{k}\leq\frac{e_{0}}{1-q}+\frac{1}{1-q}\sum^{K-1}_{k=0}b_{k}+\frac{1}{1-q}\sum^{K}_{k=0}c_{k}+\frac{CK}{1-q}.$
(27)

Similarly, if $e_{k+1}\leq qe_{k}+b_{k-1}+C,\forall k\geq 1$, then, we have
for all $K\geq 2$,

$\sum_{k=1}^{K}e_{k}\leq\frac{e_{1}}{1-q}+\frac{1}{1-q}\sum^{K-2}_{k=0}b_{k}+\frac{CK}{1-q}.$
(28)

Please follow the proof of Lemma 6 in (Xin, Khan, and Kar 2021) for more
detail. We next derive the relevant recursions for the gradient variances
$\mathbb{E}[\|\bar{\mathbf{u}}_{k}-\overline{\nabla\mathbf{J}}(\mathbf{x}_{k})\|^{2}]$
and $\mathbb{E}[\|\mathbf{u}_{k}-\nabla\mathbf{J}(\mathbf{x}_{k})\|^{2}]$.
Essentially the former can be treated intuitively as the average of the latter
up to a factor of $\frac{1}{N^{2}}$, but they are proved separately.

######  Lemma 7.

Let $\bar{\mathbf{u}}_{k}$ and $\mathbf{x}_{k}$ be generated by MDPGT. Then
for all $k\geq 1$, we have

$\begin{split}&\mathbb{E}[\|\bar{\mathbf{u}}_{k}-\overline{\nabla\mathbf{J}}(\mathbf{x}_{k})\|^{2}]\leq(1-\beta)^{2}\mathbb{E}[\|\bar{\mathbf{u}}_{k-1}-\overline{\nabla\mathbf{J}}(\mathbf{x}_{k-1})\|^{2}]\\\
&+\frac{12(L^{2}+G^{2}C_{\upsilon})\eta^{2}(1-\beta)^{2}}{N}\mathbb{E}[\|\bar{\mathbf{u}}_{k-1}\|^{2}]+\frac{2\beta^{2}\bar{\sigma}^{2}}{N}\\\
&+\frac{12(L^{2}+G^{2}C_{\upsilon})(1-\beta)^{2}}{N^{2}}\mathbb{E}[\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}\\\
&+\|\mathbf{x}_{k-1}-\Lambda\mathbf{x}_{k-1}\|^{2}],\end{split}$ (29)

and

$\begin{split}&\mathbb{E}[\|\mathbf{u}_{k}-\nabla\mathbf{J}(\mathbf{x}_{k})\|^{2}]\leq(1-\beta)^{2}\mathbb{E}[\|\mathbf{u}_{k-1}-\nabla\mathbf{J}(\mathbf{x}_{k-1})\|^{2}]\\\
&+12N(L^{2}+G^{2}C_{\upsilon})\eta^{2}(1-\beta)^{2}\mathbb{E}[\|\bar{\mathbf{u}}_{k-1}\|^{2}]+2\beta^{2}\bar{\sigma}^{2}N\\\
&+12(L^{2}+G^{2}C_{\upsilon})(1-\beta)^{2}\mathbb{E}[\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}\\\
&+\|\mathbf{x}_{k-1}-\Lambda\mathbf{x}_{k-1}\|^{2}].\end{split}$ (30)

###### Proof.

Recalling

$\mathbf{u}^{i}_{k}=\beta\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})+(1-\beta)[\mathbf{u}_{k-1}^{i}+\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})-\upsilon_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{x}^{i}_{k})\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1})],$
(31)

we then have

$\begin{split}\mathbf{u}^{i}_{k}-\nabla
J_{i}(\mathbf{x}^{i}_{k})&=\beta\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})+(1-\beta)[\mathbf{u}_{k-1}^{i}+\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})-\upsilon_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{x}^{i}_{k})\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1})]\\\
&-\beta\nabla J_{i}(\mathbf{x}^{i}_{k})-(1-\beta)\nabla
J_{i}(\mathbf{x}^{i}_{k})\\\
&=\beta(\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})-\nabla
J_{i}(\mathbf{x}^{i}_{k}))+(1-\beta)[\nabla
J_{i}(\mathbf{x}^{i}_{k-1})+\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})\\\
&-\upsilon_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{x}^{i}_{k})\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1})-\nabla
J_{i}(\mathbf{x}^{i}_{k})]+(1-\beta)(\mathbf{u}^{i}_{k}-\nabla
J_{i}(\mathbf{x}^{i}_{k-1})).\end{split}$ (32)

As $\mathbb{E}[\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})-\nabla
J_{i}(\mathbf{x}^{i}_{k})]=\mathbf{0}_{d/N}$ and
$\mathbb{E}[\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})-\upsilon_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{x}^{i}_{k})\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1})+\nabla
J_{i}(\mathbf{x}^{i}_{k-1})-\nabla
J_{i}(\mathbf{x}^{i}_{k})]=\mathbf{0}_{d/N}$. Applying Eq. 32 from 1 to $N$
and taking the ensemble average results in

$\begin{split}\bar{\mathbf{u}}_{k}-\overline{\nabla\mathbf{J}}(\mathbf{x}_{k})&=(1-\beta)(\bar{\mathbf{u}}_{k-1}-\overline{\nabla\mathbf{J}}(\mathbf{x}_{k-1}))+\beta\frac{1}{N}\sum^{N}_{i=1}(\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})-\nabla
J_{i}(\mathbf{x}^{i}_{k}))\\\
&+(1-\beta)\frac{1}{N}\sum^{N}_{i=1}(\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})-\upsilon_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{x}^{i}_{k})\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1})+\nabla
J_{i}(\mathbf{x}^{i}_{k-1})-\nabla J_{i}(\mathbf{x}^{i}_{k}))\end{split}$ (33)

Let

$\mathbf{n}_{k}=\frac{1}{N}\sum^{N}_{i=1}(\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})-\nabla
J_{i}(\mathbf{x}^{i}_{k}))$

and

$\mathbf{z}_{k}=\frac{1}{N}\sum^{N}_{i=1}(\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})-\upsilon_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{x}^{i}_{k})\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1})+\nabla
J_{i}(\mathbf{x}^{i}_{k-1})-\nabla J_{i}(\mathbf{x}^{i}_{k})).$

We have

$\begin{split}\mathbb{E}[\|\bar{\mathbf{u}}_{k}-\overline{\nabla\mathbf{J}}(\mathbf{x}_{k})\|^{2}]&=(1-\beta)^{2}\|\bar{\mathbf{u}}_{k-1}-\overline{\nabla\mathbf{J}}(\mathbf{x}_{k-1})\|^{2}+\mathbb{E}[\|\beta\mathbf{n}_{k}+(1-\beta)\mathbf{z}_{k}\|^{2}]\\\
&+2\mathbb{E}[\langle(1-\beta)(\bar{\mathbf{u}}_{k-1}-\overline{\nabla\mathbf{J}}(\mathbf{x}_{k-1})),\beta\mathbf{n}_{k}+(1-\beta)\mathbf{z}_{k}\rangle]\\\
&=(1-\beta)^{2}\|\bar{\mathbf{u}}_{k-1}-\overline{\nabla\mathbf{J}}(\mathbf{x}_{k-1})\|^{2}+\mathbb{E}[\|\beta\mathbf{n}_{k}+(1-\beta)\mathbf{z}_{k}\|^{2}]\\\
&\leq(1-\beta)^{2}\|\bar{\mathbf{u}}_{k-1}-\overline{\nabla\mathbf{J}}(\mathbf{x}_{k-1})\|^{2}+2\beta^{2}\mathbb{E}[\|\mathbf{n}_{k}\|^{2}]+2(1-\beta)^{2}\mathbb{E}[\|\mathbf{z}_{k}\|^{2}].\end{split}$
(34)

The second equality follows from that
$\mathbb{E}[\langle(1-\beta)(\bar{\mathbf{u}}_{k-1}-\overline{\nabla\mathbf{J}}(\mathbf{x}_{k-1})),\beta\mathbf{n}_{k}+(1-\beta)\mathbf{z}_{k}\rangle]=0$
as $\mathbb{E}[\mathbf{n}_{k}]=\mathbf{0}_{d/N}$ and
$\mathbb{E}[\mathbf{z}_{k}]=\mathbf{0}_{d/N}$. We next bound the terms
$\mathbb{E}[\|\mathbf{n}_{k}\|^{2}]$ and $\mathbb{E}[\|\mathbf{z}_{k}\|^{2}]$.
As for $k\geq 1$,

$\begin{split}\mathbb{E}[\|\mathbf{n}_{k}\|^{2}]&=\frac{1}{N^{2}}\sum^{N}_{i=1}\mathbb{E}[\|\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})-\nabla
J_{i}(\mathbf{x}^{i}_{k})\|^{2}]+\frac{1}{N^{2}}\sum_{i\neq
j}\mathbb{E}[\langle\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})-\nabla
J_{i}(\mathbf{x}^{i}_{k}),\mathbf{g}_{j}(\tau^{j}_{k}|\mathbf{x}^{j}_{k})-\nabla
J_{j}(\mathbf{x}^{j}_{k})\rangle]\\\
&=\frac{1}{N^{2}}\sum^{N}_{i=1}\mathbb{E}[\|\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})-\nabla
J_{i}(\mathbf{x}^{i}_{k})\|^{2}]\leq\frac{\bar{\sigma}^{2}}{N}.\end{split}$
(35)

The second equality follows from that

$\begin{split}&\mathbb{E}[\langle\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})-\nabla
J_{i}(\mathbf{x}^{i}_{k}),\mathbf{g}_{j}(\tau^{j}_{k}|\mathbf{x}^{j}_{k})-\nabla
J_{j}(\mathbf{x}^{j}_{k})\rangle]\\\
&=\mathbb{E}[\langle\mathbb{E}[\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})]-\nabla
J_{i}(\mathbf{x}^{i}_{k}),\mathbf{g}_{j}(\tau^{j}_{k}|\mathbf{x}^{j}_{k})-\nabla
J_{j}(\mathbf{x}^{j}_{k})\rangle]\\\ &=0.\end{split}$ (36)

For the term $\mathbb{E}[\|\mathbf{z}_{k}\|^{2}]$, we have for all $k\geq 1$,

$\begin{split}\mathbb{E}[\|\mathbf{z}_{k}\|^{2}]&=\mathbb{E}[\|\frac{1}{N}\sum^{N}_{i=1}(\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})-\upsilon_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{x}^{i}_{k})\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1})+\nabla
J_{i}(\mathbf{x}^{i}_{k-1})-\nabla J_{i}(\mathbf{x}^{i}_{k}))\|^{2}]\\\
&=\frac{1}{N^{2}}\sum^{N}_{i=1}\mathbb{E}[\|\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})-\upsilon_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{x}^{i}_{k})\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1})+\nabla
J_{i}(\mathbf{x}^{i}_{k-1})-\nabla J_{i}(\mathbf{x}^{i}_{k})\|^{2}]\\\
&+\frac{1}{N^{2}}\sum_{i\neq
j}\mathbb{E}[\langle\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})-\upsilon_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{x}^{i}_{k})\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1})+\nabla
J_{i}(\mathbf{x}^{i}_{k-1})-\nabla J_{i}(\mathbf{x}^{i}_{k}),\\\
&\mathbf{g}_{j}(\tau^{j}_{k}|\mathbf{x}^{j}_{k})-\upsilon_{j}(\tau^{j}_{k}|\mathbf{x}^{j}_{k-1},\mathbf{x}^{j}_{k})\mathbf{g}_{j}(\tau^{j}_{k}|\mathbf{x}^{j}_{k-1})+\nabla
J_{j}(\mathbf{x}^{j}_{k-1})-\nabla J_{j}(\mathbf{x}^{j}_{k})\rangle]\\\
&=\frac{1}{N^{2}}\sum^{N}_{i=1}\mathbb{E}[\|\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})-\upsilon_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{x}^{i}_{k})\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1})+\nabla
J_{i}(\mathbf{x}^{i}_{k-1})-\nabla J_{i}(\mathbf{x}^{i}_{k})\|^{2}]\\\
&\leq\frac{1}{N^{2}}\sum^{N}_{i=1}\mathbb{E}[\|\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})-\upsilon_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{x}^{i}_{k})\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1})\|^{2}].\end{split}$
(37)

The third equality follows from the same argument we have for the
$\mathbb{E}[\langle\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})-\nabla
J_{i}(\mathbf{x}^{i}_{k}),\mathbf{g}_{j}(\tau^{j}_{k}|\mathbf{x}^{j}_{k})-\nabla
J_{j}(\mathbf{x}^{j}_{k})\rangle]$. The last inequality follows from the fact
that

$\mathbb{E}[\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})-\upsilon_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{x}^{i}_{k})\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1})]=\nabla
J_{i}(\mathbf{x}^{i}_{k})-\nabla J_{i}(\mathbf{x}^{i}_{k-1})$

and the variance decomposition, i.e., for any vector
$\mathbf{a}\in\mathbb{R}^{p}$,

$\mathbb{E}[\|\mathbf{a}-\mathbb{E}[\mathbf{a}]\|^{2}]=\mathbb{E}[\|\mathbf{a}\|^{2}]-\|\mathbb{E}[\mathbf{a}]\|^{2}.$

Hence, we have now

$\begin{split}\mathbb{E}[\|\mathbf{z}_{k}\|^{2}]&\leq\frac{1}{N^{2}}\sum^{N}_{i=1}\mathbb{E}[\|\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})-\upsilon_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{x}^{i}_{k})\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1})\|^{2}]\\\
&=\frac{1}{N^{2}}\sum^{N}_{i=1}\mathbb{E}[\|\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})-\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1})+\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1})-\upsilon_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{x}^{i}_{k})\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1})\|^{2}]\\\
&=\frac{1}{N^{2}}\sum^{N}_{i=1}\mathbb{E}[2\|\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})-\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1})\|^{2}+2\|\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1})-\upsilon_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{x}^{i}_{k})\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1})\|^{2}]\\\
&\leq\frac{2L^{2}}{N^{2}}\sum_{i=1}^{N}\mathbb{E}[\|\mathbf{x}^{i}_{k}-\mathbf{x}^{i}_{k-1}\|^{2}]+\frac{2}{N^{2}}\sum^{N}_{i=1}\mathbb{E}[\|1-\upsilon_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{x}^{i}_{k})\|^{2}\|\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1})\|^{2}]\\\
&\leq\frac{2L^{2}}{N^{2}}\sum_{i=1}^{N}\mathbb{E}[\|\mathbf{x}^{i}_{k}-\mathbf{x}^{i}_{k-1}\|^{2}]+\frac{2G^{2}}{N^{2}}\mathbb{E}[\|1-\upsilon_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{x}^{i}_{k})\|^{2}]\\\
&=\frac{2L^{2}}{N^{2}}\sum_{i=1}^{N}\mathbb{E}[\|\mathbf{x}^{i}_{k}-\mathbf{x}^{i}_{k-1}\|^{2}]+\frac{2G^{2}}{N^{2}}\sum^{N}_{i=1}\mathbb{V}(\upsilon_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{x}^{i}_{k}))\\\
&\leq\frac{2L^{2}}{N^{2}}\sum_{i=1}^{N}\mathbb{E}[\|\mathbf{x}^{i}_{k}-\mathbf{x}^{i}_{k-1}\|^{2}]+\frac{2G^{2}C_{\upsilon}}{N^{2}}\sum^{N}_{i=1}\mathbb{E}[\|\mathbf{x}^{i}_{k}-\mathbf{x}^{i}_{k-1}\|^{2}]\\\
&=\frac{2L^{2}+2G^{2}C_{\upsilon}}{N^{2}}\mathbb{E}[\|\mathbf{x}_{k}-\mathbf{x}_{k-1}\|^{2}]\\\
&=\frac{2L^{2}+2G^{2}C_{\upsilon}}{N^{2}}\mathbb{E}[\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}+\Lambda\mathbf{x}_{k}-\Lambda\mathbf{x}_{k-1}+\Lambda\mathbf{x}_{k-1}-\mathbf{x}_{k-1}\|^{2}]\\\
&\leq\frac{6L^{2}+6G^{2}C_{\upsilon}}{N}\mathbb{E}[\|\bar{\mathbf{u}}_{k-1}\|^{2}]+\frac{6L^{2}+6G^{2}C_{\upsilon}}{N^{2}}\mathbb{E}[\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}+\|\Lambda\mathbf{x}_{k-1}-\mathbf{x}_{k-1}\|^{2}].\end{split}$
(38)

The second inequality follows from the Cauthy-Schwartz inequality. The third
inequality follows from Lemma 1 while the fourth inequality follows from Lemma
3. Hence, the following relationship can be obtained

$\begin{split}\mathbb{E}[\|\bar{\mathbf{u}}_{k}-\overline{\nabla\mathbf{f}}(\mathbf{x}_{k})\|^{2}]&\leq(1-\beta)^{2}\mathbb{E}[\|\bar{\mathbf{u}}_{k-1}-\overline{\nabla\mathbf{f}}(\mathbf{x}_{k-1})\|^{2}]+\frac{12L^{2}+12G^{2}C_{\upsilon}}{N}\eta^{2}(1-\beta)^{2}\mathbb{E}[\|\bar{\mathbf{u}}_{k-1}\|^{2}]\\\
&+\frac{2\beta^{2}\bar{\sigma}^{2}}{N}+\frac{12L^{2}+12G^{2}C_{\upsilon}}{N^{2}}(1-\beta)^{2}\mathbb{E}[\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}+\|\Lambda\mathbf{x}_{k-1}-\mathbf{x}_{k-1}\|^{2}],\end{split}$
(39)

which completes the first part of proof for Lemma 7.

Recalling Eq. 32, we have

$\begin{split}\mathbf{u}^{i}_{k}-\nabla
J_{i}(\mathbf{x}^{i}_{k})&=\beta(\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})-\nabla
J_{i}(\mathbf{x}^{i}_{k}))+(1-\beta)[\nabla
J_{i}(\mathbf{x}^{i}_{k-1})+\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})\\\
&-\upsilon_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{x}^{i}_{k})\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1})-\nabla
J_{i}(\mathbf{x}^{i}_{k})]+(1-\beta)(\mathbf{u}^{i}_{k}-\nabla
J_{i}(\mathbf{x}^{i}_{k-1})).\end{split}$ (40)

Based on the analysis above, it is immediately obtained that the expectations
of the first and second terms are all zero. Following the similar proof in Eq.
34, we have for all $k\geq 1$,

$\begin{split}&\mathbb{E}[\|\mathbf{u}^{i}_{k}-\nabla
J_{i}(\mathbf{x}^{i}_{k})\|^{2}]\leq(1-\beta)^{2}\|\mathbf{u}^{i}_{k-1}-\nabla
J_{i}(\mathbf{x}^{i}_{k-1})\|^{2}+2\beta^{2}\mathbb{E}[\|\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})-\nabla
J_{i}(\mathbf{x}^{i}_{k})\|^{2}]\\\ &+2(1-\beta)^{2}\mathbb{E}[\|\nabla
J_{i}(\mathbf{x}^{i}_{k-1})+\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})-\upsilon_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{x}^{i}_{k})\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1})-\nabla
J_{i}(\mathbf{x}^{i}_{k})\|^{2}]\\\
&\leq(1-\beta)^{2}\|\mathbf{u}^{i}_{k-1}-\nabla
J_{i}(\mathbf{x}^{i}_{k-1})\|^{2}+2\beta^{2}\mathbb{E}[\|\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})-\nabla
J_{i}(\mathbf{x}^{i}_{k})\|^{2}]\\\
&+2(1-\beta)^{2}\mathbb{E}[\|\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k})-\upsilon_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{x}^{i}_{k})\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1})\|^{2}]\end{split}$
(41)

According to Lemma 1 and Lemma 3, the following relationship can be attained

$\begin{split}\mathbb{E}[\|\mathbf{u}^{i}_{k}-\nabla
J_{i}(\mathbf{x}^{i}_{k})\|^{2}]&\leq(1-\beta)^{2}\|\mathbf{u}^{i}_{k-1}-\nabla
J_{i}(\mathbf{x}^{i}_{k-1})\|^{2}+2\beta^{2}\sigma^{2}_{i}+4(1-\beta)^{2}(L^{2}+C_{\upsilon}G^{2})\mathbb{E}[\|\mathbf{x}^{i}_{k}-\mathbf{x}^{i}_{k-1}\|^{2}]\\\
&\leq(1-\beta)^{2}\|\mathbf{u}^{i}_{k-1}-\nabla
J_{i}(\mathbf{x}^{i}_{k-1})\|^{2}+2\beta^{2}\sigma^{2}_{i}\\\
&+12(1-\beta)^{2}(L^{2}+C_{\upsilon}G^{2})(\mathbb{E}[\|\mathbf{x}^{i}_{k}-\bar{\mathbf{x}}_{k}\|^{2}]+\|\bar{\mathbf{x}}_{k}-\bar{\mathbf{x}}_{k-1}\|^{2}+\|\bar{\mathbf{x}}_{k-1}-\mathbf{x}^{i}_{k-1}\|^{2})\\\
&=(1-\beta)^{2}\|\mathbf{u}^{i}_{k-1}-\nabla
J_{i}(\mathbf{x}^{i}_{k-1})\|^{2}+2\beta^{2}\sigma^{2}_{i}+12(1-\beta)^{2}(L^{2}+C_{\upsilon}G^{2})\eta^{2}\mathbb{E}[\|\bar{\mathbf{u}}_{k-1}\|^{2}]\\\
&+12(1-\beta)^{2}(L^{2}+C_{\upsilon}G^{2})(\mathbb{E}[\|\mathbf{x}^{i}_{k}-\bar{\mathbf{x}}_{k}\|^{2}+\|\mathbf{x}^{i}_{k-1}-\bar{\mathbf{x}}_{k-1}\|^{2}]).\end{split}$
(42)

Applying Eq. 42 over $i$ from 1 to $N$ completes the second part of the proof
for Lemma 7. ∎

With Lemma 6 and Lemma 7 in hand, we now are ready to present the upper bounds
for both
$\sum_{k=1}^{K}\mathbb{E}[\|\bar{\mathbf{u}}_{k}-\overline{\nabla\mathbf{J}}(\mathbf{x}_{k})\|^{2}]$
and
$\sum_{k=1}^{K}\mathbb{E}[\|\mathbf{u}_{k}-\nabla\mathbf{J}(\mathbf{x}_{k})\|^{2}]$
in the following lemma.

######  Lemma 8.

Let $\bar{\mathbf{u}}_{k}$ and $\mathbf{x}_{k}$ be generated by MDPGT
initialized with a mini-batch of trajectories $\mathcal{B}$. Then for any
$\beta\in(0,1)$, $\forall K\geq 1$, we have

$\begin{split}&\sum_{k=0}^{K}\mathbb{E}[\|\bar{\mathbf{u}}_{k}-\overline{\nabla\mathbf{J}}(\mathbf{x}_{k})\|^{2}]\leq\frac{\bar{\sigma}^{2}}{|\mathcal{B}|N\beta}+\frac{12(L^{2}+G^{2}C_{\upsilon})\eta^{2}}{N\beta}\times\\\
&\sum_{k=0}^{K-1}\mathbb{E}[\|\bar{\mathbf{u}}_{k}\|^{2}]+\frac{24(L^{2}+G^{2}C_{\upsilon})}{N^{2}\beta}\sum_{k=0}^{K}\mathbb{E}[\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}]\\\
&+\frac{2\beta\bar{\sigma}^{2}K}{N},\end{split}$ (43)

and

$\begin{split}&\sum_{k=0}^{K}\mathbb{E}[\|\mathbf{u}_{k}-\nabla\mathbf{J}(\mathbf{x}_{k})\|^{2}]\leq\frac{N\bar{\sigma}^{2}}{|\mathcal{B}|\beta}+\frac{12N(L^{2}+G^{2}C_{\upsilon})\eta^{2}}{\beta}\times\\\
&\sum_{k=0}^{K-1}\mathbb{E}[\|\bar{\mathbf{u}}_{k}\|^{2}]+\frac{24(L^{2}+G^{2}C_{\upsilon})}{\beta}\sum_{k=0}^{K}\mathbb{E}[\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}]\\\
&+2N\beta\bar{\sigma}^{2}K.\end{split}$ (44)

###### Proof.

We apply the conclusions from Lemma 6 to Lemma 7. We first show the upper
error bound for
$\sum_{k=0}^{K}\mathbb{E}[\|\bar{\mathbf{u}}_{k}-\overline{\nabla\mathbf{J}}(\mathbf{x}_{k})\|^{2}]$.
Substituting Eq. 29 into Eq. 27 leads to the following inequality

$\begin{split}&\sum_{k=0}^{K}\mathbb{E}[\|\bar{\mathbf{u}}_{k}-\overline{\nabla\mathbf{J}}(\mathbf{x}_{k})\|^{2}]\leq\frac{1}{1-(1-\beta)^{2}}\mathbb{E}[\|\bar{\mathbf{u}}_{0}-\overline{\nabla\mathbf{J}}(\mathbf{x}_{0})\|^{2}]+\frac{12(L^{2}+G^{2}C_{\upsilon})\eta^{2}(1-\beta)^{2}}{N(1-(1-\beta)^{2})}\sum^{K-1}_{k=0}\mathbb{E}[\|\bar{\mathbf{u}}_{k-1}\|^{2}]\\\
&+\frac{2\beta^{2}\bar{\sigma}^{2}K}{N(1-(1-\beta)^{2})}+\frac{12(L^{2}+G^{2}C_{\upsilon})(1-\beta)^{2}}{N^{2}(1-(1-\beta)^{2})}\mathbb{E}\sum^{K-1}_{k=0}[\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}+\|\mathbf{x}_{k-1}-\Lambda\mathbf{x}_{k-1}\|^{2}]\\\
&\leq\frac{1}{1-(1-\beta)^{2}}\mathbb{E}[\|\bar{\mathbf{u}}_{0}-\overline{\nabla\mathbf{J}}(\mathbf{x}_{0})\|^{2}]+\frac{12(L^{2}+G^{2}C_{\upsilon})\eta^{2}(1-\beta)^{2}}{N(1-(1-\beta)^{2})}\sum^{K-1}_{k=0}\mathbb{E}[\|\bar{\mathbf{u}}_{k-1}\|^{2}]\\\
&+\frac{2\beta^{2}\bar{\sigma}^{2}K}{N(1-(1-\beta)^{2})}+\frac{24(L^{2}+G^{2}C_{\upsilon})(1-\beta)^{2}}{N^{2}(1-(1-\beta)^{2})}\sum^{K}_{k=0}\mathbb{E}[\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}].\end{split}$
(45)

We now process the first term on the right hand side of the last inequality.

$\mathbb{E}[\|\bar{\mathbf{u}}_{0}-\overline{\nabla\mathbf{J}}(\mathbf{x}_{0})\|^{2}]=\mathbb{E}[\|\frac{1}{N}\sum_{i=1}^{N}\frac{1}{|\mathcal{B}|}\sum^{|\mathcal{B}|}_{m=1}(\mathbf{g}_{i}(\tau^{i,m}_{0}|\mathbf{x}_{0}^{i})-\nabla
J_{i}(\mathbf{x}_{0}^{i}))\|^{2}].$

It is immediately obtained that

$\mathbb{E}[\|\bar{\mathbf{u}}_{0}-\overline{\nabla\mathbf{J}}(\mathbf{x}_{0})\|^{2}]=\frac{1}{N^{2}|\mathcal{B}|^{2}}\sum^{N}_{i=1}\sum^{|\mathcal{B}|}_{m=1}\mathbb{E}[\|\mathbf{g}_{i}(\tau^{i,m}_{0}|\mathbf{x}_{0}^{i})-\nabla
J_{i}(\mathbf{x}_{0}^{i})\|^{2}]\leq\frac{\bar{\sigma}^{2}}{N|\mathcal{B}|}.$

Observing that $\frac{1}{1-(1-\beta)^{2}}\leq\frac{1}{\beta}$, we can obtain
the first conclusion in Lemma 8. Likewise, if we apply Eq. 27 to Eq. 30, the
following relationship can be attained

$\begin{split}\sum^{K}_{k=1}\mathbb{E}[\|\mathbf{u}_{k}-\nabla\mathbf{J}(\mathbf{x}_{k})\|^{2}]&\leq\frac{\mathbb{E}[\|\mathbf{u}_{0}-\nabla\mathbf{J}(\mathbf{x}_{0})\|^{2}]}{\beta}+\frac{12N(L^{2}+G^{2}C_{\upsilon})\eta^{2}}{\beta}\sum^{K-1}_{k=0}\mathbb{E}[\|\bar{\mathbf{u}}_{k}\|^{2}]\\\
&+\frac{24(L^{2}+G^{2}C_{\upsilon})}{\beta}\sum_{k=0}^{K}\mathbb{E}[\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}]+2N\beta\bar{\sigma}^{2}K.\end{split}$
(46)

As

$\begin{split}\mathbb{E}[\|\mathbf{u}_{0}-\nabla\mathbf{J}(\mathbf{x}_{0})\|^{2}]&=\sum_{i=1}^{N}\mathbb{E}[\|\frac{1}{|\mathcal{B}|}\sum^{|\mathcal{B}|}_{m=1}(\mathbf{g}_{i}(\tau^{i,m}_{0}|\mathbf{x}_{0}^{i})-\nabla
J_{i}(\mathbf{x}_{0}^{i}))\|^{2}]\\\
&=\frac{1}{|\mathcal{B}|^{2}}\sum^{N}_{i=1}\sum^{|\mathcal{B}|}_{m=1}\mathbb{E}[\|\mathbf{g}_{i}(\tau^{i,m}_{0}|\mathbf{x}_{0}^{i})-\nabla
J_{i}(\mathbf{x}_{0}^{i})\|^{2}]\leq\frac{N\bar{\sigma}^{2}}{|\mathcal{B}|}.\end{split}$
(47)

Substituting Eq. 47 into Eq. 46 completes the proof for the second conclusion
in Lemma 8. ∎

It is clearly observed from Eq. 43 that the error bound is dependent of the
consensus error. To bound the consensus, we first present a lemma to bound the
gradient tracking errors. As the initialization approach can be to use either
a single trajectory or a mini-batch of trajectories, the lemma also consists
of the explicit bound for the initial gradient tracking error. Before that, we
present a fact that establishes the relationship between
$\|\mathbf{v}_{k+1}-\Lambda\mathbf{v}_{k+1}\|^{2}$ and
$\|\mathbf{x}_{k+1}-\Lambda\mathbf{x}_{k+1}\|^{2}$.

######  Fact 1.

For all $k\geq 0$, we have

$\begin{split}\|\mathbf{x}_{k+1}-\Lambda\mathbf{x}_{k+1}\|^{2}&\leq\frac{1+\lambda^{2}}{2}\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}\\\
&+\frac{2\eta^{2}\lambda^{2}}{1-\lambda^{2}}\|\mathbf{v}_{k+1}-\Lambda\mathbf{v}_{k+1}\|^{2},\end{split}$
(48)

and

$\begin{split}\|\mathbf{x}_{k+1}-\Lambda\mathbf{x}_{k+1}\|^{2}&\leq
2\lambda^{2}\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}\\\
&+2\eta^{2}\lambda^{2}\|\mathbf{v}_{k+1}-\Lambda\mathbf{v}_{k+1}\|^{2}.\end{split}$
(49)

This fact can be easily shown true by using the update laws of MDPGT and the
Young’s inequality.

######  Lemma 9.

Let $\mathbf{v}_{k}$ be generated by MDPGT initialized with a mini-batch of
trajectories $\mathcal{B}$. We have the following relationships:

$\mathbb{E}[\|\mathbf{v}_{1}-\Lambda\mathbf{v}_{1}\|^{2}]\leq\frac{N\bar{\sigma}^{2}}{|\mathcal{B}|}+\|\nabla\mathbf{J}(\mathbf{x}_{0})\|^{2},$
(50)

and if
$0<\eta\leq\frac{1-\lambda^{2}}{8\lambda\sqrt{9L^{2}+3G^{2}C_{\upsilon}}}$,
then for all $k\leq 1$,

$\begin{split}&\mathbb{E}[\|\mathbf{v}_{k+1}-\Lambda\mathbf{v}_{k+1}\|^{2}]\leq\frac{3+\lambda^{2}}{4}\mathbb{E}[\|\mathbf{v}_{k}-\Lambda\mathbf{v}_{k}\|^{2}]\\\
&+\frac{48L^{2}+24(L^{2}+G^{2}C_{\upsilon})}{1-\lambda^{2}}N\eta^{2}\mathbb{E}[\|\bar{\mathbf{u}}_{k-1}\|^{2}]+4N\beta^{2}\bar{\sigma}^{2}\\\
&+\frac{114L^{2}+72(L^{2}+G^{2}C_{\upsilon})}{1-\lambda^{2}}\mathbb{E}[\|\mathbf{x}_{k-1}-\Lambda\mathbf{x}_{k-1}\|^{2}]\\\
&+\frac{10\beta^{2}}{1-\lambda^{2}}\mathbb{E}[\|\mathbf{u}_{k-1}-\nabla\mathbf{J}(\mathbf{x}_{k-1})\|^{2}].\end{split}$
(51)

###### Proof.

According to the update law, we have

$\mathbf{v}_{1}-\Lambda\mathbf{v}_{1}=\underline{\mathbf{W}}\mathbf{v}_{0}+\mathbf{u}_{0}-\mathbf{u}_{-1}-\Lambda\underline{\mathbf{W}}\mathbf{v}_{0}-\Lambda\mathbf{u}_{0}+\Lambda\mathbf{u}_{-1}=(\mathbf{I}-\Lambda)\mathbf{u}_{0}.$
(52)

It follows from $\mathbf{v}_{0}=\mathbf{0}$ and $\mathbf{u}_{-1}=\mathbf{0}$.
Hence, we have

$\mathbb{E}[\|\mathbf{v}_{1}-\Lambda\mathbf{v}_{1}\|^{2}]=\mathbb{E}[\|(\mathbf{I}-\Lambda)\mathbf{u}_{0}\|^{2}]\leq\mathbb{E}[\|\mathbf{u}_{0}-\nabla\mathbf{J}(\mathbf{x}_{0})+\nabla\mathbf{J}(\mathbf{x}_{0})\|^{2}]$
(53)

Since
$\mathbf{u}^{i}_{0}=\frac{1}{|\mathcal{B}|}\sum^{|\mathcal{B}|}_{m=1}\mathbf{g}_{i}(\tau^{i,m}_{0}|\mathbf{x}^{i}_{0})$,
we then obtain

$\begin{split}&\mathbb{E}[\|\mathbf{u}_{0}-\nabla\mathbf{J}(\mathbf{x}_{0})+\nabla\mathbf{J}(\mathbf{x}_{0})\|^{2}]=\sum^{N}_{i=1}\mathbb{E}[\|\mathbf{u}^{i}_{0}-\nabla
J_{i}(\mathbf{x}^{i}_{0})\|^{2}]+\|\nabla\mathbf{J}(\mathbf{x}_{0})\|^{2}\\\
&=\sum_{i=1}^{N}\mathbb{E}[\|\frac{1}{|\mathcal{B}|}\sum^{|\mathcal{B}|}_{m=1}(\mathbf{g}_{i}(\tau^{i,m}_{0}|\mathbf{x}_{0}^{i})-\nabla
J_{i}(\mathbf{x}_{0}^{i}))\|^{2}]+\|\nabla\mathbf{J}(\mathbf{x}_{0})\|^{2}\\\
&=\frac{1}{|\mathcal{B}|^{2}}\sum^{N}_{i=1}\sum^{|\mathcal{B}|}_{m=1}\mathbb{E}[\|\mathbf{g}_{i}(\tau^{i,m}_{0}|\mathbf{x}_{0}^{i})-\nabla
J_{i}(\mathbf{x}_{0}^{i})\|^{2}]+\|\nabla\mathbf{J}(\mathbf{x}_{0})\|^{2}\\\
&\leq\frac{N\bar{\sigma}^{2}}{|\mathcal{B}|}+\|\nabla\mathbf{J}(\mathbf{x}_{0})\|^{2},\end{split}$
(54)

which completes the first part of the proof.

We now proceed to prove the second part. As

$\begin{split}\mathbf{v}_{k+1}-\Lambda\mathbf{v}_{k+1}&=\underline{\mathbf{W}}\mathbf{v}_{k}+\mathbf{u}_{k}-\mathbf{u}_{k-1}-\Lambda(\underline{\mathbf{W}}\mathbf{v}_{k}+\mathbf{u}_{k}-\mathbf{u}_{k-1})\\\
&=(\underline{\mathbf{W}}-\Lambda)\mathbf{v}_{k}+(\mathbf{I}-\Lambda)(\mathbf{u}_{k}-\mathbf{u}_{k-1}),\end{split}$
(55)

conducting the squared norm results in the following relationship

$\begin{split}\|\mathbf{v}_{k+1}-\Lambda\mathbf{v}_{k+1}\|^{2}&=\|(\underline{\mathbf{W}}-\Lambda)\mathbf{v}_{k}+(\mathbf{I}-\Lambda)(\mathbf{u}_{k}-\mathbf{u}_{k-1})\|^{2}\\\
&=\|(\underline{\mathbf{W}}-\Lambda)\mathbf{v}_{k}\|^{2}+2\langle(\underline{\mathbf{W}}-\Lambda)\mathbf{v}_{k},(\mathbf{I}-\Lambda)(\mathbf{u}_{k}-\mathbf{u}_{k-1})\rangle\\\
&+\|(\mathbf{I}-\Lambda)(\mathbf{u}_{k}-\mathbf{u}_{k-1})\|^{2}\\\
&\leq\lambda^{2}\|\mathbf{v}_{k}-\Lambda\mathbf{v}_{k}\|^{2}+2\langle(\underline{\mathbf{W}}-\Lambda)\mathbf{v}_{k},(\mathbf{I}-\Lambda)(\mathbf{u}_{k}-\mathbf{u}_{k-1})\rangle+\|\mathbf{u}_{k}-\mathbf{u}_{k-1}\|^{2}.\end{split}$
(56)

We next investigate the second term on the right hand side of the last
inequality. Based on the update for the policy gradient surrogate, we have

$\mathbf{u}^{i}_{k}=\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}_{k}^{i})+(1-\beta)\mathbf{u}^{i}_{k-1}-(1-\beta)\upsilon_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{x}^{i}_{k})\mathbf{g}_{i}(\mathbf{x}^{i}_{k-1}).$

Thus,

$\begin{split}\mathbf{u}^{i}_{k}-\mathbf{u}^{i}_{k-1}&=\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}_{k}^{i})-\beta\mathbf{u}^{i}_{k-1}-(1-\beta)\upsilon_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{x}^{i}_{k})\mathbf{g}_{i}(\mathbf{x}^{i}_{k-1})\\\
&=(1-\beta)(\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}_{k}^{i})-\upsilon_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{x}^{i}_{k})\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}_{k-1}^{i}))+\beta(\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}_{k}^{i})-\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}_{k-1}^{i}))\\\
&+\beta(\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}_{k-1}^{i})-\nabla
J(\mathbf{x}^{i}_{k-1})+\nabla
J(\mathbf{x}^{i}_{k-1})-\mathbf{u}^{i}_{k-1}).\end{split}$ (57)

Hence, taking the expectation on both sides of the last equality leads to

$\begin{split}\mathbb{E}[\mathbf{u}_{k}-\mathbf{u}_{k-1}]=(1-\beta)(\nabla\mathbf{J}(\mathbf{x}_{k})-\nabla\mathbf{J}(\mathbf{x}_{k-1}))+\beta\mathbb{E}[\mathbf{g}(\tau_{k}|\mathbf{x}_{k})-\mathbf{g}(\tau_{k}|\mathbf{x}_{k-1})]-\beta(\mathbf{u}_{k-1}-\nabla\mathbf{J}(\mathbf{x}_{k-1})).\end{split}$
(58)

We now investigate the term of
$2\langle(\underline{\mathbf{W}}-\Lambda)\mathbf{v}_{k},(\mathbf{I}-\Lambda)(\mathbf{u}_{k}-\mathbf{u}_{k-1})\rangle$
and denote it as $B_{k}$. Consequently, the following relationship is attained

$\begin{split}\mathbb{E}[B_{k}]&=2\langle(\underline{\mathbf{W}}-\Lambda)\mathbf{v}_{k},(\mathbf{I}-\Lambda)\mathbb{E}[(\mathbf{u}_{k}-\mathbf{u}_{k-1})]\rangle\\\
&=2\langle(\underline{\mathbf{W}}-\Lambda)\mathbf{v}_{k},(\mathbf{I}-\Lambda)[(1-\beta)(\nabla\mathbf{J}(\mathbf{x}_{k})-\nabla\mathbf{J}(\mathbf{x}_{k-1}))+\beta\mathbb{E}[\mathbf{g}(\tau_{k}|\mathbf{x}_{k})-\mathbf{g}(\tau_{k}|\mathbf{x}_{k-1})]\\\
&-\beta(\mathbf{u}_{k-1}-\nabla\mathbf{J}(\mathbf{x}_{k-1}))]\rangle\\\ &\leq
2\lambda\|\mathbf{v}_{k}-\Lambda\mathbf{v}_{k}\|\|(1-\beta)(\nabla\mathbf{J}(\mathbf{x}_{k})-\nabla\mathbf{J}(\mathbf{x}_{k-1}))+\beta\mathbb{E}[\mathbf{g}(\tau_{k}|\mathbf{x}_{k})-\mathbf{g}(\tau_{k}|\mathbf{x}_{k-1})]\\\
&-\beta(\mathbf{u}_{k-1}-\nabla\mathbf{J}(\mathbf{x}_{k-1}))\|\\\
&\leq\frac{1-\lambda^{2}}{2}\|\mathbf{v}_{k}-\Lambda\mathbf{v}_{k}\|^{2}+\frac{6\lambda^{3}(1-\beta)^{2}}{1-\lambda^{2}}\|\nabla\mathbf{J}(\mathbf{x}_{k})-\nabla\mathbf{J}(\mathbf{x}_{k-1})\|^{2}\\\
&+\frac{6\lambda^{3}\beta^{2}}{1-\lambda^{2}}\mathbb{E}[\|\mathbf{g}(\tau_{k}|\mathbf{x}_{k})-\mathbf{g}(\tau_{k}|\mathbf{x}_{k-1})\|^{2}]+\frac{6\lambda^{3}\beta^{2}}{1-\lambda^{2}}\|\mathbf{u}_{k-1}-\nabla\mathbf{J}(\mathbf{x}_{k-1})\|^{2}\\\
&\leq\frac{1-\lambda^{2}}{2}\|\mathbf{v}_{k}-\Lambda\mathbf{v}_{k}\|^{2}+(\beta^{2}+(1-\beta)^{2})\frac{6\lambda^{3}L^{2}}{1-\lambda^{2}}\mathbb{E}[\|\mathbf{x}_{k}-\mathbf{x}_{k-1}\|^{2}]+\frac{6\lambda^{3}\beta^{2}}{1-\lambda^{2}}\|\mathbf{u}_{k-1}-\nabla\mathbf{J}(\mathbf{x}_{k-1})\|^{2}.\end{split}$
(59)

The first inequality is due to Cauthy-Schwartz inequality, while the second
inequality follows from basic inequalities $2ab\leq ea^{2}+b^{2}/e$, with
$e=\frac{1-\lambda^{2}}{2\lambda}$ for all $a,b\in\mathbb{R}$ and
$\|\mathbf{a}+\mathbf{b}+\mathbf{c}\|^{2}\leq
3\|\mathbf{a}\|^{2}+3\|\mathbf{b}\|^{2}+3\|\mathbf{c}\|^{2},\forall\mathbf{a},\mathbf{b},\mathbf{c}\in\mathbb{R}^{d}$.
The last inequality is due to the smoothness property. Using Eq. 57 yields the
upper bound of the second moment of $\mathbf{u}^{i}_{k}-\mathbf{u}^{i}_{k-1}$,

$\begin{split}&\mathbb{E}[\|\mathbf{u}^{i}_{k}-\mathbf{u}^{i}_{k-1}\|^{2}]\leq
4(1-\beta)^{2}\mathbb{E}[\|\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}_{k}^{i})-\upsilon_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{x}^{i}_{k})\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}_{k-1}^{i})\|^{2}]\\\
&+4\beta^{2}\mathbb{E}[\|\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}_{k}^{i})-\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}_{k-1}^{i})\|^{2}]+4\beta^{2}\mathbb{E}[\|\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}_{k-1}^{i})-\nabla
J(\mathbf{x}^{i}_{k-1})\|^{2}]\\\ &+4\beta^{2}\mathbb{E}[\|\nabla
J(\mathbf{x}^{i}_{k-1})-\mathbf{u}^{i}_{k-1}\|^{2}]\\\
&=4(1-\beta)^{2}\mathbb{E}[\|\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}_{k}^{i})-\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}_{k-1}^{i})+\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}_{k-1}^{i})-\upsilon_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{x}^{i}_{k})\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}_{k-1}^{i})\|^{2}]\\\
&+4\beta^{2}\mathbb{E}[\|\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}_{k}^{i})-\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}_{k-1}^{i})\|^{2}]+4\beta^{2}\mathbb{E}[\|\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}_{k-1}^{i})-\nabla
J(\mathbf{x}^{i}_{k-1})\|^{2}]\\\ &+4\beta^{2}\mathbb{E}[\|\nabla
J(\mathbf{x}^{i}_{k-1})-\mathbf{u}^{i}_{k-1}\|^{2}]\\\ &\leq
4(1-\beta)^{2}\mathbb{E}[2\|\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}_{k}^{i})-\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}_{k-1}^{i})\|^{2}+2\|\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}_{k-1}^{i})-\upsilon_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{x}^{i}_{k})\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}_{k-1}^{i})\|^{2}]\\\
&+4\beta^{2}\mathbb{E}[\|\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}_{k}^{i})-\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}_{k-1}^{i})\|^{2}]+4\beta^{2}\mathbb{E}[\|\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}_{k-1}^{i})-\nabla
J(\mathbf{x}^{i}_{k-1})\|^{2}]\\\ &+4\beta^{2}\mathbb{E}[\|\nabla
J(\mathbf{x}^{i}_{k-1})-\mathbf{u}^{i}_{k-1}\|^{2}]\\\
&=8(1-\beta)^{2}\mathbb{E}[\|\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}_{k}^{i})-\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}_{k-1}^{i})\|^{2}]\\\
&+8(1-\beta)^{2}\mathbb{E}[\|\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}_{k-1}^{i})-\upsilon_{i}(\tau^{i}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{x}^{i}_{k})\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}_{k-1}^{i})\|^{2}]\\\
&+4\beta^{2}\mathbb{E}[\|\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}_{k}^{i})-\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}_{k-1}^{i})\|^{2}]+4\beta^{2}\mathbb{E}[\|\mathbf{g}_{i}(\tau^{i}_{k}|\mathbf{x}_{k-1}^{i})-\nabla
J(\mathbf{x}^{i}_{k-1})\|^{2}]\\\ &+4\beta^{2}\mathbb{E}[\|\nabla
J(\mathbf{x}^{i}_{k-1})-\mathbf{u}^{i}_{k-1}\|^{2}]\\\ &\leq
8(1-\beta)^{2}L^{2}\mathbb{E}[\|\mathbf{x}^{i}_{k}-\mathbf{x}^{i}_{k-1}\|^{2}]+8(1-\beta)^{2}G^{2}C_{\upsilon}\mathbb{E}[\|\mathbf{x}^{i}_{k}-\mathbf{x}^{i}_{k-1}\|^{2}]\\\
&+4\beta^{2}L^{2}\mathbb{E}[\|\mathbf{x}^{i}_{k}-\mathbf{x}^{i}_{k-1}\|^{2}]+4\beta^{2}\sigma^{2}_{i}+4\beta^{2}\mathbb{E}[\|\nabla
J(\mathbf{x}^{i}_{k-1})-\mathbf{u}^{i}_{k-1}\|^{2}]\\\
&=[8(1-\beta)^{2}L^{2}+8(1-\beta)^{2}G^{2}C_{\upsilon}+4\beta^{2}L^{2}]\mathbb{E}[\|\mathbf{x}^{i}_{k}-\mathbf{x}^{i}_{k-1}\|]\\\
&+4\beta^{2}\sigma^{2}_{i}+4\beta^{2}\mathbb{E}[\|\nabla
J(\mathbf{x}^{i}_{k-1})-\mathbf{u}^{i}_{k-1}\|^{2}].\end{split}$ (60)

The last inequality is based on the smoothness property, Lemma 1 and Lemma 2.
We then have that

$\begin{split}\mathbb{E}[\|\mathbf{u}_{k}-\mathbf{u}_{k-1}\|^{2}]&\leq[8(1-\beta)^{2}L^{2}+8(1-\beta)^{2}G^{2}C_{\upsilon}+4\beta^{2}L^{2}]\mathbb{E}[\|\mathbf{x}_{k}-\mathbf{x}_{k-1}\|^{2}]+\\\
&4\beta^{2}N\bar{\sigma}^{2}+4\beta^{2}\mathbb{E}[\|\nabla\mathbf{J}(\mathbf{x}_{k-1}-\mathbf{u}_{k-1})\|^{2}].\end{split}$
(61)

Taking the expectation on both sides of Eq. 56 and substituting Eq. 59 and Eq.
61 into it, the following relationship can be obtained

$\begin{split}&\mathbb{E}[\|\mathbf{v}_{k+1}-\Lambda\mathbf{v}_{k+1}\|^{2}]\leq\lambda^{2}\mathbb{E}[\|\mathbf{v}_{k}-\Lambda\mathbf{v}_{k}\|^{2}]+\frac{1-\lambda^{2}}{2}\mathbb{E}[\|\mathbf{v}_{k}-\Lambda\mathbf{v}_{k}\|^{2}]\\\
&+(\beta^{2}+(1-\beta)^{2})\frac{6\lambda^{3}L^{2}}{1-\lambda^{2}}\mathbb{E}[\|\mathbf{x}_{k}-\mathbf{x}_{k-1}\|^{2}]\\\
&+\frac{6\lambda^{3}\beta^{2}}{1-\lambda^{2}}\mathbb{E}[\|\mathbf{u}_{k-1}-\nabla\mathbf{J}(\mathbf{x}_{k-1})\|^{2}]+4\beta^{2}N\bar{\sigma}^{2}\\\
&+[8(1-\beta)^{2}L^{2}+8(1-\beta)^{2}G^{2}C_{\upsilon}+4\beta^{2}L^{2}]\mathbb{E}[\|\mathbf{x}_{k}-\mathbf{x}_{k-1}\|^{2}]\\\
&+4\beta^{2}\mathbb{E}[\|\mathbf{u}_{k-1}-\nabla\mathbf{J}(\mathbf{x}_{k-1})\|^{2}]\\\
&=\frac{1+\lambda^{2}}{2}\lambda^{2}\mathbb{E}[\|\mathbf{v}_{k}-\Lambda\mathbf{v}_{k}\|^{2}]\\\
&+\bigg{[}(\beta^{2}+(1-\beta)^{2})\frac{6\lambda^{3}L^{2}}{1-\lambda^{2}}+8(1-\beta)^{2}L^{2}+8(1-\beta)^{2}G^{2}C_{\upsilon}+4\beta^{2}L^{2}\bigg{]}\mathbb{E}[\|\mathbf{x}_{k}-\mathbf{x}_{k-1}\|^{2}]\\\
&+\bigg{(}4\beta^{2}+\frac{6\lambda^{3}\beta^{2}}{1-\lambda^{2}}\bigg{)}\mathbb{E}[\|\mathbf{u}_{k-1}-\nabla\mathbf{J}(\mathbf{x}_{k-1})\|^{2}]+4\beta^{2}N\bar{\sigma}^{2}\\\
&\leq\frac{1+\lambda^{2}}{2}\lambda^{2}\mathbb{E}[\|\mathbf{v}_{k}-\Lambda\mathbf{v}_{k}\|^{2}]+\bigg{(}\frac{12\lambda^{3}L^{2}}{1-\lambda^{2}}+8(L^{2}+G^{2}C_{\upsilon})+4L^{2}\bigg{)}\mathbb{E}[\|\mathbf{x}_{k}-\mathbf{x}_{k-1}\|^{2}]\\\
&+4N\beta^{2}\bar{\sigma}^{2}+\frac{10\beta^{2}}{1-\lambda^{2}}\mathbb{E}[\|\mathbf{u}_{k-1}-\nabla\mathbf{J}(\mathbf{x}_{k-1})\|^{2}]\\\
&\leq\frac{1+\lambda^{2}}{2}\lambda^{2}\mathbb{E}[\|\mathbf{v}_{k}-\Lambda\mathbf{v}_{k}\|^{2}]+\frac{16L^{2}+8(L^{2}+G^{2}C_{\upsilon})}{1-\lambda^{2}}\mathbb{E}[\|\mathbf{x}_{k}-\mathbf{x}_{k-1}\|^{2}]+4\beta^{2}N\bar{\sigma}^{2}\\\
&+\frac{10\beta^{2}}{1-\lambda^{2}}\mathbb{E}[\|\mathbf{u}_{k-1}-\nabla\mathbf{J}(\mathbf{x}_{k-1})\|^{2}].\end{split}$
(62)

Since

$\|\mathbf{x}_{k}-\mathbf{x}_{k-1}\|^{2}\leq
3\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}+3N\eta^{2}\|\bar{\mathbf{u}}_{k-1}\|^{2}+3\|\mathbf{x}_{k-1}-\Lambda\mathbf{x}_{k-1}\|^{2},$

combining Fact 1, the following is obtained

$\|\mathbf{x}_{k}-\mathbf{x}_{k-1}\|^{2}\leq
6\lambda^{2}\eta^{4}\|\mathbf{v}_{k}-\Lambda\mathbf{v}_{k}\|^{2}+3N\eta^{2}\|\bar{\mathbf{u}}_{k-1}\|^{2}+9\|\mathbf{x}_{k-1}-\Lambda\mathbf{x}_{k-1}\|^{2}.$
(63)

Combining the last two inequalities, one can attain

$\begin{split}&\mathbb{E}[\|\mathbf{v}_{k+1}-\Lambda\mathbf{v}_{k+1}\|^{2}]\leq\bigg{(}\frac{1+\lambda^{2}}{2}+\frac{96L^{2}+48(L^{2}+G^{2}C_{\upsilon})}{1-\lambda^{2}}\lambda^{2}\eta^{2}\bigg{)}\mathbb{E}[\|\mathbf{v}_{k}-\Lambda\mathbf{v}_{k}\|^{2}]\\\
&+\frac{48L^{2}+24(L^{2}+G^{2}C_{\upsilon})}{1-\lambda^{2}}N\eta^{2}\mathbb{E}[\|\bar{\mathbf{u}}_{k-1}\|^{2}]+\frac{144L^{2}+72(L^{2}+G^{2}C_{\upsilon})}{1-\lambda^{2}}\mathbb{E}[\|\mathbf{x}_{k-1}-\Lambda\mathbf{x}_{k-1}\|^{2}]\\\
&+4\beta^{2}N\bar{\sigma}^{2}+\frac{10\beta^{2}}{1-\lambda^{2}}\mathbb{E}[\|\mathbf{u}_{k-1}-\nabla\mathbf{J}(\mathbf{x}_{k-1})\|^{2}].\end{split}$
(64)

Combining the fact that
$0<\eta\leq\frac{1-\lambda^{2}}{8\lambda\sqrt{9L^{2}+3G^{2}C_{\upsilon}}}$
yields the desirable result. ∎

According to Lemma 8 and Lemma 9, the remaining step to obtain the explicitly
accurate error bound is to inaugurate the correlation between
$\sum_{k=0}^{K}\mathbb{E}[\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}]$ and
$\sum_{k=0}^{K}\mathbb{E}[\|\bar{\mathbf{u}}_{k}\|^{2}]$. Thus, the following
lemma is constructed for this purpose.

######  Lemma 10.

Let $\mathbf{x}_{k}$ be generated by MDPGT initialized with a mini-batch of
trajectories $\mathcal{B}$. If
$0<\eta\leq\frac{(1-\lambda^{2})^{2}}{\lambda\sqrt{14592L^{2}+9984G^{2}C_{\upsilon}}}$
and $\beta\in(0,1)$, then for $K\geq 2$, we have

$\begin{split}&\sum_{k=0}^{K}\mathbb{E}[\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}]\leq\Bigg{[}\frac{1436L^{2}+4608(L^{2}+G^{2}C_{\upsilon})}{(1-\lambda^{2})^{4}}\Bigg{]}\cdot\\\
&\lambda^{2}N\eta^{4}\sum_{k=0}^{K-2}\mathbb{E}[\|\bar{\mathbf{u}}_{k}\|^{2}]+\frac{32\lambda^{2}N\bar{\sigma}^{2}\eta^{2}}{(1-\lambda^{2})^{3}|\mathcal{B}|}\bigg{(}1+\frac{10\beta}{1-\lambda^{2}}\bigg{)}\\\
&+\frac{32\lambda^{2}\eta^{2}}{(1-\lambda^{2})^{3}}\|\nabla\mathbf{J}(\mathbf{x}_{0})\|^{2}+\frac{128\lambda^{2}N\beta^{2}K\bar{\sigma}^{2}\eta^{2}}{(1-\lambda^{2})^{3}}\bigg{(}1+\frac{5\beta}{1-\lambda^{2}}\bigg{)}.\end{split}$
(65)

###### Proof.

Applying Eq. 27 to the fact that

$\|\mathbf{x}_{k+1}-\Lambda\mathbf{x}_{k+1}\|^{2}\leq\frac{1+\lambda^{2}}{2}\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}+\frac{2\eta^{2}\lambda^{2}}{1-\lambda^{2}}\|\mathbf{v}_{k+1}-\Lambda\mathbf{v}_{k+1}\|^{2}$

leads to

$\sum_{k=0}^{K}\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}\leq\frac{4\lambda^{2}\eta^{2}}{(1-\lambda^{2})^{2}}\sum^{K}_{k=1}\|\mathbf{v}_{k}-\Lambda\mathbf{v}_{k}\|^{2}.$

Similarly, applying Eq. 28 to Eq. 51 results in

$\begin{split}&\sum_{k=1}^{K}\mathbb{E}[\|\mathbf{v}_{k}-\Lambda\mathbf{v}_{k}\|^{2}]\leq\frac{4}{1-\lambda^{2}}\mathbb{E}[\|\mathbf{v}_{1}-\Lambda\mathbf{v}_{1}\|^{2}]+\frac{192L^{2}+96(L^{2}+G^{2}C_{\upsilon})}{1-\lambda^{2}}N\eta^{2}\sum_{k=0}^{K-2}\mathbb{E}[\|\bar{\mathbf{u}}_{k-1}\|^{2}]\\\
&+\frac{576L^{2}+288(L^{2}+G^{2}C_{\upsilon})}{(1-\lambda^{2})^{2}}\sum_{k=0}^{K-2}\mathbb{E}[\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}]+\frac{40\beta^{2}}{(1-\lambda^{2})^{2}}\sum_{k=0}^{K-2}\mathbb{E}[\|\mathbf{u}_{k}-\nabla\mathbf{J}(\mathbf{x}_{k})\|^{2}]\\\
&+\frac{16\beta^{2}N\bar{\sigma}^{2}K}{1-\lambda^{2}}\\\
&\leq\frac{4N\bar{\sigma}^{2}}{(1-\lambda^{2})|\mathcal{B}|}+\frac{4\|\nabla\mathbf{J}(\mathbf{x}_{0})\|^{2}}{1-\lambda^{2}}+\frac{192L^{2}+96(L^{2}+G^{2}C_{\upsilon})}{1-\lambda^{2}}N\eta^{2}\sum_{k=0}^{K-2}\mathbb{E}[\|\bar{\mathbf{u}}_{k-1}\|^{2}]\\\
&+\frac{576L^{2}+288(L^{2}+G^{2}C_{\upsilon})}{(1-\lambda^{2})^{2}}\sum_{k=0}^{K-2}\mathbb{E}[\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}]+\frac{40\beta^{2}}{(1-\lambda^{2})^{2}}\sum_{k=0}^{K-2}\mathbb{E}[\|\mathbf{u}_{k}-\nabla\mathbf{J}(\mathbf{x}_{k})\|^{2}]\\\
&+\frac{16\beta^{2}N\bar{\sigma}^{2}K}{1-\lambda^{2}}.\end{split}$ (66)

The last inequality follows from Lemma 9. As

$\begin{split}\frac{40\beta^{2}}{(1-\lambda^{2})^{2}}&\sum_{k=0}^{K-2}\mathbb{E}[\|\mathbf{u}_{k}-\nabla\mathbf{J}(\mathbf{x}_{k})\|^{2}]\leq\frac{40N\bar{\sigma}^{2}\beta}{(1-\lambda^{2})^{2}|\mathcal{B}|}+\frac{480\beta
N(L^{2}+G^{2}C_{\upsilon})\eta^{2}}{(1-\lambda^{2})^{2}}\\\
&\sum_{k=0}^{K-1}\mathbb{E}[\|\bar{\mathbf{u}}_{k}\|^{2}]+\frac{960\beta(L^{2}+G^{2}C_{\upsilon})}{(1-\lambda^{2})^{2}}\sum_{k=0}^{K}\mathbb{E}[\|\mathbf{x}_{k}-\Lambda(\mathbf{x}_{k})\|^{2}]\\\
&+\frac{80N\beta^{3}\bar{\sigma}^{2}K}{(1-\lambda^{2})^{2}},\end{split}$ (67)

we then have

$\begin{split}\sum_{k=1}^{K}\mathbb{E}[\|\mathbf{v}_{k}-\Lambda\mathbf{v}_{k}\|^{2}]&\leq\frac{192L^{2}+576(L^{2}+G^{2}C_{\upsilon})}{(1-\lambda^{2})^{2}}N\eta^{2}\sum_{k=0}^{K-2}\mathbb{E}[\|\bar{\mathbf{u}}_{k}\|^{2}]\\\
&+\frac{576L^{2}+1248(L^{2}+G^{2}C_{\upsilon})}{(1-\lambda^{2})^{2}}\sum^{k=0}_{K-1}\mathbb{E}[\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}]\\\
&\frac{4N\bar{\sigma}^{2}}{(1-\lambda^{2})|\mathcal{B}|}\bigg{(}1+\frac{10\beta}{1-\lambda^{2}}\bigg{)}+\frac{16\beta^{2}N\bar{\sigma}^{2}K}{1-\lambda^{2}}\bigg{(}1+\frac{5\beta}{1-\lambda^{2}}\bigg{)}+\frac{4\|\nabla\mathbf{J}(\mathbf{x}_{0})\|^{2}}{1-\lambda^{2}}.\end{split}$
(68)

Thus, we have

$\begin{split}&\sum_{k=0}^{K}\mathbb{E}[\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}]\leq\frac{768L^{2}+2304(L^{2}+G^{2}C_{\upsilon})\lambda^{2}N\eta^{4}}{(1-\lambda^{2})^{4}}\sum^{K-2}_{k=0}\mathbb{E}[\|\bar{\mathbf{u}}_{k}\|^{2}]\\\
&+\frac{2304L^{2}+4992(L^{2}+G^{2}C_{\upsilon})}{(1-\lambda^{2})^{4}}\lambda^{2}\eta^{2}\sum^{K-1}_{k=0}\mathbb{E}[\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}]+\frac{16N\bar{\sigma}^{2}\lambda^{2}\eta^{2}}{(1-\lambda^{2})^{3}|\mathcal{B}|}\bigg{(}1+\frac{10\beta}{1-\lambda^{2}}\bigg{)}\\\
&+\frac{64\beta^{2}N\bar{\sigma}^{2}K\eta^{2}\lambda^{2}}{(1-\lambda^{2})^{3}}\bigg{(}1+\frac{5\beta}{1-\lambda^{2}}\bigg{)}+\frac{16\lambda^{2}\eta^{2}\|\nabla\mathbf{J}(\mathbf{x}_{0})\|^{2}}{(1-\lambda^{2})^{3}}.\end{split}$
(69)

With simple mathematical manipulations, based on the conditions
$0<\eta\leq\frac{(1-\lambda^{2})^{2}}{\lambda\sqrt{14592L^{2}+9984G^{2}C_{\upsilon}}}$
and $\beta\in(0,1)$, it is obtained that

$1-\frac{2304L^{2}+4992(L^{2}+G^{2}C_{\upsilon})}{(1-\lambda^{2})^{4}}\lambda^{2}\eta^{2}\leq\frac{1}{2}.$

Hence the proof is completed by adopting this inequality. ∎

Hence, with Lemma 10, it suffices to show Theorem 1, whose proof is presented
next. While for Theorem 2, the same proof ideas can be applied and the
corresponding results are obtained by setting $|\mathcal{B}|=1$. Thus, we are
not going to repeat statements for the auxiliary lemmas, instead using the
conclusions from the lemmas and adjusting slightly the constants in the error
bounds.

### Proof of Theorem 1

###### Proof.

Recall the conclusion of Lemma 5, we have

$\sum_{k=0}^{K}\|\nabla
J(\bar{\mathbf{x}}_{k})\|^{2}\leq\frac{2(J^{*}-J(\bar{\mathbf{x}}_{0}))}{\eta}-\frac{1}{2}\sum_{k=0}^{K}\|\bar{\mathbf{u}}_{k}\|^{2}+2\sum_{k=0}^{K}\|\bar{\mathbf{u}}_{k}-\overline{\nabla\mathbf{J}}(\mathbf{x}_{k})\|^{2}+\frac{2L^{2}}{N}\sum_{k=0}^{K}\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}.$

Substituting Eq. 43 into the last inequality yields

$\begin{split}&\sum_{k=0}^{K}\|\nabla
J(\bar{\mathbf{x}}_{k})\|^{2}\leq\frac{2(J^{*}-J(\bar{\mathbf{x}}_{0}))}{\eta}-\frac{1}{2}\sum_{k=0}^{K}\mathbb{E}[\|\bar{\mathbf{u}}_{k}\|^{2}]+\frac{2L^{2}}{N}\sum^{K}_{k=0}\mathbb{E}[\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}]\\\
&+\frac{2\bar{\sigma}^{2}}{N\beta|\mathcal{B}|}+\frac{24(L^{2}+G^{2}C_{\upsilon})\eta^{2}}{N\beta}\sum_{k=0}^{K-1}\mathbb{E}[\|\bar{\mathbf{u}}_{k}\|^{2}]+\frac{48(L^{2}+G^{2}C_{\upsilon})}{N^{2}\beta}\sum^{K}_{k=0}\mathbb{E}[\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}]+\frac{4\beta\bar{\sigma}^{2}K}{N}\\\
&=\frac{2(J^{*}-J(\bar{\mathbf{x}}_{0}))}{\eta}-\frac{1}{4}\sum_{k=0}^{K}\mathbb{E}[\|\bar{\mathbf{u}}_{k}\|^{2}]+\frac{2}{N}\bigg{(}L^{2}+\frac{24(L^{2}+G^{2}C_{\upsilon})}{N\beta}\bigg{)}\sum^{K}_{k=0}\mathbb{E}[\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}]\\\
&+\frac{2\bar{\sigma}^{2}}{N\beta|\mathcal{B}|}+\frac{4\beta\bar{\sigma}^{2}K}{N}-\bigg{(}\frac{1}{4}-\frac{24(L^{2}+G^{2}C_{\upsilon})\eta^{2}}{N\beta}\bigg{)}\sum_{k=0}^{K}\mathbb{E}[\|\bar{\mathbf{u}}_{k}\|^{2}].\end{split}$
(70)

To get rid of the last term on the right hand side of the above inequality,
the following relationship can be obtained

$\frac{1}{4}-\frac{24(L^{2}+G^{2}C_{\upsilon})\eta^{2}}{N\beta}\geq
0\Rightarrow\frac{96(L^{2}+G^{2}C_{\upsilon})\eta^{2}}{N}\leq\beta<1\Rightarrow
0<\eta<\frac{1}{6\sqrt{6(L^{2}+G^{2}C_{\upsilon})}}.$

It is easily to verified the above relationship as we have that
$\beta=\frac{96(L^{2}+G^{2}C_{\upsilon})\eta^{2}}{N}$. Thus,

$\begin{split}&\sum_{k=0}^{K}\|\nabla
J(\bar{\mathbf{x}}_{k})\|^{2}\leq\frac{2(J^{*}-J(\bar{\mathbf{x}}_{0}))}{\eta}-\frac{1}{4}\sum_{k=0}^{K}\mathbb{E}[\|\bar{\mathbf{u}}_{k}\|^{2}]+\frac{2\bar{\sigma}^{2}}{N\beta|\mathcal{B}|}+\frac{4\beta\bar{\sigma}^{2}K}{N}\\\
&+\frac{2}{N}\bigg{(}L^{2}+\frac{24(L^{2}+G^{2}C_{\upsilon})}{N\beta}\bigg{)}\sum^{K}_{k=0}\mathbb{E}[\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}].\end{split}$
(71)

As

$\begin{split}&\frac{1}{N}\sum_{i=1}^{N}\sum_{k=0}^{K}\mathbb{E}[\|\nabla
J(\mathbf{x}^{i}_{k})\|^{2}]\leq\frac{2}{N}\sum^{N}_{i=1}\sum^{K}_{k=0}\mathbb{E}[\|\nabla
J(\mathbf{x}^{i}_{k})-\nabla J(\bar{\mathbf{x}}_{k})\|^{2}+\|\nabla
J(\bar{\mathbf{x}}_{k})\|^{2}]\\\
&\leq\frac{2L^{2}}{N}\sum^{K}_{k=0}\mathbb{E}[\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}]+2\sum^{K}_{k=0}\mathbb{E}[\|\nabla
J(\bar{\mathbf{x}}_{k})\|^{2}],\end{split}$ (72)

we have

$\begin{split}&\frac{1}{N}\sum_{i=1}^{N}\sum_{k=0}^{K}\mathbb{E}[\|\nabla
J(\mathbf{x}^{i}_{k})\|^{2}]\leq\frac{2L^{2}}{N}\sum^{K}_{k=0}\mathbb{E}[\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}]+\frac{4(J^{*}-J(\bar{\mathbf{x}}_{0}))}{\eta}\\\
&-\frac{1}{2}\sum_{k=0}^{K}\mathbb{E}[\|\bar{\mathbf{u}}_{k}\|^{2}]+\frac{4\bar{\sigma}^{2}}{N\beta|\mathcal{B}|}+\frac{8\beta\bar{\sigma}^{2}K}{N}+\frac{4}{N}\bigg{(}L^{2}+\frac{24(L^{2}+G^{2}C_{\upsilon})}{N\beta}\bigg{)}\sum^{K}_{k=0}\mathbb{E}[\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}]\\\
&=\frac{6}{N}\bigg{(}L^{2}+\frac{16(L^{2}+G^{2}C_{\upsilon})}{N\beta}\bigg{)}\sum^{K}_{k=0}\mathbb{E}[\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}]+\frac{4(J^{*}-J(\bar{\mathbf{x}}_{0}))}{\eta}-\frac{1}{2}\sum_{k=0}^{K}\mathbb{E}[\|\bar{\mathbf{u}}_{k}\|^{2}]\\\
&+\frac{4\bar{\sigma}^{2}}{N\beta|\mathcal{B}|}+\frac{8\beta\bar{\sigma}^{2}K}{N}.\end{split}$
(73)

As
$\frac{6}{N}\bigg{(}L^{2}+\frac{16(L^{2}+G^{2}C_{\upsilon})}{N\beta}\bigg{)}=\frac{6}{N}\bigg{(}L^{2}+\frac{1}{6\eta^{2}L^{2}}\bigg{)}$
and $\eta^{2}L^{2}\leq\eta^{2}(L^{2}+G^{2}C_{\upsilon})\leq\frac{1}{96}$ based
on $\eta<\frac{1}{6\sqrt{6(L^{2}+G^{2}C_{\upsilon})}}$, we have

$\frac{6}{N}\bigg{(}L^{2}+\frac{16(L^{2}+G^{2}C_{\upsilon})}{N\beta}\bigg{)}<\frac{17}{16N\eta^{2}}.$

Thus, the following relationship can be attained

$\begin{split}\frac{1}{N}\sum_{i=1}^{N}\sum_{k=0}^{K}\mathbb{E}[\|\nabla
J(\mathbf{x}^{i}_{k})\|^{2}]&\leq\frac{17}{16\eta^{2}N}\mathbb{E}[\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}]+\frac{4(J^{*}-J(\bar{\mathbf{x}}_{0}))}{\eta}-\frac{1}{2}\sum_{k=0}^{K}\mathbb{E}[\|\bar{\mathbf{u}}_{k}\|^{2}]\\\
&+\frac{4\bar{\sigma}^{2}}{N\beta|\mathcal{B}|}+\frac{8\beta\bar{\sigma}^{2}K}{N}.\end{split}$
(74)

We investigate the combined term
$-\frac{1}{2}\sum_{k=0}^{K}\mathbb{E}[\|\bar{\mathbf{u}}_{k}\|^{2}]+\frac{17}{16\eta^{2}N}\mathbb{E}[\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}]$.
Since

$\begin{split}&-\frac{1}{2}\sum_{k=0}^{K}\mathbb{E}[\|\bar{\mathbf{u}}_{k}\|^{2}]+\frac{17}{16\eta^{2}N}\mathbb{E}[\|\mathbf{x}_{k}-\Lambda\mathbf{x}_{k}\|^{2}]\leq-\frac{1}{2}\sum_{k=0}^{K}\mathbb{E}[\|\bar{\mathbf{u}}_{k}\|^{2}]\\\
&+\frac{17}{16\eta^{2}N}\bigg{[}\bigg{(}\frac{1436L^{2}+4608(L^{2}+G^{2}C_{\upsilon})}{(1-\lambda^{2})^{4}}\bigg{)}\lambda^{2}N\eta^{4}\sum_{k=0}^{K-2}\mathbb{E}[\|\bar{\mathbf{u}}_{k}\|^{2}]+\frac{32\lambda^{2}N\bar{\sigma}^{2}\eta^{2}}{(1-\lambda^{2})^{3}|\mathcal{B}|}\\\
&\bigg{(}1+\frac{10\beta}{1-\lambda^{2}}\bigg{)}+\frac{32\lambda^{2}\eta^{2}}{(1-\lambda^{2})^{3}}\|\nabla\mathbf{J}(\mathbf{x}_{0})\|^{2}+\frac{128\lambda^{2}N\beta^{2}K\bar{\sigma}^{2}\eta^{2}}{(1-\lambda^{2})^{3}}\bigg{(}1+\frac{5\beta}{1-\lambda^{2}}\bigg{)}\bigg{]}\\\
&=-\frac{1}{2}\sum_{k=0}^{K}\mathbb{E}[\|\bar{\mathbf{u}}_{k}\|^{2}]+\frac{1526L^{2}+4896(L^{2}+G^{2}C_{\upsilon})}{(1-\lambda^{2})^{4}}\lambda^{2}\eta^{2}\sum_{k=0}^{K-2}\mathbb{E}[\|\bar{\mathbf{u}}_{k}\|^{2}]\\\
&+\frac{34\lambda^{2}}{N(1-\lambda^{2})^{3}}\|\nabla\mathbf{J}(\mathbf{x}_{0})\|^{2}+\frac{34\lambda^{2}\bar{\sigma}^{2}}{(1-\lambda^{2})^{3}|\mathcal{B}|}\bigg{(}1+\frac{10\beta}{1-\lambda^{2}}\bigg{)}+\frac{136\lambda^{2}\beta^{2}K\bar{\sigma}^{2}}{(1-\lambda^{2})^{3}}\bigg{(}1+\frac{5\beta}{1-\lambda^{2}}\bigg{)}\\\