# A divide-and-conquer approach for
sparse recovery of high dimensional signals

Aron Bevelander, Kim Batselier, Nitin Jonathan Myers
Delft Center for Systems and Control, Delft University of Technology, The
Netherlands
Email<EMAIL_ADDRESS>{K.Batselier<EMAIL_ADDRESS>
###### Abstract

Compressed sensing (CS) techniques demand significant storage and
computational resources, when recovering high-dimensional sparse signals.
Block CS (BCS), a special class of CS, addresses both the storage and
complexity issues by partitioning the sparse recovery problem into several
sub-problems. In this paper, we derive a Welch bound-based guarantee on the
reconstruction error with BCS. Our guarantee reveals that the reconstruction
quality with BCS monotonically reduces with an increasing number of
partitions. To alleviate this performance loss, we propose a sparse recovery
technique that exploits correlation across the partitions of the sparse
signal. Our method outperforms BCS in the moderate SNR regime, for a modest
increase in the storage and computational complexities.

###### Index Terms:

Block compressed sensing, computational complexity, memory limitation

##  I Introduction

Compressed sensing (CS) is a technique to recover sparse high-dimensional
signals from their low-dimensional representation. The dimension of sparse
signals in typical CS applications is rapidly increasing with advances in
sensing technology. For instance, sparse signals in a million dimension space
have been considered in the context of fluorescence microscopy [1], imaging
[2], and wireless channel estimation using terahertz arrays [3]. Solving for
such high-dimensional signals in real-time requires intense memory and
computational resources that may often be impractical or unavailable.

Block compressed sensing (BCS), a class of CS, adopts a divide-and-conquer
approach to reduce the storage and computational complexities needed for
sparse recovery [4]. In BCS, the sparse recovery problem is partitioned into
sub-problems, where each sub-problem solves for a block within the sparse
vector. The sub-problems have a lower complexity than the original problem and
they can be solved in parallel. BCS was applied in [5] to reduce the storage
needed for sparse recovery. Standard BCS methods typically obtain an equal
number of CS measurements for each block, assuming a uniform distribution of
sparse non-zero components across the blocks. In [6], however, varying numbers
of measurements were acquired per block based on its expected sparsity level.
In [5], a permutation technique was proposed to “equally” distribute the non-
zero components over different blocks. An important question in BCS is if
partitioning the problem impacts reconstruction performance. An empirical
study conducted in [7] showed that partitioning results in a higher mean
squared error in the reconstruction. To the best of our knowledge, however,
the fundamental limits on reconstruction performance with BCS as a function of
the number of partitions or blocks have not been studied.

To overcome the performance loss with standard BCS, the sub-problems of BCS
can be solved in series. Assuming that the signal is correlated across the
blocks, information from one reconstructed block can be used as a prior in
reconstructing other correlated blocks. Prior work on dynamical CS has
extensively looked at exploiting such correlation information in a different
context than BCS. For instance, the techniques in [8, 9] exploit temporal
correlation in the signal by incorporating signal dynamics within the sparse
recovery problem. Here, each temporal snapshot can be interpreted as a block
within a high dimensional sparse signal. In [10], sparse recovery algorithms
were developed to exploit known probabilistic information on the support of
the sparse signal. These probabilities can potentially be derived from the
correlated blocks within BCS. Sequentially solving the sub-problems in BCS by
exploiting block correlations can enhance reconstruction performance, albeit
at the cost of losing the ability to solve the sub-problems in parallel.

In this paper, we study the fundamental limits of BCS and discuss how to
exploit correlation across different blocks for sparse recovery. First, we
derive a Welch bound-based guarantee on the reconstruction error with BCS.
Then, we develop a data-driven approach to learn the correlation across
blocks. The learned correlation is then used for sparse recovery in our serial
BCS method. Finally, we consider a near-field channel estimation problem to
show that our data-driven serial BCS technique outperforms standard BCS at a
moderate SNR.

Notation: $a\in\mathbb{R}$ denotes a scalar, $\textbf{a}\in\mathbb{R}^{n}$ a
vector, $\textbf{A}\in\mathbb{R}^{n_{1}\times n_{2}}$ a matrix and
${\boldsymbol{\mathcal{A}}}\in\mathbb{R}^{n_{1}\times n_{2}\times\dots\times
n_{d}}$ a $d$-th order tensor. The convolution of ${\boldsymbol{\mathcal{A}}}$
and ${\boldsymbol{\mathcal{B}}}\in\mathbb{R}^{m_{1}\times
m_{2}\times\dots\times m_{d}}$ is denoted by
${\boldsymbol{\mathcal{C}}}={\boldsymbol{\mathcal{A}}}\circledast{\boldsymbol{\mathcal{B}}}\in\mathbb{R}^{k_{1}\times
k_{2}\times\dots\times k_{d}}$ where
$k_{i}=\max(n_{i},m_{i})\;i\in\\{1,2,\dots,d\\}$, such that its entry
${\boldsymbol{\mathcal{C}}}(i_{1},i_{2},\ldots,i_{d})=\sum_{j_{1}=1}^{n_{1}}\sum_{j_{2}=1}^{n_{2}}\ldots\sum_{j_{d}=1}^{n_{d}}{\boldsymbol{\mathcal{A}}}(j_{1},j_{2},\ldots,j_{d}){\boldsymbol{\mathcal{B}}}(i_{1}-j_{1},i_{2}-j_{2},\ldots,i_{d}-j_{d})$.
$\langle\mathbf{a},\mathbf{b}\rangle$ is the inner product and
$\|\mathbf{a}\|_{2}$ is the $\ell_{2}$ norm of $\mathbf{a}$. $[N]$ denotes the
set $\\{1,2,\cdots N\\}$.

##  II Analysis of block compressed sensing

###  II-A Motivation for BCS

In classical CS, a sparse signal $\textbf{x}\in\mathbb{R}^{n}$ is projected on
an $m$-dimensional space using a CS matrix $\textbf{A}\in\mathbb{R}^{m\times
n}$ with $m\ll n$. The projections are usually perturbed by noise
$\textbf{v}\in\mathbb{R}^{m}$, where each entry has a variance of
$\sigma^{2}$. The vector containing the noisy compressed representation of x
is defined as

$\textbf{y}=\textbf{A}\textbf{x}+\textbf{v}.$ (1)

CS algorithms estimate a sparse vector that best explains the observed
measurements. Computationally tractable algorithms based on $\ell_{1}$ norm
minimization and greedy methods have been developed to solve for x from (1).
These algorithms can find a sparse solution when $\mathbf{x}$ is sufficiently
sparse and when the CS matrix $\mathbf{A}$ satisfies the restricted isometry
property or has a low mutual coherence [11].

The computational complexity of CS algorithms in [12, 13, 14] is
$\mathcal{O}(mn)$ for $o(1)$ sparse signals, due to matrix-vector
multiplications with $\mathbf{A}$ and $\mathbf{A}^{T}$. For high dimensional
signals, both $n$ and $m$ can be significant so that real-time sparse recovery
becomes infeasible under practical constraints.

###  II-B Storage and computational complexity of BCS

BCS uses a block-diagonal structure for the CS matrix to alleviate the heavy
storage demanded by generic CS matrices. We define $\beta$ as the number of
blocks or partitions in BCS, and assume that the blocks are of equal size. To
model discontiguous blocks in BCS, we define a permutation matrix
$\boldsymbol{\Pi}\in\mathbb{R}^{n\times n}$. The CS matrix in BCS is then

$\mathbf{A}=\begin{bmatrix}\textbf{A}_{1}&\textbf{0}&\dots&\textbf{0}\\\
\textbf{0}&\textbf{A}_{2}&\ddots&\vdots\\\ \vdots&\ddots&\ddots&\textbf{0}\\\
\textbf{0}&\dots&\textbf{0}&\textbf{A}_{\beta}\end{bmatrix}\boldsymbol{\Pi},$
(2)

where $\mathbf{A}_{b}\in\mathbb{R}^{\frac{m}{\beta}\times\frac{n}{\beta}}$.
Setting $\boldsymbol{\Pi}$ to an identity matrix I corresponds to the standard
BCS scenario where the signal is partitioned into contiguous blocks. When the
BCS matrix in (2) is used in (1), we observe that any measurement is a
projection of a block within the sparse signal. We define
$\textbf{z}=\boldsymbol{\Pi}\textbf{x}$. The $b^{\mathrm{th}}$ block of
$\boldsymbol{\Pi}\textbf{x}$, defined as $\textbf{z}_{b}$, is a vector that
contains entries of $\boldsymbol{\Pi}\textbf{x}$ indexed from $(b-1)n/\beta+1$
to $bn/\beta$. The dimension of the each signal block is $n/\beta$, and
$m/\beta$ compressed measurements are acquired per block.

In the BCS model, we define $\textbf{y}_{b}\in\mathbb{R}^{\frac{m}{\beta}}$ as
the measurements of the $b^{\mathrm{th}}$ block and $\textbf{v}_{b}$ as the
associated noise. Then,
$\textbf{y}=(\textbf{y}_{1};\textbf{y}_{2};\cdots;\textbf{y}_{\beta})$ is a
column stacked version of the block measurements and
$\textbf{v}=(\textbf{v}_{1};\textbf{v}_{2};\cdots;\textbf{v}_{\beta})$. The
measurements in (1) can be split in BCS as

$\textbf{y}_{b}=\textbf{A}_{b}\textbf{z}_{b}+\textbf{v}_{b}\;\;\forall
b\in\\{1,2,\cdots,\beta\\}.$ (3)

To estimate x, a straightforward approach is to first solve for
$\\{\textbf{z}_{b}\\}_{b=1}^{\beta}$ from (3). These problems can be solved in
parallel. Then, the reconstructed blocks can be stacked together to obtain z.
Finally, the permuntation operation in $\textbf{z}=\boldsymbol{\Pi}\textbf{x}$
can be inverted to estimate x.

Due to the block structure in BCS matrices, it is sufficient to store
$\mathcal{O}(\beta{mn}/{\beta^{2}})$ entries of (2). This requirement is
$\beta$ times lower compared to the standard CS approach, which necessitates
storing $\mathcal{O}({mn})$ entries. Furthermore, the complexity of matrix-
vector multiplications in solving for one block in (3) is
$\mathcal{O}({mn}/{\beta^{2}})$, when compared to $\mathcal{O}(mn)$ associated
with (1). When the $\beta$ blocks are independently recovered, BCS requires a
computational complexity that is $1/\beta^{2}$ lower than that of standard CS.
When the $\beta$ blocks are estimated sequentially, BCS has a computational
complexity that is $1/\beta$ lower than standard CS.

###  II-C Mutual coherence-based error bounds for BCS

In this section, we derive fundamental limits on parallel BCS-based recovery
using the notion of mutual coherence.

The mutual coherence of $\textbf{A}\in\mathbb{R}^{m\times n}$ is defined as
[11]

$\mu(\textbf{A})=\max_{(k,\ell):k\neq\ell}\frac{|{\langle\textbf{a}_{k},\textbf{a}_{\ell}\rangle}|}{{\|\textbf{a}_{k}\|_{2}\|\textbf{a}_{\ell}\|_{2}}},$
(4)

where $\textbf{a}_{k}$ is the $k^{\mathrm{th}}$ column of $\mathbf{A}$. A
small mutual coherence allows for better sparse reconstruction with the
orthogonal matching pursuit (OMP)[11]. For a CS matrix of size $m\times n$,
the smallest possible mutual coherence that can be achieved is given by the
Welch bound [15], i.e.,

$\mu(\textbf{A})\geq\sqrt{\frac{n-m}{m(n-1)}}.$ (5)

The bound in (5) only depends on the size of the CS matrix.

We derive a new lower bound on the mutual coherence of the BCS matrix in (2).
To this end, we first notice that permuting the columns of A does not impact
the mutual coherence and assume $\boldsymbol{\Pi}=\textbf{I}$ in our analysis.
Next, we observe that the columns of A corresponding to different blocks are
orthogonal. For instance, the first column of A in (2) is orthogonal to all
the columns indexed from $n/\beta+1$ to $n$, due to the block structure.
Therefore, for a BCS matrix, the mutual coherence in (4) is not determined by
a pair of columns chosen from two different blocks. Applying the definition in
(4) for the BCS matrix in (2) results in

$\displaystyle\mu(\textbf{A})$
$\displaystyle=\max\\{\mu(\textbf{A}_{1}),\mu(\textbf{A}_{2}),\cdots,\mu(\textbf{A}_{\beta})\\}.$
(6)

As each of the $\beta$ mutual coherences in (6) are lower bounded by the Welch
bound associated with an $m/\beta\times n/\beta$ matrix, we can write that

$\displaystyle\mu(\textbf{A})\geq\sqrt{\frac{\frac{n}{\beta}-\frac{m}{\beta}}{\frac{m}{\beta}(\frac{n}{\beta}-1)}}=\sqrt{\frac{n-m}{m(\frac{n}{\beta}-1)}}.$
(7)

For $\beta>1$, we observe that our lower bound in (7) for BCS is larger than
the Welch bound in (5). The key takeaway here is that BCS matrices cannot
achieve the Welch bound in (5) due to the block diagonal constraint on the CS
matrix. For $\beta$ partitions, the smallest possible mutual coherence that
BCS matrices can achieve is given by our bound in (7).

The coherence bound in (7) can be used to find limits on the mean squared
error (MSE) with the OMP algorithm. For an $s$-sparse signal, this error bound
is given by [11]

$\|\hat{\textbf{x}}-\textbf{x}\|_{2}^{2}\leq\frac{2(1+\alpha)}{(1-(s-1)\mu)^{2}}s\sigma^{2}\log{m},$
(8)

where $\alpha$ is a constant independent of $m$ and $n$. We notice from (8)
that a small mutual coherence results in a tight upper bound on the
reconstruction error. When a BCS matrix in (2) is used for sparse recovery
with the OMP, the tightest upper bound on the MSE is obtained when the mutual
coherence is equal to the lower limit in (7). A sketch of the MSE bound for
OMP-based BCS is shown in Fig. 1 for different $\beta$. We observe that the
error bound monotonically increases with $\beta$, indicating that partitioning
the sparse recovery problem can lead to poor reconstruction.

Figure 1: The plot shows the upper bound on the MSE in (8) achieved by our
coherence bound in (2). The bound increases with $\beta$, the number of
blocks. Here, $s=50$, $n=10^{4}$, $\alpha=0.5$ and $\sigma=10^{-2}$.

##  III Data-driven serial BCS

To alleviate the MSE loss due to partitioning in BCS, we exploit correlation
across different blocks within the sparse signal during reconstruction. We
demonstrate our approach by modeling and exploiting these correlations, using
an example of a clustered sparse signal illustrated in Fig. 2.

A good partitioning strategy is key to the success of both standard BCS and
our data-driven serial BCS. For instance, contiguous partitioning, shown in
Fig. 2a, that groups contiguous indices into a block does not work well for
clustered sparse signals. This is because the non-zero coefficients are likely
to be concentrated in a few blocks, resulting in poor reconstruction when an
equal number of CS measurements are allocated per block. Comb-like periodic
partitioning, shown in Fig. 2b, is better suited to recover clustered sparse
signals through BCS as it results in almost the same number of non-zero
coefficients per block. In this paper, we consider clustered sparse signals
because spatial wireless channels are known to exhibit such property in the
Fourier dictionary [16]. Furthermore, comb-like partitions can be realized in
hardware using the codes designed in [17]. We observe from Fig. 2b that a non-
zero entry reconstructed in a comb-like block (e.g., entry at $\Delta$)
provides useful side-information on the support of its neighbours for
clustered sparse signals. This information can be used as a prior to
reconstruct other blocks.

(a) Contiguous partitioning

(b) Comb-like partitioning

Figure 2: Examples of partitioning strategies for $\beta=4$. The non-zero
components of the 2D sparse signal here are shaded. For clustered sparse
signals partitioned with a comb-like pattern, reconstructing one block
provides side information on the support of other blocks.

Now, we discuss an offline algorithm to learn the side-information provided by
a non-zero component on the support of its neighbours. Our algorithm considers
clustered sparse signals and comb-like partitions. Using a dataset of sparse
signals $\\{{\boldsymbol{\mathcal{X}}}^{(j)}\\}_{j=1}^{J}$, our algorithm
estimates the probability that an immediate neighbour is non-zero conditioned
on an entry being non-zero. This probability is determined empirically by
counting the number of samples within the dataset in which the neighbouring
entry is non-zero. For a 2D signal, a naive approach finds $\sim 8n$ such
probabilities as there are $8$ neighbours around each entry except for those
at the edges. For a $d^{\mathrm{th}}$ order tensor, this number increases to
$\sim(3^{d}-1)n$. As storing all these probabilities consumes significant
memory, we assume that the probability that a neighbour is non-zero is
spatially invariant. Under this assumption, our algorithm only learns the
support correlation kernel with $3^{d}-1$ entries for a $d^{\mathrm{th}}$
order tensor. Our procedure to learn this kernel
${\boldsymbol{\mathcal{\boldsymbol{\Theta}}}}$ from a dataset of sparse
signals is summarized in Algorithm 1. We explain this algorithm for the tensor
case as we consider a fourth order tensor estimation problem in our
simulations.

Algorithm 1 Construction of the support correlation kernel

Input: Sparse signals $\\{{\boldsymbol{\mathcal{X}}}^{(j)}\\}^{J}_{j=1}$ in
$\mathbb{R}^{n_{1}\times n_{2}\times\dots\times n_{d}}$.

Define: Set $\mathcal{S}_{j}$ has indices of non-zero entries in
${\boldsymbol{\mathcal{X}}}^{(j)}$.

for $j:1\rightarrow J$ do

Initialize $d^{\mathrm{th}}$ order kernel
${\boldsymbol{\mathcal{\boldsymbol{\kappa}}}}=\mathbf{0}\in\mathbb{R}^{3\times
3\times\cdots(d\,\mathrm{times})}$.

for index $(\omega_{1},\omega_{2},\cdots,\omega_{d})$ in $\mathcal{S}_{j}$ do

Find immediate neighbours around this index:

$\Omega=\\{(\zeta_{1},\zeta_{2},\cdots,\zeta_{d}):|\omega_{i}-\zeta_{i}|\leq
1\forall i\in[d]\\}$.

Scan through all the entries of ${\boldsymbol{\mathcal{X}}}^{(j)}$ at $\Omega$

If ${\boldsymbol{\mathcal{X}}}^{(j)}(\zeta_{1},\zeta_{2},\cdots,\zeta_{d})\neq
0:$

Add $1$ to
${\boldsymbol{\mathcal{\boldsymbol{\kappa}}}}(\zeta_{1}-\omega_{1},\zeta_{2}-\omega_{2},\cdots,\zeta_{d}-\omega_{d})$.

Set ${\boldsymbol{\mathcal{\boldsymbol{\kappa}}}}(0,0,\cdots,0)=0$.

end for

${\boldsymbol{\mathcal{\boldsymbol{\Theta}}}}^{(j)}={\boldsymbol{\mathcal{\boldsymbol{\kappa}}}}/{\text{cardinality}(\mathcal{S}_{j})}$

end for

return
${\boldsymbol{\mathcal{\boldsymbol{\Theta}}}}=\sum_{j=1}^{J}{\boldsymbol{\mathcal{\boldsymbol{\Theta}}}}^{(j)}/J$

Figure 3: Our proposed data-driven serial BCS method first reconstructs one of
the signal blocks (partitions). Then, the reconstructed signal together with
the learned support correlation ${\boldsymbol{\mathcal{\boldsymbol{\Theta}}}}$
is used to find a prior on the support associated with the other blocks. This
prior is used together with the BCS measurements to reconstruct the next
block, and the process is repeated until all the blocks are recovered.

We discuss notation involved in our data-driven BCS method to recover a sparse
tensor ${\boldsymbol{\mathcal{X}}}\in\mathbb{R}^{n_{1}\times
n_{2}\times\dots\times n_{d}}$ from its compressed representation. Here, the
dimension of the signal is $n=n_{1}n_{2}\cdots n_{d}$. Assuming that the
$\ell^{\mathrm{th}}$ mode is partitioned by a factor of $\beta_{\ell}$, the
total number of partitions is $\beta=\beta_{1}\beta_{2}\cdots\beta_{d}$. The
tensor
${\boldsymbol{\mathcal{X}}}_{b}\in\mathbb{R}^{\frac{n_{1}}{\beta_{1}}\times\frac{n_{2}}{\beta_{2}}\times\dots\times\frac{n_{d}}{\beta_{d}}}$
comprises the entries of ${\boldsymbol{\mathcal{X}}}$ at the indices in
partition $b$. Further,
${\boldsymbol{\mathcal{X}}}_{b,\mathrm{ext}}\in\mathbb{R}^{n_{1}\times
n_{2}\times\dots\times n_{d}}$ has the entries of ${\boldsymbol{\mathcal{X}}}$
at the indices in partition $b$ and zeros at all other indices. Finally, we
define ${\boldsymbol{\mathcal{P}}}\in\mathbb{R}^{n_{1}\times
n_{2}\times\dots\times n_{d}}$ such that
${\boldsymbol{\mathcal{P}}}(\zeta_{1},\zeta_{2},\cdots,\zeta_{d})$ is
proportional to the probability that
${\boldsymbol{\mathcal{X}}}(\zeta_{1},\zeta_{2},\cdots,\zeta_{d})$ is non-
zero. After each stage of our algorithm that solves for a block within
${\boldsymbol{\mathcal{X}}}$, ${\boldsymbol{\mathcal{P}}}$ is updated using
the reconstructed block and the support correlation kernel
${\boldsymbol{\mathcal{\boldsymbol{\Theta}}}}$.

Our data-driven BCS technique, summarized in Algorithm 2, first solves for one
of the $\beta$ blocks within the sparse signal from its compressed
representation. Without loss of generality, our algorithm first reconstructs
${\boldsymbol{\mathcal{X}}}_{1}$, i.e., block $1$, from

$\textbf{y}_{1}=\mathsf{A}_{1}({\boldsymbol{\mathcal{X}}}_{1})+\mathbf{v}_{1},$
(9)

where $\mathsf{A}_{1}(\cdot)$ denotes a linear compression operator. To solve
for this first block, we set ${\boldsymbol{\mathcal{P}}}=s\mathbf{1}/n$, where
$s$ is the average sparsity level estimated from the dataset. We use the
logit-weighted OMP (LW-OMP) algorithm [10], that exploits the support prior
${\boldsymbol{\mathcal{P}}}$ together with the compressed measurements, to
reconstruct a signal. The reconstructed block
$\hat{{\boldsymbol{\mathcal{X}}}_{1}}$ is extended with zeros to obtain
$\hat{{\boldsymbol{\mathcal{X}}}}_{1,\mathrm{ext}}$, which is then convolved
with ${\boldsymbol{\mathcal{\boldsymbol{\Theta}}}}$ to update the support
prior. The updated prior ${\boldsymbol{\mathcal{P}}}$ is then used to
determine the block for which the most side-information on its support is
available, i.e., the block for which
$\mathrm{sum}\\{\mathrm{vec}({\boldsymbol{\mathcal{P}}}_{b})\\}$ is maximum.
This block is recovered next and the procedure is repeated until all the
blocks within ${\boldsymbol{\mathcal{X}}}$ are recovered. Our algorithm for
the matrix case ($d=2$) is illustrated in Fig. 3.

Algorithm 2 Proposed data-driven serial BCS algorithm

Input: BCS measurements $\\{\textbf{y}_{b}\\}^{\beta}_{b=1}$, CS operators
$\\{\mathsf{A}_{b}\\}^{\beta}_{b=1}$,

Learned support correlation kernel
${\boldsymbol{\mathcal{\boldsymbol{\Theta}}}}$.

Define: Set $\nu=[\beta]$, Support prior tensor
${\boldsymbol{\mathcal{P}}}=s\mathbf{1}/n$ ($s$ is

the average sparsity level), $\hat{{\boldsymbol{\mathcal{X}}}}=\mathbf{0}$ and
$b=1$.

for $j:1\rightarrow\beta$ do

$\hat{{\boldsymbol{\mathcal{X}}}}_{b}=\text{LW-
OMP}(\textbf{y}_{b},\mathsf{A}_{b},{\boldsymbol{\mathcal{P}}}_{b},\sigma^{2})$
$\\#\mathrm{Solveforblock}b$

$\hat{{\boldsymbol{\mathcal{X}}}}=\hat{{\boldsymbol{\mathcal{X}}}}+\hat{{\boldsymbol{\mathcal{X}}}}_{b,\mathrm{ext}}$

${\boldsymbol{\mathcal{P}}}={\boldsymbol{\mathcal{P}}}+|\hat{{\boldsymbol{\mathcal{X}}}}_{b,\mathrm{ext}}\circledast{\boldsymbol{\mathcal{\boldsymbol{\Theta}}}}|$
$\\#$Update prior using estimate

$\nu\leftarrow\nu\setminus\\{b\\}$

$b=\underset{k\in\nu}{\arg\max}[\,\mathrm{sum}\\{\mathrm{vec}({\boldsymbol{\mathcal{P}}}_{k})\\}]$
$\\#$ Next block to solve

end for

return $\hat{{\boldsymbol{\mathcal{X}}}}$

Our data-driven serial BCS with LW-OMP exploits structure across partitions
but incurs higher computational complexity than standard BCS with OMP. This is
because our method serially solves for the $\beta$ partitions, unlike standard
BCS which can be parallelized. Further, support prior update in our method
requires convolution, which adds a complexity of $\mathcal{O}(3^{d}n/\beta)$
per partition. Considering $\beta$ partitions, the increase in complexity with
our method over serial BCS is $\mathcal{O}(3^{d}n)$, still modest relative to
unpartitioned standard CS.

##  IV Simulation results

We consider near-field spatial channel estimation between a $16\times 16$
transmitter (TX) and a $8\times 8$ receiver (RX) at a carrier frequency of
$300\,\mathrm{GHz}$. The TX and the RX use half-wavelength spaced uniform
planar arrays. The RX is placed on a plane at a distance of $30\,\mathrm{cm}$
from the TX. We generate several channel realizations according to the
propagation model in [18], by translating the RX along the plane and also
rotating it at random. The kernel
${\boldsymbol{\mathcal{\boldsymbol{\Theta}}}}$ in our method was learned with
Algorithm 1 using these realizations. Each channel realization is a complex-
valued tensor in $\mathbb{C}^{16\times 16\times 8\times 8}$. We use
${\boldsymbol{\mathcal{X}}}\in\mathbb{C}^{16\times 16\times 8\times 8}$ to
denote the 4D-discrete Fourier transform (4D-DFT) of the channel. As the
channel exhibits clustered sparsity in the 4D-DFT,
${\boldsymbol{\mathcal{X}}}$ is a sparse tensor.

Figure 4: NMSE with SNR for different partitions $\beta$. Here, the
subsampling ratio is $40\%$. Our method outperforms standard BCS in the
moderate SNR regime, at the expense of extra computations.

The measurements in BCS were acquired by applying 2D-codes in [17] at both the
TX and the RX. These codes facilitate partitioning by implementing
complementary comb-like masks over the sparse signal. By adjusting the
periodicity within the comb, the number of partitions $\beta$ can be
configured. In our simulations, we consider $\beta\in\\{1,16,64\\}$ to compare
the performance of the proposed approach against standard BCS. Note that
$\beta=1$ corresponds to the unpartitioned problem, which is solved using the
standard OMP algorithm. Our approach employs the LW-OMP [10] within Algorithm
2, whereas standard BCS uses the classical OMP algorithm to solve for each
partition. The transmit power is scaled to achieve the desired SNR in the
compressed measurements.

Figure 5: For all the methods, NMSE monotonically decreases with the
subsampling ratio for each $\beta$. Here, the SNR was set to
$30\,\mathrm{dB}$. Our data-driven serial BCS method outperforms standard BCS.

To evaluate the performance of BCS algorithms, we use the normalized mean
squared error (NMSE) between the sparse channel and its estimate, i.e.,
$\mathbb{E}[\|{\boldsymbol{\mathcal{X}}}-\hat{{\boldsymbol{\mathcal{X}}}}\|_{\mathrm{F}}^{2}]/\mathbb{E}[\|{\boldsymbol{\mathcal{X}}}\|_{\mathrm{F}}^{2}]$,
where the expectation is taken across all channel realizations. The dimension
of the sparse vector solved for $\beta=1$ is $16\times 16\times\times 8\times
8=16384$. With BCS, however, this dimension reduces to $16384/\beta$ per
partition. From Fig. 4, we observe that BCS yields higher NMSE than standard
CS (equivalent to BCS with $\beta=1$). Moreover, the NMSE deteriorates for
increasing $\beta$. These findings align with our analysis in Sec. II-C, where
higher values of $\beta$ led to an increased error bound.

Number of blocks $\beta$ | Standard BCS | Proposed method
---|---|---
1 | 18.2 $\times$ $10^{8}$ ms | -
16 | 5910 ms | 9827 ms
64 | 476 ms | 1469 ms
TABLE I: Computation time with standard BCS and the proposed method. BCS with
$\beta=1$ is same as classical CS.

We notice from Fig. 4 that the proposed data-driven serial BCS outperforms
standard BCS in the moderate SNR regime. At low SNR, the performance of both
the methods is almost the same. This is because of the noise in the
reconstruction, which impacts the support prior estimated, i.e.,
${\boldsymbol{\mathcal{P}}}$, in Algorithm 2. As the support prior estimate is
not reliable at low SNR, it does not improve the reconstruction with LW-OMP.
At high SNR, the measurements in each block provide reliable information to
determine the support even without any prior information. Next, we observe
from Fig. 5 that the proposed method outperforms standard BCS for different
subsampling ratios, i.e., $m/n$, at an SNR of $30\,\mathrm{dB}$. This
performance improvement comes at the expense of an increased computational
complexity compared to standard BCS, owing to the support prior calculation in
our approach. The increase in complexity, however, is small when compared to
the complexity of the unpartitioned CS problem for $\beta=1$. The computation
times of all the methods, on a desktop computer, for $n=16384$ and a
subsampling ratio of $\sim 40\%$ is summarized in Table I.

##  V Conclusions

In this paper, we studied block compressed sensing (BCS) for high-dimensional
sparse recovery. We proved that the lower bound on the mutual coherence of a
BCS matrix is higher than the Welch bound associated with a standard CS matrix
of the same dimensions. We also proposed a data-driven serial BCS method,
which uses reconstructed blocks to estimate priors on the support of
subsequent blocks. Using simulations for the channel estimation problem, we
showed that our method outperforms standard BCS at moderate SNR.

## References

  * [1] G. Calisesi, A. Ghezzi, D. Ancora, C. D’Andrea, G. Valentini, A. Farina, and A. Bassi, “Compressed sensing in fluorescence microscopy,” _Progress in Biophysics and Molec. Biol._ , vol. 168, pp. 66–80, 2022.
  * [2] W. L. Chan, K. Charan, D. Takhar, K. F. Kelly, R. G. Baraniuk, and D. M. Mittleman, “A single-pixel terahertz imaging system based on compressed sensing,” _Applied Physics Letters_ , vol. 93, no. 12, 2008.
  * [3] S. Nie and I. F. Akyildiz, “Deep kernel learning-based channel estimation in ultra-massive MIMO communications at 0.06-10 THz,” in _2019 IEEE Globecom Workshops (GC Wkshps)_ , 2019, pp. 1–6.
  * [4] L. Gan, “Block compressed sensing of natural images,” in _2007 15th Intl. Conf. on Digital Signal Process._ , 2007, pp. 403–406.
  * [5] B. Zhang, Y. Liu, J. Zhuang, K. Wang, and Y. Cao, “Matrix permutation meets block compressed sensing,” _Journal of Visual Communication and Image Representation_ , vol. 60, pp. 69–78, 2019.
  * [6] J. Chen, X. Zhang, and H. Meng, “Self-adaptive sampling rate assignment and image reconstruction via combination of structured sparsity and non-local total variation priors,” _Digital Signal Process._ , vol. 29, pp. 54–66, 2014.
  * [7] R. Pournaghshband and M. Modarres-Hashemi, “A novel block compressive sensing algorithm for sar image formation,” _Signal Processing_ , vol. 210, p. 109053, 2023.
  * [8] J. Ziniel and P. Schniter, “Dynamic compressive sensing of time-varying signals via approximate message passing,” _IEEE Transactions on Signal Processing_ , vol. 61, no. 21, pp. 5270–5284, 2013.
  * [9] M. S. Asif and J. Romberg, “Dynamic updating for sparse time varying signals,” in _Annual Conf. on Inform. Sciences and Sys._ , 2009, pp. 3–8.
  * [10] J. Scarlett, J. Evans, and S. Dey, “Compressed sensing with prior information: Information-theoretic limits and practical decoders,” _IEEE Transactions on Signal Processing_ , vol. 61, p. 427, 01 2013.
  * [11] Z. Ben-Haim, Y. C. Eldar, and M. Elad, “Coherence-based performance guarantees for estimating a sparse vector under random noise,” _IEEE Trans. on Signal Process._ , vol. 58, no. 10, pp. 5030–5043, 2010.
  * [12] J. A. Tropp and A. C. Gilbert, “Signal recovery from random measurements via orthogonal matching pursuit,” _IEEE Trans. on Inform. theory_ , vol. 53, no. 12, pp. 4655–4666, 2007.
  * [13] T. Blumensath and M. E. Davies, “Iterative hard thresholding for compressed sensing,” _Appl. and comput. harmonic anal._ , vol. 27, no. 3, pp. 265–274, 2009.
  * [14] D. L. Donoho, A. Maleki, and A. Montanari, “Message-passing algorithms for compressed sensing,” _Proceedings of the National Academy of Sciences_ , vol. 106, no. 45, pp. 18 914–18 919, 2009.
  * [15] L. Welch, “Lower bounds on the maximum cross correlation of signals,” _IEEE Trans. on Inform. theory_ , vol. 20, no. 3, pp. 397–399, 1974.
  * [16] C. K. Anjinappa, Y. Zhou, Y. Yapici, D. Baron, and I. Guvenc, “Channel estimation in mmWave hybrid MIMO system via off-grid dirichlet kernels,” in _IEEE Global Commun. Conf. (GLOBECOM)_ , 2019, pp. 1–6.
  * [17] H. Masoumi, M. Verhaegen, and N. J. Myers, “In-sector compressive beam alignment for mmWave and THz radios,” _arXiv preprint arXiv:2308.13268_ , 2023.
  * [18] E. Torkildson, U. Madhow, and M. Rodwell, “Indoor millimeter wave MIMO: Feasibility and performance,” _IEEE Trans. on Wireless Commun._ , vol. 10, pp. 4150–4160, 12 2011.


# Rates of Convergence in the Central Limit Theorem for Markov Chains, with an
Application to TD Learning

R. Srikant
ECE and CSL, UIUC
<EMAIL_ADDRESS>
###### Abstract

We prove a non-asymptotic central limit theorem for vector-valued martingale
differences using Stein’s method, and use Poisson’s equation to extend the
result to functions of Markov Chains. We then show that these results can be
applied to establish a non-asymptotic central limit theorem for Temporal
Difference (TD) learning with averaging.

##  1 Introduction

TD learning is the most widely studied algorithm for evaluating the
performance of a given policy from data in Markov Decision Processes [1]. The
convergence of TD learning with decaying stepsizes was proved in [2]. Over the
last few years, there has been a resurgent interest in understanding the non-
asymptotic convergence behavior of TD learning and further using it to study
other other reinforcement algorithms which use TD learning within their
framework. The first result of this type was obtained in [3] who used a
projection step in their version of the TD learning algorithm. The first
finite-time bounds for unprojected TD learning were obtained in [4], who
studied both fixed step-sizes and decaying step-sizes. TD learning can be
viewed as a special case of linear stochastic approximation with Markovian and
multiplicative noise. The results in [4] were extended to nonlinear,
contractive stochastic approximation algorithms in [5], who also derived
finite-time performance bounds for specific choices of decaying step-size
rules, which also allowed them to study the robustness of different step-size
choices. Finite-time performance bounds for variants of TD learning and other
schemes have been studied in [6, 7, 8].

While all of the above papers study standard TD learning, in practice, it is
well-known that averaging the iterates can help improve the asymptotic
variance. This idea was originally studied for stochastic approximation and
stochastic gradient descent in the seminal works in [9, 10, 11]; see [12] for
a finite-time analysis. More recently, the asymptotic behavior of TD learning
and other algorithms with averaging has been studied in [13, 14, 15]. As in
the papers on stochastic approximation, the conclusion in [13] is that, for a
suitable choice of decaying stepsizes, averaging the output leads to efficient
estimators in the sense that the variance is minimized. Further, [13] also
established functional central limit theorems and central limit theorems for
TD learning and other reinforcement learning algorithms. The work in [14]
shows order-optimal rate of convergence for the mean-squared error, while [15]
obtains higher moment bounds and high probability deviation bounds. In
contemporaneous work to this paper, finite-time performance bounds on the
covariance of the averaged TD learning iterates have been obtained in [16]
using more robust step-size rules. Their model is more general; they study
two-time scale linear stochastic approximation algorithms, with TD learning
with averaging being a special case of the models that they study. A central
limit theorem for two time-scale stochastic approximation has also been
established for two time-scale stochastic approximation in [17]. Related work
also includes a central limit theorem for Polyak-Ruppert averaged Q-learning
[18].

Our goal in this paper is to obtain bounds on the rate of convergence to
normality to complement the central limit theorem in [13]. To this end, we
first establish a result for martingale differences by combining the Lindeberg
decomposition along with Stein’s method, which was first used in [19] for
scalar valued martingale differences. This key idea in [19] was extended to
the case of vector-valued martingales in [20] who also used it to obtain a
non-asymptotic version of the central limit theorem in [11]. Our approach
builds upon the ideas in [19, 20], but we also introduce new techniques in
order to obtain stronger results and to deal with the differences in our model
compare to [20]:

  * •

Unlike [20], we would like our results to hold when the deviation from
normality is measured using the Wasserstein distance. The results in [20]
provide bounds only for weaker distances. Therefore, we have to obtain results
for more general test functions than the ones considered in [20]. To do so, we
leverage the results on the regularity of solutions to Stein’s equation for
vector-valued random variables in [21, 22], which makes our analysis
considerably different.

  * •

Since our main intended application is TD learning, our noise sequence is
Markovian. Our approach is to convert the Markovian noise into a martingale
difference noise using Poisson’s equation; see [23, 24, 25], for example. To
the best of our knowledge, this idea has not been used to derive rates of
convergence for the central limit theorem for functions of Markov chains.

  * •

In addition to the Markovian nature of our noise, we have to deal with the
fact that the noise multiplies the TD learning parameters, unlike the case of
stochastic gradient descent, where the martingale difference noise is
additive.

A notable difference between our results and the results in [19, 20] for
martingales is the fact that we explicitly relate our convergence bounds to
the rate of convergence of conditional covariances of martingale differences
to the asymptotic covariance (see item (2) in Assumption 2 and Theorem 1).
This is important to us since our goal is relate our results on TD learning to
the asymptotic covariance results in [13]. Therefore, we cannot use the
assumptions on the variance in [19, Theorem 1] or the covariance assumption in
[20, Equation 3.1]. It is also worth remarking that Stein’s method to derive
rates of convergence for the martingale central limit theorem is relatively
new. The early papers on obtaining rates of convergence for the martingale
central limit theorem use other techniques; see, for example, [26] and other
references in [19]. Here, we use Stein’s method because we found it easier to
quantify the impact of the asymptotic variance on the rate of convergence
using Stein’s method. Specifically, we show that the regularity results in
[21, 22] and the calculations in [19, 21] can be used, along with insights
from Markov chain theory, to understand the impact of the asymptotic
covariance term on the rates of convergence in the martingale and Markov chain
central limit theorems.

Our main results in the paper are as follows:

  * •

We derive non-asymptotic central limit theorems for martingales and Markov
chains by exploiting ideas from [19, 21, 22]. One key difference compared to
[19, 20] is that we write our results explicitly in terms of the asymptotic
covariance of the martingale, which is useful when we study Markov chains.

  * •

We then extend the results to Markov chains using Poisson’s equation to relate
Markov chains to martingales. In this case, the asymptotic covariance can be
characterized as in well known in the literature; see [25], for example. The
error bound on the deviation from normality in the case of martingale
difference sequences appears to be complicated. Interestingly however, under
appropriate assumptions, the expressions simplify remarkably in the case of
Markov chains and leads to a simple convergence rate result.

  * •

Finally, we apply the Markov chain results, along with the ideas in [11], to
estimate the rate of convergence in the central limit theorem for TD learning
with Polyak-Ruppert averaging. An advantage of the approach here is that our
result directly establishes a central limit theorem for the scaled averaged of
the iterates, instead of first establishing a function central limit theorem
for the entire trajectory of the iterates and then establishing the central
limit theorem.

It is important to note that, unlike [11], our result only holds for averaging
along with an appropriate decaying step-size rule. Specifically, the result
does not hold for fixed step-size TD learning with averaging due to the bias
issue pointed out in [13, 27, 28]); see also related observations [29, 30] for
other problems with Markov noise. A nonasymptotic central limit theorem has
been established for constant step-size stochastic gradient descent in [31],
but the model considered here is different from the problem in [31] (for
example, Markovian noise, constant step-size vs decaying step-size) and
therefore, the techniques there do not seem to apply to our problem.

The rest of the paper is organized as follows. In Section 2, we establish an
upper bound on the rate of convergence in the martingale central limit theorem
under some assumptions, and in Section 3, we apply these results to bound the
rate of convergence in the central limit theorem for functions of Markov
chains using Poisson’s equation. The Markov chain results are then applied to
study the central limit for averaged iterates of TD leaning in Section 4.
Concluding remarks and suggestions for future work are provided in Section 5.

##  2 Martingale Central Limit Theorem

In this section, we present a bound on the rate of convergence for the
martingale central limit theorem in terms of the Wasserstein distance using
Stein’s method.

###  2.1 Preliminaries

In this subsection, we will describe the notation we use to establish the rate
of convergence in the central limit theorem for vector-valued martingale
differences. Let $\\{m_{k}\\}_{k\geq 1}$ be a $d$-dimensional martingale
difference sequence with respect to a filtration $\\{\mathcal{F}_{k}\\}_{k\geq
0}$ i.e., $m_{k}$ is $\mathcal{F}_{k}$-adapted, $E(||m_{k}||)<\infty,$ where
$||\cdot||$ is the standard Euclidean norm in $\Re^{d}$ and

$E(m_{k}|\mathcal{F}_{k-1})=0\qquad\forall k\geq 1.$

######  Assumption 1.

We assume that the martingale difference sequence satisfies the following
properties:

  1. 1.

$E(||m_{k}||^{2+\beta})$ exists a.s. for all $k\geq 1$ and some
$\beta\in(0,1).$

  2. 2.

$\Sigma_{k}=E(m_{k}m_{k}^{T}|\mathcal{F}_{k-1})$ exists and further assume
that

$\lim_{n\rightarrow\infty}\frac{\Sigma_{1}+\ldots+\Sigma_{n}}{n}=\Sigma_{\infty}$

almost surely for some $\Sigma_{\infty}\succ 0,$ where $\succ 0$ indicates
positive definiteness.

$\diamond$

Item (1) in the assumption is motivated by the fact that we want to extend the
techniques in [21] to martingales and item (2) is motivated by the fact that
one can show that $\Sigma_{\infty}$ exists for certain martingales derived
from some Markov chains, as we will see in the next section.

Define $S_{n}=\sum_{k=1}^{n}m_{k}.$ We are interested in whether the following
normalized sum

$W_{n}:=\frac{S_{n}}{\sqrt{n}}$ (1)

converges to $N(0,\Sigma_{\infty})$, and if so, we are interested in the rate
of convergence. We will use Stein’s method to establish such a convergence and
to estimate the rate of convergence [32, 33].

The Wasserstein distance between two $d$-dimensional random vectors $X$ and
$Y$ is defined as

$d_{\mathcal{W}}(X,Y)=\sup_{h\in Lip_{1}}E[h(X)-h(Y)],$

where $Lip_{L}$, the class of Lipschitz functions from $\Re^{d}$ to $\Re$ with
Lipschitz constant $L,$ is defined as

$Lip_{L}=\\{h:|h(x)-h(y)|\leq L||x-y||\\}.$

For later reference, we use the notation $||A||_{op}$ to denote the operator
norm of a matrix $A,$ i.e.,

$||A||_{op}=\sup_{x:||x||=1}||Ax||.$

The following results from [21, Propositions 2.2 and 2.3] and [22, Theorem
3.1] will be useful to us and so we present them together as a lemma for
future reference.

######  Lemma 1.

Let $g\in Lip_{L}$ and $f$ be the solution to Stein’s equation:

$g(u)-E(g(Z))=\Delta f(u)-u^{T}\nabla f(u),$ (2)

where $\Delta$ and $\nabla$ are the Laplacian and gradient operators,
respectively, and $Z\sim N(0,I)$. Then, $f$ has the following regularity
properties for any $\beta\in(0,1)$:

  1. 1.

$||\nabla^{2}f(x)-\nabla^{2}f(y)||_{op}\leq\tilde{C}_{1}(d,\beta)||x-y||^{\beta}L,$
where $\nabla^{2}$ is the Hessian operator.

  2. 2.

$|\Delta f(x)-\Delta f(y)|\leq\tilde{C_{2}}(d,\beta)||x-y||^{\beta}L$, where

$\tilde{C}_{1}=C_{1}(d)+\frac{2}{1-\beta},\quad\tilde{C}_{2}=C_{2}(d)+\frac{2d}{1-\beta},\quad
C_{1}=2^{3/2}\frac{1+2d\Gamma((1+d)/2)}{d\Gamma(d/2)},\quad
C_{2}=2\sqrt{\frac{2d}{\pi}}.$

  3. 3.

$\langle\nabla^{2}f(x),u_{1}u_{2}^{T}\rangle_{HS}\leq
CL||u_{1}||\cdot||u_{2}||\,\forall x,u_{1},u_{2}\in\Re^{d},$

where $C$ is a constant (independent of $d$ and $\beta$) and
$\langle\cdot,\cdot\rangle_{HS}$ is the Hilbert-Schmidt inner product (also
called the Frobenius inner product), i.e., the sum of the element-by-element
product of two matrices.

$\diamond$

Next, we state a well-known result regarding Ornstein-Uhlenbeck processes
which will be useful later to establish the nonasymptotic version of the
martingale central limit theorem (for example, the result can be easily
inferred from the discussions on generators and linear stochastic differential
equations in [34]).

######  Lemma 2.

Consider the O-U process

$dX_{t}=AX_{t}dt+BdW_{t},$

where $W$ is the standard $d$-dimensional Wiener process, and $A$ and $B$ are
$d\times d$ matrices. Then, the generator of the process $\mathcal{A}$ is the
operator given by

$\mathcal{A}f(x):=\lim_{\delta\rightarrow
0}E\left(\frac{X_{t+\delta}-X_{t}|X_{t}=x}{\delta}\right)=\nabla^{T}f(x)Ax+\frac{1}{2}Tr(\nabla^{2}f(x)BB^{T}),$

where $f:\Re^{d}\rightarrow\Re$ is any twice differentiable function. If the
eigenvalues of $A$ have strictly negative real parts, then the stationary
distribution is Gaussian with mean zero and covariance $\Sigma,$ which
satisfies the Lyapunov equation

$\Sigma A^{T}+A\Sigma+BB^{T}=0.$

Thus, if $\tilde{Z}_{k}\sim N(0,\Sigma),$ then
$E(\mathcal{A}f(\tilde{Z}_{k}))=0$ for any twice differentiable function $f.$
$\diamond$

It is worth noting that Stein’s equation (2) is motivated by the above lemma
when $A=-I$ and $BB^{T}=2I,$ in which case the right-hand side of (2) is the
same as $\mathcal{A}f(u)$. Thus, if $u$ in (2) is replaced by a random vector,
then the left-hand side of (2) is a measure of the difference between the
distribution of that random vector and a standard Gaussian distribution while
the right-hand side will hopefully be small if the distribution of the random
vector is close to a standard Gaussian distribution (see, for example, [35,
36, 37]).

###  2.2 Rate of Convergence

We now state and prove the main result of this section.

######  Theorem 1.

Let $Z\sim N(0,I)$ and $W_{n}$ be defined as in (1). Then, under Assumption 1,
we have

$\displaystyle d_{\mathcal{W}}(W_{n},\Sigma_{\infty}^{1/2}Z)$
$\displaystyle\leq\frac{1}{\sqrt{n}}\sum_{k=1}^{n}\Biggl{[}\frac{\tilde{C}_{1}(d,\beta)||\Sigma_{\infty}||_{op}}{(n-k+1)^{(1+\beta)/2}}E\biggl{(}||\Sigma_{\infty}^{-1/2}m_{k}||^{2+\beta}\biggr{)}+\frac{\tilde{C}_{2}(d,\beta)||\Sigma_{\infty}||_{op}}{(n-k+1)^{(1+\beta)/2}}E\biggl{(}||\Sigma_{\infty}^{-1/2}m_{k}||^{\beta}\biggr{)}$
$\displaystyle\quad-\frac{1}{n-k+1}Tr\Bigl{(}A_{k}\Bigl{(}\Sigma_{\infty}^{-1/2}E(\Sigma_{k})\Sigma_{\infty}^{-1/2}-I\Bigr{)}\Bigr{)}\Biggr{]},$

where $A_{k}$ is a matrix with the property

$||A_{k}||_{op}\leq C(n-k+1)||\Sigma_{\infty}^{1/2}||_{op},$

where $C,\tilde{C}_{1},\tilde{C}_{2}$ are the constants in statement (3) of
Lemma 1.

###### Proof.

We first note that

$d_{\mathcal{W}}(W_{n},\Sigma_{\infty}^{1/2}Z)=\sup_{h\in
Lip_{1}}E[h(W_{n})-h(\Sigma_{\infty}^{1/2}Z)]=\sup_{h\in
Lip_{1}}\frac{1}{\sqrt{n}}E[h(S_{n})-h(T_{1})],$ (3)

where $T_{k}=\Sigma_{\infty}^{1/2}(Z_{k}+\ldots+Z_{n}),$ where $Z_{i}$ are
i.i.d. $N(0,I)$ random vectors, independent of the martingale difference
sequence. Then, using the Lindeberg decomposition [19], we have

$E\Bigl{(}h(S_{n})-h(T_{1})\Bigr{)}=\sum_{k=1}^{n}E\Bigl{(}h(S_{k}+T_{k+1})-h(S_{k-1}+T_{k})\Bigr{)},$
(4)

where we have defined $T_{n+1}=S_{0}=0.$

The rest of the proof uses the vector extensions of the techniques used in
[19] for scalar random variables. As in well known in the Stein’s method
literature, going from scalar-valued random variables to random vectors is
non-trivial. To deal with random vectors, we use the ideas in [21, Theorem
3.1] which were developed for i.i.d. random variables. One of our
contributions is to show that these ideas can be used for martingales using
the Lindeberg decomposition. While doing so, we also have to deal with the
asymptotic covariance $\Sigma_{\infty},$ which does not appear in [19, 20] due
to their assumptions, but will be important to us when we study Markov chains.

Motivated by the form of each term in the Lindeberg decomposition, we consider
the difference

$h(x)-E(h(a+QZ)),$

where $x,a\in\Re^{d}$ and $Q\in\Re^{d\times d}$ is a symmetric positive
definite matrix. The idea is to study this expression and then later, after
appropriate conditioning, substitute

$x=S_{k}+T_{k+1},\quad a=S_{k-1},\quad\mathrm{and}\quad
Q=\sqrt{n-k+1}\Sigma_{\infty}^{1/2}.$

Note that the choice of $Q$ is motivated by the fact that $T_{k}$ has the same
distribution as

$\sqrt{n-k+1}\Sigma_{\infty}^{1/2}Z.$

Next define the function $\tilde{h}$ by $\tilde{h}(u)=h(a+Qu).$ Since $h\in
Lip_{1},$ $\tilde{h}\in Lip_{L}$ with $L=||Q||_{op}.$ Further, from Lemma 1,
we know that the solution $f$ to the Stein equation

$\tilde{h}(u)-E(\tilde{h}(Z))=\Delta f(u)-u^{T}\nabla f(u)$

has nice regularity properties. Note that Stein’s equation can be rewritten in
terms of $h$ as follows:

$\displaystyle h(x)-E(h(a+QZ))$ $\displaystyle=$
$\displaystyle\tilde{h}(Q^{-1}(x-a))-E(\tilde{h}(Z))$ $\displaystyle=$
$\displaystyle\Delta f(Q^{-1}(x-a))-(Q^{-1}(x-a))^{T}\nabla f(Q^{-1}(x-a)).$

Thus,

$\displaystyle
E\Bigl{(}h(S_{k}+T_{k+1})-h(S_{k-1}+T_{k})\Bigr{|}\mathcal{F}_{k}\Bigr{)}$
$\displaystyle=E\biggl{(}\Delta
f\Bigl{(}\frac{\Sigma_{\infty}^{-1/2}}{\sqrt{n-k+1}}(m_{k}+T_{k+1})\Bigr{)}$
$\displaystyle-(m_{k}+T_{k+1})^{T}\frac{\Sigma_{\infty}^{-1/2}}{\sqrt{n-k+1}}\nabla
f\Bigl{(}\frac{\Sigma_{\infty}^{-1/2}}{\sqrt{n-k+1}}(m_{k}+T_{k+1})\Bigr{)}|\mathcal{F}_{k}\biggr{)}$
$\displaystyle=E\biggl{(}\Delta
f\Bigl{(}\frac{\Sigma_{\infty}^{-1/2}}{\sqrt{n-k+1}}m_{k}+\tilde{Z}_{k}\Bigr{)}-(\frac{\Sigma_{\infty}^{-1/2}}{\sqrt{n-k+1}}m_{k}+\tilde{Z}_{k})^{T}\nabla
f\Bigl{(}\frac{\Sigma_{\infty}^{-1/2}}{\sqrt{n-k+1}}m_{k}+\tilde{Z}_{k}\Bigr{)}|\mathcal{F}_{k}\biggr{)},$

where $\tilde{Z}_{k}\sim N(0,\sqrt{\frac{n-k}{n-k+1}}I)$ has the same
distribution as $\frac{1}{\sqrt{n-k+1}}\Sigma_{\infty}^{-1/2}T_{k+1}$. Next,
we rewrite the expression as a sum of four terms as follows:

$\displaystyle
E\Bigl{(}h(S_{k}+T_{k+1})-h(S_{k-1}+T_{k})\Bigr{|}\mathcal{F}_{k-1}\Bigr{)}$
$\displaystyle=E\biggl{(}\frac{n-k}{n-k+1}\Delta
f\Bigl{(}\frac{\Sigma_{\infty}^{-1/2}}{\sqrt{n-k+1}}m_{k}+\tilde{Z}_{k}\Bigr{)}-\tilde{Z}_{k}^{T}\nabla
f\Bigl{(}\frac{\Sigma_{\infty}^{-1/2}}{\sqrt{n-k+1}}m_{k}+\tilde{Z}_{k}\Bigr{)}|\mathcal{F}_{k-1}\biggr{)}$
$\displaystyle+E\biggl{(}\frac{1}{n-k+1}\Delta
f\Bigl{(}\frac{\Sigma_{\infty}^{-1/2}}{\sqrt{n-k+1}}m_{k}+\tilde{Z}_{k}\Bigr{)}|\mathcal{F}_{k-1}\biggr{)}$
$\displaystyle-E\Bigl{(}(\frac{\Sigma_{\infty}^{-1/2}}{\sqrt{n-k+1}}m_{k})^{T}\nabla
f(\tilde{Z}_{k})|\mathcal{F}_{k-1}\Bigr{)}$
$\displaystyle-E\biggl{(}(\frac{\Sigma_{\infty}^{-1/2}}{\sqrt{n-k+1}}m_{k})^{T}\Bigl{(}\nabla
f\Bigl{(}\frac{\Sigma_{\infty}^{-1/2}}{\sqrt{n-k+1}}m_{k}+\tilde{Z}_{k}\Bigr{)}-\nabla
f(\tilde{Z}_{k})\Bigr{)}|\mathcal{F}_{k-1}\biggr{)}$

The first term on the right-hand side of the previous equation is zero by
applying Lemma 2 with

$A=-I,\qquad BB^{T}=2\frac{n-k}{n-k+1}I,$

and using the fact that $\Delta f=Tr(\nabla^{2}f).$ The third term is zero by
the martingale difference property. Next, let $\Theta\sim Uniform[0,1]$ be a
random variable, independent of all other random variables, and using the
fundamental theorem of calculus, we have

$\displaystyle
E\Bigl{(}h(S_{k}+T_{k+1})-h(S_{k-1}+T_{k})\Bigr{|}\mathcal{F}_{k-1}\Bigr{)}$
$\displaystyle=E\biggl{(}\frac{1}{n-k+1}\Delta
f\Bigl{(}\frac{\Sigma_{\infty}^{-1/2}}{\sqrt{n-k+1}}m_{k}+\tilde{Z}_{k}\Bigr{)}|\mathcal{F}_{k-1}\biggr{)}$
$\displaystyle-E\biggl{(}(\frac{\Sigma_{\infty}^{-1/2}}{\sqrt{n-k+1}}m_{k})^{T}\Bigl{[}\nabla^{2}f\Bigl{(}\Theta\frac{\Sigma_{\infty}^{-1/2}}{\sqrt{n-k+1}}m_{k}+\tilde{Z}_{k}\Bigr{)}-\nabla^{2}f(\tilde{Z}_{k})\Bigr{]}(\frac{\Sigma_{\infty}^{-1/2}}{\sqrt{n-k+1}}m_{k})|\mathcal{F}_{k-1}\biggr{)}$
$\displaystyle-E\biggl{(}(\frac{\Sigma_{\infty}^{-1/2}}{\sqrt{n-k+1}}m_{k})^{T}\nabla^{2}f(\tilde{Z}_{k})(\frac{\Sigma_{\infty}^{-1/2}}{\sqrt{n-k+1}}m_{k})|\mathcal{F}_{k-1}\biggr{)}.$

Next, we perform the following sequence of manipulations: first, we apply
statement (1) in Lemma 1 to the second term above, use the definition of
$\Sigma_{k}$ (in Assumption 1) in the fourth term, use the fact that $\Delta
f=Tr(\nabla^{2}f),$ and finally apply statement (2) in Lemma 1, while noting
that $|\Theta|\leq 1$ a.s., to get

$\displaystyle
E\Bigl{(}h(S_{k}+T_{k+1})-h(S_{k-1}+T_{k})\Bigr{|}\mathcal{F}_{k-1}\Bigr{)}$
$\displaystyle=E\biggl{(}\frac{1}{n-k+1}\Delta
f\Bigl{(}\frac{\Sigma_{\infty}^{-1/2}}{\sqrt{n-k+1}}m_{k}+\tilde{Z}_{k}\Bigr{)}|\mathcal{F}_{k-1}\biggr{)}$
$\displaystyle-\frac{1}{n-k+1}Tr\Bigl{(}E(\nabla^{2}f(\tilde{Z}_{k}))\Sigma_{\infty}^{-1/2}\Sigma_{k}\Sigma_{\infty}^{-1/2}\Bigr{)}$
$\displaystyle\qquad+\frac{\tilde{C}_{1}(d,\beta)||\Sigma_{\infty}||_{op}}{(n-k+1)^{(1+\beta)/2}}E\biggl{(}||\Sigma_{\infty}^{-1/2}m_{k}||^{2+\beta}|\mathcal{F}_{k-1}\biggr{)}$
$\displaystyle=E\biggl{(}\frac{1}{n-k+1}\Bigl{(}\Delta
f\Bigl{(}\frac{\Sigma_{\infty}^{-1/2}}{\sqrt{n-k+1}}m_{k}+\tilde{Z}_{k}\Bigr{)}-\Delta
f(\tilde{Z}_{k})\Bigr{)}|\mathcal{F}_{k-1}\biggr{)}$
$\displaystyle\quad-\frac{1}{n-k+1}Tr\Bigl{(}\Sigma_{\infty}^{-1/2}\Bigl{(}\Sigma_{k}-\Sigma_{\infty}\Bigr{)}\Sigma_{\infty}^{-1/2}E(\nabla^{2}f(\tilde{Z}_{k}))\Bigr{)}$
$\displaystyle+\frac{\tilde{C}_{1}(d,\beta)||\Sigma_{\infty}||_{op}}{(n-k+1)^{(1+\beta)/2}}E\biggl{(}||\Sigma_{\infty}^{-1/2}m_{k}||^{2+\beta}|\mathcal{F}_{k-1}\biggr{)}$
$\displaystyle=\frac{\tilde{C}_{2}(d,\beta)||\Sigma_{\infty}||_{op}}{(n-k+1)^{(1+\beta)/2}}E\biggl{(}||\Sigma_{\infty}^{-1/2}m_{k}||^{\beta}|\mathcal{F}_{k-1}\biggr{)}$
$\displaystyle\quad-\frac{1}{n-k+1}Tr\Bigl{(}\Sigma_{\infty}^{-1/2}\Bigl{(}\Sigma_{k}-\Sigma_{\infty}\Bigr{)}\Sigma_{\infty}^{-1/2}E(\nabla^{2}f(\tilde{Z}_{k}))\Bigr{)}$
$\displaystyle+\frac{\tilde{C}_{1}(d,\beta)||\Sigma_{\infty}||_{op}}{(n-k+1)^{(1+\beta)/2}}E\biggl{(}||\Sigma_{\infty}^{-1/2}m_{k}||^{2+\beta}|\mathcal{F}_{k-1}\biggr{)}.$

We obtain the result by taking an expectation to remove the conditioning,
defining $A_{k}$ to be $E(\nabla^{2}f(\tilde{Z})),$ then substituting the
above expression in (4) and (3). To obtain the upper bound on $A_{k},$ we use
statement (3) in Lemma 1 and the fact that

$||A_{k}||_{op}=\sup_{||u_{1}||=||u_{2}||=1}|\langle
A_{k},u_{1}u_{2}^{T}\rangle_{HS}|.$

∎

The above theorem is useful if one can show that the upper bound on the
Wasserstein distance goes to zero as $n\rightarrow\infty,$ in which case, the
rate at which the upper bound goes to zero provides some information on the
rate of convergence. In the next section, we will show that such a result can
be obtained for a certain martingale derived from the solution of Poisson’s
equation for Markov chains and obtain a bound on the rate of convergence for
the central limit theorem for functions of Markov chains.

##  3 Markov Chain Central Limit Theorem

Consider a time-homogeneous Markov chain $\\{X_{k}\\}_{k\geq 0}$ on a state-
space $\mathcal{S}$ with some initial distribution $\eta.$ For simplicity, we
will assume that the state space is discrete, although that assumption may be
relaxed by using the theory for general state spaces in [25, 38] for some of
our results. Let $r(X_{k})\in\Re^{d}$ denote a vector of rewards accrued when
in state $X_{k}.$ We are interested in rates of convergence to normality of

$\frac{\sum_{k=1}^{n}r(X_{k})}{\sqrt{n}}.$

As is common in the literature, we will study the above scaled sum by relating
the reward process to a martingale via Poisson’s equation [25, 39, 40, 41].

######  Assumption 2.

We make the following assumptions about the Markov chain.

  1. 1.

The Markov chain has a unique stationary distribution $\pi$ and is uniformly
geometrically ergodic, i.e., there exists $K_{3}>0$ and $\rho\in[0,1)$ such
that $\forall x\in\mathcal{S},$ we have

$||P^{n}(x,y)-\pi(y)||_{TV}\leq K_{1}\rho^{n}.$

  2. 2.

There exists a solution $V:S\rightarrow\Re^{d}$ to the following Poisson’s
equation:

$\bar{r}=r(x)+E(V(X_{1})|X_{0}=x)-V(x)\qquad\forall x\in S,$ (5)

where $\bar{r}=E_{y\sim\pi}(r(y)),$ such that $E(||V(X_{k})||)\leq
M_{V},\forall k\geq 0$ for some $M\leq\infty.$

  3. 3.

Define $m_{k}=E(V(X_{k})|X_{k-1})-V(X_{k})$ for $k\geq 1.$ Clearly,
$E(m_{k}|\mathcal{F}_{k-1})=0$ from Poisson’s equation. Further, we assume
that $\\{m_{k}\\}_{k\geq 1}$ is a martingale difference sequence with respect
to the filtration $\\{\mathcal{F}_{k}\\}_{k\geq 0}$ satisfying Assumption 1,
with $\mathcal{F}_{k}=\sigma(X_{0},X_{1},X_{2},\ldots,X_{k})$ and
$\Sigma_{\infty}$ given by

$\Sigma_{\infty}=\sum_{i,j\in S}\pi_{i}P_{ij}(V(j)-V(i))(V(j)-V(i))^{T},$ (6)

where $P_{ij}$ is the $(i,j)^{\rm th}$ element of the probability transition
matrix.

$\diamond$

Clearly, under Assumption 2 item (2), there exists constants $M_{2+\beta}$ and
$M_{\beta}$ such that

$E(||m_{k}||^{2+\beta})\leq M_{2+\beta}\mbox{ and }E(||m_{k}||^{\beta})\leq
M_{\beta}\qquad\forall k\geq 1.$

Assumption 2 item (2) on the solution to Poisson’s equation is made for
simplicity. It can be weakened considerably as one can see by following the
proof of the subsequent theorem. The assumption here would be sufficient for
our application to TD learning in the next section, so we choose to make this
assumption.

We now state the main result of this section, and note that the constants
hidden in the $\mathcal{O}$-notation can be made explicit by applying the
bounds assumed in Assumption 2 to the right-hand side of the bound in Theorem
1 as outlined in the proof of the theorem. However, this would lead to messy
expressions.

######  Theorem 2.

Let

$U_{n}=\frac{1}{\sqrt{n}}\sum_{k=1}^{n}(r(X_{k-1})-\bar{r}).$

Then, under Assumption 2,

$d_{\mathcal{W}}(U_{n},\Sigma_{\infty}^{1/2}Z)=\mathcal{O}\Bigl{(}\frac{1}{n^{\beta/2}}\Bigr{)},$

where $Z\sim N(0,I).$

###### Proof.

We start with the following standard use of Poisson’s equation and the
definition of $m_{k}$ in Assumption 2 (see, for example, [25]) to write

$\displaystyle r(X_{k-1})-\bar{r}$
$\displaystyle=V(X_{k-1})-E(V(X_{k})|X_{k-1})$
$\displaystyle=V(X_{k-1})-V(X_{k})+V(X_{k})-E(V(X_{k})|X_{k-1})$
$\displaystyle=V(X_{k-1})-V(X_{k})+m_{k}$

Thus,

$U_{n}=W_{n}+\frac{V(X_{0})-V(X_{n})}{\sqrt{n}},\quad\mathrm{where}\quad
W_{n}=\frac{1}{\sqrt{n}}\sum_{k=1}^{n}m_{k}.$

For any $h\in Lip_{1},$

$h(U_{n})-h(Z)\leq
h(W_{n})-h(Z)+||U_{n}-W_{n}||=h(W_{n})-h(Z)+\frac{2M_{V}}{\sqrt{n}}.$

This implies

$d_{\mathcal{W}}(U_{n},\Sigma_{\infty}^{1/2}Z)\leq
d_{\mathcal{W}}(W_{n},\Sigma_{\infty}^{1/2}Z)+\frac{2M_{V}}{\sqrt{n}}.$

We can now apply Theorem 1 to bound
$d_{\mathcal{W}}(W_{n},\Sigma_{\infty}^{1/2}Z).$ By Assumption 2, the first
two terms in the bound in Theorem 1 go to zero at rate

$\frac{1}{\sqrt{n}}\sum_{k=1}^{n}\frac{1}{(n-k+1)^{(1+\beta)/2}}=\frac{1}{\sqrt{n}}\sum_{j=1}^{n}\frac{1}{j^{(1+\beta)/2}}\leq\frac{K_{2}}{n^{\beta/2}},$

for some universal constant $K_{2}.$ For the third term in the bound in
Theorem 1,

$\displaystyle\frac{1}{n-k+1}Tr\Bigl{(}A_{k}\Sigma_{\infty}^{-1/2}\Bigl{(}E(\Sigma_{k})-\Sigma_{\infty}\Bigr{)}\Sigma_{\infty}^{-1/2}\Bigr{)}$
$\displaystyle\leq\frac{1}{n-k+1}||A_{k}||_{HS}||\Sigma_{\infty}^{-1/2}E(\Sigma_{k})\Sigma_{\infty}^{-1/2}-I||_{HS}$
$\displaystyle\leq\frac{1}{n-k+1}\sqrt{d}||A_{k}||_{op}||\Sigma_{\infty}^{-1/2}E(\Sigma_{k})\Sigma_{\infty}^{-1/2}-I||_{HS}$
$\displaystyle\leq
C\sqrt{d}||\Sigma_{\infty}^{1/2}||_{op}||\Sigma_{\infty}^{-1/2}E(\Sigma_{k})\Sigma_{\infty}^{-1/2}-I||_{HS},$
(7)

where the last step follows from Theorem 1.

Next, letting $p_{i}(k)=Prob(X_{k}=i),$ we get

$\displaystyle E(\Sigma_{k})$ $\displaystyle=\sum_{i,j\in
S}p_{i}(k)P_{ij}(V(j)-V(i))(V(j)-V(i))^{T}$ $\displaystyle=\sum_{i,j\in
S}(p_{i}(k)-\pi_{i})P_{ij}(V(j)-V(i))(V(j)-V(i))^{T}+\Sigma_{\infty},$

which implies

$||\Sigma_{\infty}^{-1/2}E(\Sigma_{k})\Sigma_{\infty}^{-1/2}-I||_{HS}\leq||\Sigma_{\infty}^{-1/2}||_{HS}^{2}\sum_{i,j\in
S}|p_{i}(k)-\pi_{i}|\cdot||P_{ij}(V(j)-V(i))(V(j)-V(i))^{T}||_{HS}.$

By our assumption of uniform geometric ergodicity, the above expression goes
to zero geometrically fast and thus,

$\frac{1}{\sqrt{n}}\sum_{k=1}^{n}\frac{1}{n-k+1}Tr\Bigl{(}A_{k}\Sigma_{\infty}^{-1/2}\Bigl{(}E(\Sigma_{k})-\Sigma_{\infty}\Bigr{)}\Sigma_{\infty}^{-1/2}\Bigr{)}\leq\frac{K_{3}}{\sqrt{n}},$

where $K_{3}$ is a problem-dependent constant. ∎

Our proof relied on assumptions about the solution to Poisson’s equation (5).
The Poisson equation is the same as the average-reward relative value function
equation, and there are simple conditions for the existence of a solution for
this equation in the case of finite state-space Markov chains [42]. Examples
of models with more general state-spaces where one can obtain solutions to
Poisson’s equation satisying our assumptions can be found in [43].

As in [21], we can obtain the following corollary by choosing $\beta=1-2/\log
n$ and noting that $\tilde{C}_{1}$ and $\tilde{C}_{2}$ become $O(\log n)$ in
that case.

######  Corollary 1.

Under Assumption 3, the following rate of convergence holds:

$d_{\mathcal{W}}(\sqrt{n}(U_{n},\Sigma_{\infty}^{1/2}Z)=\mathcal{O}\Bigl{(}\frac{\log
n}{\sqrt{n}}\Bigr{)},$

where $\Sigma_{\infty}$ has been defined in Assumption 2, with $r(X_{k})$ as
in Assumption 3. $\diamond$

It is worth noting that, if we set $\beta=1-2/\log n,$ we require third
moments to exist for the non-asymptotic martingale central limit theorem to
hold. However, this condition is easily satisfied if we assume the Markov
chain has a finite state space and is irreducible and aperiodic, as we will do
so in the next section. If one were to apply Theorems 1 and 2 to more general
scenarios, then the third moment condition must be verified or assumed.

##  4 An Application to TD Learning

Temporal Difference (TD) Learning is a common method to learn the performance
of a policy in reinforcement learning; see [44, 1] for details. In [2], it was
shown that TD learning can be represented as a jump Markovian linear dynamical
system as follows:

$\theta_{k+1}=\theta_{k}-\epsilon_{k}(A(X_{k})\theta_{k}+b(X_{k})),$ (8)

where $\theta_{k}\in\Re^{d},$ $\epsilon_{k}$ is the stepsize,
$\\{X_{k}\\}_{k\geq 0}$ is a finite state-space Markov chain, $A(X_{k})$ is a
$d\times d$ matrix, $b(X_{k})$ is a $d$-dimensional vector and $\theta_{0}$
being some constant vector. We make the following assumptions.

######  Assumption 3.

We assume that the TD learning algorithms satisfies the following conditions:

  1. 1.

$\\{X_{k}\\}$ is an irreducible, aperioidic, finite state Markov chain. So it
satisfies Assumption 2, item (1).

  2. 2.

Let $\bar{A}:=E_{X_{k}\sim\pi}(A(X_{k})).$ Then, $-\bar{A}$ is a Hurwitz
matrix, i.e., all eigenvalues have strictly negative real parts.

$\diamond$

Define $\theta^{*}$ by $\bar{A}\theta^{*}+\bar{b}=0,$ where
$\bar{b}=E_{X_{k}\sim\pi}(b(X_{k})).$ Here, we are interested in the Polyak-
Ruppert averaged version of the TD-learning algorithm, i.e.,

$\bar{\theta}_{n}=\frac{1}{n}\sum_{k=1}^{n}\theta_{k}.$

Our main result of this section characterizes the rate of convergence of this
algorithm for an appropriately chosen $\epsilon$. As in the previous theorem,
we have chosen to write the result using the $\mathcal{O}$ notation rather
than making the constants explicit, but these constants can be made explicit
by following the details the proof of the theorem.

######  Theorem 3.

Under Assumption 3, and with $\epsilon_{k}=1/(k+1)^{\delta},$
$\delta\in(0.5,1)$ the following rate of convergence holds:

$d_{\mathcal{W}}(\sqrt{n}(\bar{\theta}_{n}-\theta^{*}),(\bar{A}^{-1}\Sigma_{\infty}\bar{A}^{-T})^{1/2}Z)=\mathcal{O}\Bigl{(}\frac{1}{n^{\delta-0.5}}\Bigr{)},$

where $\Sigma_{\infty}$ has been defined in Assumption 2, with $r(X_{k})$ as
in Assumption 3.

###### Proof.

Let $\Delta_{k}=\theta_{k}-\theta^{*}.$ It follows from (8) that

$\displaystyle\Delta_{k+1}=$
$\displaystyle(I-\epsilon_{k}\bar{A})\Delta_{k}-\epsilon_{k}(A(X_{k})-\bar{A})\Delta_{k}-\epsilon_{k}\Bigl{(}(A(X_{k})-\bar{A})\theta^{*}+(b(X_{k})-\bar{b})\Bigr{)}.$

Thus,

$\displaystyle\sqrt{n}(\bar{\theta}_{n}-\theta^{*})$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{k=1}^{n}\prod_{l=0}^{k-1}(I-\epsilon_{l}\bar{A})\Delta_{0}$
$\displaystyle-\frac{1}{\sqrt{n}}\sum_{k=1}^{n}\sum_{j=0}^{k-1}\prod_{l=j}^{k-1}(I-\epsilon_{l}\bar{A})\epsilon_{j}(A(X_{j})-\bar{A})\Delta_{j}$
$\displaystyle-\frac{1}{\sqrt{n}}\sum_{k=1}^{n}\sum_{j=0}^{k-1}\prod_{l=j}^{k-1}(I-\epsilon_{l}\bar{A})\epsilon_{j}\Bigl{(}(A(X_{j})-\bar{A})\theta^{*}+(b(X_{j})-\bar{b})\Bigr{)}$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{k=1}^{n}\prod_{l=0}^{k-1}(I-\epsilon_{l}\bar{A})\Delta_{0}$
$\displaystyle-\frac{1}{\sqrt{n}}\sum_{j=0}^{n-1}\epsilon_{j}\sum_{k=j+1}^{n}\prod_{l=j}^{k-1}(I-\epsilon_{l}\bar{A})(A(X_{j})-\bar{A})\Delta_{j}$
$\displaystyle-\frac{1}{\sqrt{n}}\sum_{j=0}^{n-1}\epsilon_{j}\sum_{k=j+1}^{n}\prod_{l=j}^{k-1}(I-\epsilon_{l}\bar{A})\Bigl{(}(A(X_{j})-\bar{A})\theta^{*}+(b(X_{j})-\bar{b})\Bigr{)}$

Define

$\Upsilon_{j}=\epsilon_{j}\sum_{k=j+1}^{n}\prod_{l=j}^{k-1}(I-\epsilon_{l}\bar{A})-\bar{A}^{-1}.$

Then, we have

$\displaystyle\sqrt{n}(\bar{\theta}_{n}-\theta^{*})$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{k=1}^{n}\prod_{l=0}^{k-1}(I-\epsilon_{l}\bar{A})\Delta_{0}$
(9)
$\displaystyle-\frac{1}{\sqrt{n}}\sum_{j=0}^{n-1}\bar{A}^{-1}(A(X_{j})-\bar{A})\Delta_{j}$
(10)
$\displaystyle-\frac{1}{\sqrt{n}}\sum_{j=0}^{n-1}\Upsilon_{j}(A(X_{j})-\bar{A})\Delta_{j}$
(11)
$\displaystyle-\frac{1}{\sqrt{n}}\sum_{j=0}^{n-1}\Upsilon_{j}\Bigl{(}(A(X_{j})-\bar{A})\theta^{*}+(b(X_{j})-\bar{b})\Bigr{)}$
(12)
$\displaystyle-\frac{1}{\sqrt{n}}\sum_{j=0}^{n-1}\bar{A}^{-1}\Bigl{(}(A(X_{j})-\bar{A})\theta^{*}+(b(X_{j})-\bar{b})\Bigr{)}$
(13)

Now, we will complete the proof in several steps by considering each of the
terms above.

#### Step 1: Equation (9).

Since

$\left|\left|\sum_{k=1}^{n}\prod_{l=0}^{k-1}(I-\epsilon_{l}\bar{A})\epsilon_{j}\right|\right|\leq||\Upsilon_{0}||_{op}+||\bar{A}^{-1}||_{op},$
(14)

and, from [11], $\Upsilon_{0}$ is bounded independent of $n$, we conclude that
this term is $\mathcal{O}(\frac{1}{\sqrt{n}}).$

#### Step 2: Expression (10).

To analyze (10), we consider the solution $V_{A}$ to the following matrix-
valued Poisson’s equation:

$A(X_{k-1})-\bar{A}=V_{A}(X_{k-1})-E(V_{A}(X_{k})|X_{k-1}).$

Let the associated martingale difference sequence $\mu$ be defined by

$\mu_{k}=V_{A}(X_{k})-E(V_{A}(X_{k})|X_{k-1}).$

Thus, up to a sign change, (10) can be written as

$\displaystyle=\frac{1}{\sqrt{n}}\sum_{j=0}^{n-1}\bar{A}^{-1}(V_{A}(X_{j})-V_{A}(X_{j+1}))\Delta_{j}+\frac{1}{\sqrt{n}}\sum_{j=0}^{n-1}\bar{A}^{-1}\mu_{j+1}\Delta_{j}$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{j=0}^{n-1}\bar{A}^{-1}(V_{A}(X_{j})\Delta_{j}-V_{A}(X_{j+1})\Delta_{j+1})$
(15)
$\displaystyle+\frac{1}{\sqrt{n}}\sum_{j=0}^{n-1}\bar{A}^{-1}V_{A}(X_{j+1})(\Delta_{j+1}-\Delta_{j})$
(16)
$\displaystyle+\frac{1}{\sqrt{n}}\sum_{j=0}^{n-1}\bar{A}^{-1}\mu_{j+1}\Delta_{j}$
(17)

After telescoping,

$E(||(\ref{eq3})||)\leq\frac{1}{\sqrt{n}}||\bar{A}^{-1}||_{op}E(||V_{A}(X_{0})||\cdot||\Delta_{0}||+||V_{A}(X_{n})||\cdot||\Delta_{n}||))=\mathcal{O}\left(\frac{1}{\sqrt{n}}\right),$

where in the last step we have used the fact that solution to the matrix
Poisson solution is bounded due to the assumptions on our Markov chain and
$E(||\Delta(n)||^{2})$ is bounded from [4]. Next, due to the fact that our
Markov chain’s state space is finite, there exists constants $\kappa_{3}$ and
$\kappa_{4}$ such that

$\displaystyle
E(||\ref{eq4})||)\leq\frac{1}{\sqrt{n}}\sum_{j=0}^{n-1}\frac{1}{(j+1)^{\delta}}||\bar{A}^{-1}||_{op}(\kappa_{3}E(||\Delta_{j}|)|+\kappa_{4})=\mathcal{O}\left(\frac{n^{1-\delta}}{\sqrt{n}}\right)=\mathcal{O}\left(\frac{1}{n^{\delta-0.5}}\right),$

where we have again used the fact $E(||\Delta(n)||^{2})$ is uniformly bounded.
To study (17), we use the martingale difference property of $\mu$ and the fact
that it is bounded to conclude that there exists a constant $\kappa_{5}$ such
that

$\displaystyle
E(||(\ref{eq5})||^{2})\leq\frac{1}{n}\kappa_{5}\sum_{j=0}^{n-1}E\left(||\Delta_{j}||^{2}\right)=\mathcal{O}\left(\frac{n^{1.5-2\delta}}{n}\right),$

where, in the last step, we have used the fact that
$E(||\Delta_{j}||^{2})=O(1/n^{2\delta-0.5})$ [4]; see Appendix A for details.
Thus, $E(||(\ref{eq5})||)$ is $\mathcal{O}(\frac{1}{n^{\delta-0.25}}).$ It
follows that $E(||(\ref{eq2})||)$ is $\mathcal{O}(\frac{1}{n^{\delta-0.5}}).$

#### Step 3: Expression (11).

We write (11) as

$\displaystyle\frac{1}{\sqrt{n}}\sum_{j=0}^{n-1}\Upsilon_{j}(A(X_{j})-\bar{A})\Delta_{j}$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{j=0}^{n-1}\Upsilon_{j}(V_{A}(X_{j})-V_{A}(X_{j+1})+\mu_{j+1})\Delta_{j}$
(18)
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{j=0}^{n-1}(\Upsilon_{j}V_{A}(X_{j})\Delta_{j}-\Upsilon_{j+1}V_{A}(X_{j+1})\Delta_{j+1})$
(19)
$\displaystyle+\frac{1}{\sqrt{n}}\sum_{j=0}^{n-1}\left(\Upsilon_{j+1}V_{A}(X_{j+1})(\Delta_{j+1}-\Delta_{j})\right)$
(20)
$\displaystyle+\frac{1}{\sqrt{n}}\sum_{j=0}^{n-1}\left((\Upsilon_{j+1}-\Upsilon_{j})V_{A}(X_{j+1})\Delta_{j}\right)$
(21)
$\displaystyle+\frac{1}{\sqrt{n}}\sum_{j=0}^{n-1}\Upsilon_{j}\mu_{j+1}\Delta_{j}$
(22)

From [11], $||\Upsilon_{j}||_{op}$ is uniformly bounded. Therefore, all the
the terms above, except for (21), can be handled in a manner similar to the
corresponding terms in (10). If we can show that
$||\Upsilon_{j+1}-\Upsilon_{j}||_{op}$ is $\mathcal{O}(\epsilon_{j}),$ then
the analysis of (21) will be similar to the analysis of (20), and we can
conclude that $E(||(\ref{step3})||$ is also $\mathcal{O}(n^{\delta/2}).$

We can write $\Upsilon_{j+1}-\Upsilon_{j}$ as follows:

$\displaystyle\epsilon_{j+1}\sum_{k=j+2}^{n}\prod_{l=j+1}^{k-1}(I-\epsilon_{l}\bar{A})-\epsilon_{j}\sum_{k=j+1}^{n}\prod_{l=j}^{k-1}(I-\epsilon_{l}\bar{A})$
$\displaystyle=\epsilon_{j+1}\sum_{k=j+2}^{n}\prod_{l=j+1}^{k-1}(I-\epsilon_{l}\bar{A})-\epsilon_{j+1}\sum_{k=j+1}^{n}\prod_{l=j}^{k-1}(I-\epsilon_{l}\bar{A})+(\epsilon_{j+1}-\epsilon_{j})\sum_{k=j+1}^{n}\prod_{l=j}^{k-1}(I-\epsilon_{l}\bar{A})$
$\displaystyle=\epsilon_{j+1}\sum_{k=j+2}^{n}\prod_{l=j+1}^{k-1}(I-\epsilon_{l}\bar{A})-\epsilon_{j+1}\sum_{k=j+2}^{n}\prod_{l=j}^{k-1}(I-\epsilon_{l}\bar{A})$
$\displaystyle+\epsilon_{j+1}\sum_{k=j+2}^{n}\prod_{l=j}^{k-1}(I-\epsilon_{l}\bar{A})-\epsilon_{j+1}\sum_{k=j+1}^{n}\prod_{l=j}^{k-1}(I-\epsilon_{l}\bar{A})+(\epsilon_{j+1}-\epsilon_{j})\sum_{k=j+1}^{n}\prod_{l=j}^{k-1}(I-\epsilon_{l}\bar{A})$
$\displaystyle=\epsilon_{j+1}\epsilon_{j}\sum_{k=j+2}^{n}\prod_{l=j+1}^{k-1}(I-\epsilon_{l}\bar{A})\bar{A}-\epsilon_{j+1}(1-\epsilon_{j}\bar{A})+\frac{\epsilon_{j+1}-\epsilon_{j}}{\epsilon_{j}}(\Upsilon_{j}+\bar{A}^{-1})$
$\displaystyle=\epsilon_{j}(\Upsilon_{j+1}+\bar{A}^{-1})\bar{A}-\epsilon_{j+1}(1-\epsilon_{j}\bar{A})+\frac{\epsilon_{j+1}-\epsilon_{j}}{\epsilon_{j}}(\Upsilon_{j}+\bar{A}^{-1}).$

Substitute the above expression in (21) and use the following facts: (i) the
uniform boundedness $||\Upsilon_{j}||_{op}$ and $E(||\Delta_{j}||)$, (ii) the
fact that $\frac{\epsilon_{j+1}-\epsilon_{j}}{\epsilon_{j}}$ is
${o}(\epsilon_{j})$ as mentioned in [11], and (iii)
$\sum_{j=1}^{n}\epsilon_{j}=\mathcal{O}(n^{1-\delta}).$ It follows that
$E(||(\ref{eq:3})||)$ is $O(1/n^{\delta-0.5}).$

#### Step 4: Expression (12).

We note that with $r(X_{k})=A(X_{k})\theta^{*}+b(X_{k}),$ the solution of
Poisson’s equation (5) satisfies the conditions in items (2)-(4) of Assumption
2 because our Markov chain is finite-state, irreducible and aperiodic. Using
the notation in items (2)-(3) of Assumption 2, (12) can be written as

$\displaystyle\frac{1}{\sqrt{n}}\sum_{j=0}^{n-1}\Upsilon_{j}\Bigl{(}(A(X_{j})-\bar{A})\theta^{*}+(b(X_{j})-\bar{b})\Bigr{)}$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{j=0}^{n-1}\Upsilon_{j}\Bigl{(}(V(X_{j})-V(X_{j+1})+m_{j+1}))\Bigr{)}$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{j=0}^{n-1}(\Upsilon_{j}V(X_{j})-\Upsilon_{j+1}V(X_{j+1}))$
$\displaystyle+\frac{1}{\sqrt{n}}\sum_{j=0}^{n-1}(\Upsilon_{j+1}-\Upsilon_{j})V_{j+1}$
$\displaystyle+\frac{1}{\sqrt{n}}\sum_{j=0}^{n-1}\Upsilon_{j}m_{j+1}$ (23)

The uniform boundedness of $\Upsilon_{j}$ and $V(x)$ proves that the first two
terms after the last equality above go to zero at rate
$\mathcal{O}(\frac{1}{\sqrt{n}}).$ Next, using the fact that $\\{m_{j}\\}$ is
a bounded martingale difference sequence, we get

$E(\left(||\frac{1}{\sqrt{n}}\sum_{j=0}^{n-1}\Upsilon_{j}m_{j+1}||^{2}\right)=\frac{1}{n}E\left(\sum_{j=0}^{n-1}||\Upsilon_{j}m_{j+1}||^{2}\right)\leq\kappa_{6}\frac{1}{n}E\left(\sum_{j=0}^{n-1}||\Upsilon_{j}||^{2}\right)$

for some constant $\kappa_{6}.$ It is known from [11] that this term goes to
zero as $n\rightarrow\infty$. In Appendix B, we show that this term is
$\mathcal{O}(\log n/n).$

#### Conclusion of the Proof:

To finally establish the result stated in the theorem, we observe that, for
random variables $\Xi_{i},$ $i=1,2,3,4,5,$ and $h\in Lip_{1},$

$\displaystyle
E\Big{(}h(\sum_{i=1}^{5}\Xi_{i})-(h((\bar{A}^{-1}\Sigma_{\infty}\bar{A}^{-T})^{1/2}Z)\Big{)}$
$\displaystyle=\sum_{i=1}^{4}E\left(h(\sum_{k=i}^{5}\Xi_{k})-h(\sum_{j=i+1}^{5}\Xi_{j})\right)+E(h(\Xi_{5})-h((\bar{A}^{-1}\Sigma_{\infty}\bar{A}^{-T})^{1/2}Z))$
$\displaystyle\leq\sum_{i=1}^{4}E(||\Xi_{i}||)+E(h(\Xi_{5}))-h((\bar{A}^{-1}\Sigma_{\infty}\bar{A}^{-T})^{1/2}Z)$

Let $\Xi_{i},$ $i=1,2,3,4$ be the terms in (9)-(12) and $\Xi_{5}$ be (13). The
different steps of the proof have shown that each of $E(||\Xi_{i}||)$ for
$i=1,2,3,4$ is $\mathcal{O}(\frac{1}{n^{\delta-0.5}})$. Now we apply Theorem 2
to $\Xi_{5}.$ Since $\Xi_{5}$ converges to a Gaussian distribution in the
Wasserstein distance at rate $\mathcal{O}(1/n^{\beta})$ and $\beta<1$ can be
arbitrary due to our assumptions on the Markov chain, we have proved the
theorem.

∎

######  Remark 1.

We have shown an $\mathcal{O}(1/n^{\delta-0.5})$ rate of convergence for all
$\delta\in(0.5,1).$ However, this does not necessarily mean that using a
larger $\delta$ is better. To illustrate this fact, in Appendix B, we provide
bounds on the constants hidden in the $\mathcal{O}$ notation to show that they
can blow up to infinity as $\delta\rightarrow 1.$

##  5 Conclusions

We derived a non-asymptotic central limit theorem for martingales and Markov
chains, and applied the results to get rate of convergence in the central
limit theorem for TD learning with averaging. The results for martingales and
Markov chains may be more broadly applicable to study other reinforcement
learning algorithms as well as other algorithms such as nonlinear stochastic
approximation and stochastic gradient descent.

There are several possible avenues for further work:

  * •

One can try to use the martingale CLT results in this paper to strengthen the
results and allow for Markovian noise in [20]. However, these extensions would
require us to go beyond the linear setting in this paper and consider
nonlinear problems. Along similar lines, it would be also interesting to see
if our results can be used to study contractive stochastic approximation as in
[37] or other stochastic approximation schemes [45, 46].

  * •

Under some conditions, the rate of convergence in the central limit theorem
for i.i.d. random vectors is $\mathcal{O}(1/\sqrt{n})$ [22]. It would be
interesting to see if the $\mathcal{O}(\log n/\sqrt{n})$ rate of convergence
in Corollary (1) can be improved to $\mathcal{O}(1/\sqrt{n})$ under possibly
additional conditions than the ones used in this paper. The Edgeworth
expansion for scalar functions of Markov chains in [47] suggests that the
$\log n$ could possibly be eliminated for Markov chains as well. The reason
for $\log n$ in our result is due to the fact that we use the regularity
results for the solution to Stein’s equation in [21] after we do the Lindeberg
decomposition in the proof of the rate of convergence of martingale CLT. A
different approach may be required to eliminate the $\mathcal{O}(\log n)$ term
in the numerator.

  * •

The Wassertein distance used in this paper is the $1$-Wasserstein distance, it
would be interesting to obtain convergence rates in the stronger
$p$-Wasserstein distance. Work along these lines under different conditions
than in this paper can be found in [48], it would be interesting to see if
such results can be obtained for martingales and Markov chains.

  * •

It is interesting to note that, while Theorem 2 is critical to obtain
convergence to the central limit theorem for TD learning, the rate of
convergence is dominated by the parameter $\delta$ used in the algorithm. It
would be interesting to see if the convergence rate result is tight or whether
it can be improved through other techniques.

  * •

A mentioned at the end of the last section, our TD learning result does not
provide guidelines on the choice of $\delta$. It would be interesting to
investigate whether a different analysis can shed light on the optimal choice
of $\delta.$

#### Acknowledgment:

  * •

I thank Prof. Siva Theja Maguluri, Georgia Tech, for many stimulating
discussions and drawing my attention to several references during the course
of this work. I also thank Prof. Sean Meyn, University of Florida, for
pointing out the bias issue with fixed stepsize TD learning and why it cannot
be eliminated by averaging.

  * •

Research supported by AFOSR Grant FA9550-24-1-0002, ONR Grant
N00014-19-1-2566, and NSF Grants CNS 23-12714, CNS 21-06801, CCF 19-34986, and
CCF 22-07547.

## References

  * [1] R. S. Sutton and A. G. Barto, _Reinforcement learning: An introduction_. MIT press, 2018.
  * [2] J. Tsitsiklis and B. Van Roy, “Analysis of temporal-difference learning with function approximation,” _Advances in neural information processing systems_ , vol. 9, 1996.
  * [3] J. Bhandari, D. Russo, and R. Singal, “A finite time analysis of temporal difference learning with linear function approximation,” in _Conference on learning theory_. PMLR, 2018, pp. 1691–1692.
  * [4] R. Srikant and L. Ying, “Finite-time error bounds for linear stochastic approximation and TD learning,” in _Conference on Learning Theory_. PMLR, 2019, pp. 2803–2830.
  * [5] Z. Chen, S. T. Maguluri, S. Shakkottai, and K. Shanmugam, “A Lyapunov theory for finite-sample guarantees of Markovian stochastic approximation,” _Operations Research_ , 2023.
  * [6] H. Gupta, R. Srikant, and L. Ying, “Finite-time performance bounds and adaptive learning rate selection for two time-scale reinforcement learning,” _Advances in Neural Information Processing Systems_ , vol. 32, 2019.
  * [7] T. T. Doan, “Finite-time analysis and restarting scheme for linear two-time-scale stochastic approximation,” _SIAM Journal on Control and Optimization_ , vol. 59, no. 4, pp. 2798–2819, 2021.
  * [8] ——, “Nonlinear two-time-scale stochastic approximation convergence and finite-time performance,” _IEEE Transactions on Automatic Control_ , 2022.
  * [9] D. Ruppert, “Efficient estimations from a slowly convergent Robbins-Monro process,” Cornell University Operations Research and Industrial Engineering, Tech. Rep., 1988.
  * [10] B. T. Polyak, “A new method of stochastic approximation type,” _Avtomatika i telemekhanika_ , pp. 98–107, 1990.
  * [11] B. T. Polyak and A. B. Juditsky, “Acceleration of stochastic approximation by averaging,” _SIAM journal on control and optimization_ , vol. 30, no. 4, pp. 838–855, 1992.
  * [12] E. Moulines and F. Bach, “Non-asymptotic analysis of stochastic approximation algorithms for machine learning,” _Advances in neural information processing systems_ , vol. 24, 2011.
  * [13] V. Borkar, S. Chen, A. Devraj, I. Kontoyiannis, and S. Meyn, “The ODE method for asymptotic statistics in stochastic approximation and reinforcement learning,” _arXiv preprint arXiv:2110.14427_ , 2021.
  * [14] W. Mou, A. Pananjady, M. J. Wainwright, and P. L. Bartlett, “Optimal and instance-dependent guarantees for Markovian linear stochastic approximation,” _arXiv preprint arXiv:2112.12770_ , 2021.
  * [15] A. Durmus, E. Moulines, A. Naumov, and S. Samsonov, “Finite-time high-probability bounds for polyak-ruppert averaged iterates of linear stochastic approximation,” _arXiv preprint arXiv:2207.04475_ , 2022.
  * [16] S. U. Haque, S. Khodadadian, and S. T. Maguluri, “Tight finite time bounds of two-time-scale linear stochastic approximation with Markovian noise,” _arXiv preprint arXiv:2401.00364_ , 2023.
  * [17] J. Hu, V. Doshi, and D. Y. Eun, “Central limit theorem for two-timescale stochastic approximation with Markovian noise: Theory and applications,” _arXiv preprint arXiv:2401.09339_ , 2024.
  * [18] X. Li, W. Yang, J. Liang, Z. Zhang, and M. I. Jordan, “A statistical analysis of Polyak-Ruppert averaged Q-learning,” in _International Conference on Artificial Intelligence and Statistics_. PMLR, 2023, pp. 2207–2261.
  * [19] A. Röllin, “On quantitative bounds in the mean martingale central limit theorem,” _Statistics & Probability Letters_, vol. 138, pp. 171–176, 2018.
  * [20] A. Anastasiou, K. Balasubramanian, and M. A. Erdogdu, “Normal approximation for stochastic gradient descent via non-asymptotic rates of martingale CLT,” in _Conference on Learning Theory_. PMLR, 2019, pp. 115–137.
  * [21] T. Gallouët, G. Mijoule, and Y. Swan, “Regularity of solutions of the Stein equation and rates in the multivariate central limit theorem,” _arXiv preprint arXiv:1805.01720_ , 2018.
  * [22] X. Fang, Q.-M. Shao, and L. Xu, “Multivariate approximations in Wasserstein distance by Stein’s method and Bismut’s formula,” _Probability Theory and Related Fields_ , vol. 174, pp. 945–979, 2019.
  * [23] M. Metivier and P. Priouret, “Applications of a Kushner and Clark lemma to general classes of stochastic algorithms,” _IEEE Transactions on Information Theory_ , vol. 30, no. 2, pp. 140–151, 1984.
  * [24] A. Benveniste, M. Métivier, and P. Priouret, _Adaptive algorithms and stochastic approximations_. Springer Science & Business Media, 2012, vol. 22.
  * [25] R. Douc, E. Moulines, P. Priouret, and P. Soulier, _Markov Chains_. Springer, 2018.
  * [26] E. Bolthausen, “Exact convergence rates in some martingale central limit theorems,” _The Annals of Probability_ , pp. 672–688, 1982.
  * [27] C. K. Lauand and S. Meyn, “The curse of memory in stochastic approximation,” in _Proceedings of the IEEE Conference on Decision Control_. IEEE, 2023.
  * [28] D. Huo, Y. Chen, and Q. Xie, “Bias and extrapolation in Markovian linear stochastic approximation with constant stepsizes,” in _Abstract Proceedings of the 2023 ACM SIGMETRICS International Conference on Measurement and Modeling of Computer Systems_ , 2023, pp. 81–82.
  * [29] D. Nagaraj, X. Wu, G. Bresler, P. Jain, and P. Netrapalli, “Least squares regression with Markovian data: Fundamental limits and algorithms,” _Advances in neural information processing systems_ , vol. 33, pp. 16 666–16 676, 2020.
  * [30] A. Roy and K. Balasubramanian, “Online covariance estimation for stochastic gradient descent under Markovian sampling,” _arXiv preprint arXiv:2308.01481_ , 2023.
  * [31] A. Dieuleveut, A. Durmus, and F. Bach, “Bridging the gap between constant step size stochastic gradient descent and Markov chains,” _Annals of Statistics_ , 2020.
  * [32] C. Stein, “A bound for the error in the normal approximation to the distribution of a sum of dependent random variables,” in _Proceedings of the Sixth Berkeley Symposium on Mathematical Statistics and Probability, Volume 2: Probability Theory_ , vol. 6. University of California Press, 1972, pp. 583–603.
  * [33] N. Ross, “Fundamentals of Stein’s method,” _Probability Surveys_ , 2011.
  * [34] L. Arnold, “Stochastic differential equations,” _New York_ , vol. 2, 1974.
  * [35] A. D. Barbour, “Stein’s method for diffusion approximations,” _Probability theory and related fields_ , vol. 84, no. 3, pp. 297–322, 1990.
  * [36] F. Gotze, “On the rate of convergence in the multivariate CLT,” _The Annals of Probability_ , pp. 724–739, 1991.
  * [37] L. H. Chen, L. Goldstein, and Q.-M. Shao, _Normal approximation by Stein’s method_. Springer Science & Business Media, 2010.
  * [38] S. P. Meyn and R. L. Tweedie, _Markov chains and stochastic stability_. Springer Science & Business Media, 2012.
  * [39] A. M. Makowski and A. Shwartz, “The Poisson equation for countable Markov chains: probabilistic methods and interpretations,” in _Handbook of Markov Decision Processes: Methods and Applications_. Springer, 2002, pp. 269–303.
  * [40] P. W. Glynn and S. P. Meyn, “A Liapounov bound for solutions of the Poisson equation,” _The Annals of Probability_ , pp. 916–931, 1996.
  * [41] P. W. Glynn and A. Infanger, “Solution representations for Poisson’s equation, martingale structure, and the Markov chain central limit theorem,” _Stochastic Systems_ , 2023.
  * [42] S. M. Ross, _Applied probability models with optimization applications_. Courier Corporation, 2013.
  * [43] S. Meyn, _Control techniques for complex networks_. Cambridge University Press, 2008.
  * [44] D. Bertsekas and J. N. Tsitsiklis, _Neuro-dynamic programming_. Athena Scientific, 1996.
  * [45] V. S. Borkar, _Stochastic approximation: a dynamical systems viewpoint_. Springer, 2009, vol. 48.
  * [46] H. J. Kushner and G. Yin, _Stochastic approximation and recursive algorithms and applications.(2003)_. Springer, 2003.
  * [47] I. Kontoyiannis and S. P. Meyn, “Spectral theory and limit theorems for geometrically ergodic markov processes,” _The Annals of Applied Probability_ , vol. 13, no. 1, pp. 304–362, 2003.
  * [48] T. Bonis, “Stein’s method for normal approximation in wasserstein distances with application to the multivariate central limit theorem,” _Probability Theory and Related Fields_ , vol. 178, no. 3-4, pp. 827–860, 2020.

##  Appendix A Properties of $E(||\Delta_{k}||^{2}||$

The notation that we will use for the constants are just for this section of
the appendix alone and so we may occasionally reuse some letters that were
used for other constants elsewhere in the paper for constants here. We start
with the following result from [4, Theorem 12]: there exists a positive
integer $k^{*}$ such that for any $k\geq\hat{k},$ we have

$\displaystyle
E\left[\|{\Delta}_{k}\|^{2}\right]\leq\kappa_{1}\left(\prod_{j=\hat{k}}^{k-1}a_{j}\right)+\kappa_{2}\sum_{j=\hat{k}}^{k-1}b_{j}\left(\prod_{l=j+1}^{k-1}a_{l}\right),$

where $a_{j}=1-{\kappa_{3}\epsilon_{j}},$
$b_{j}=\kappa_{4}\epsilon_{j}^{2}\log(1/\epsilon_{j}),$ and
$\kappa_{1},\kappa_{2},\kappa_{3},\kappa_{4}$ are positive constants. Since
$a_{j}\leq\exp(-\kappa_{3}\epsilon_{j}),$ the first term goes to zero
exponentially fast. The second term is upper bounded by

$\displaystyle\kappa_{2}\kappa_{4}\sum_{j=0}^{k-1}\epsilon_{j}^{2}\log\left(\frac{1}{\epsilon_{j}}\right)\exp\left(-\kappa_{3}\sum_{l=j+1}^{k-1}\epsilon_{l}\right)$
$\displaystyle\leq\kappa_{2}\kappa_{4}\left(\sum_{j=0}^{k-1}\epsilon_{j}^{4}\log^{2}\left(\frac{1}{\epsilon_{j}}\right)\right)^{1/2}\left(\exp\left(-2\kappa_{3}\sum_{l=j+1}^{k-1}\epsilon_{l}\right)\right)^{1/2}=\mathcal{O}\left(\frac{\log
k}{k^{2\delta-0.5}}\right)$

##  Appendix B Properties of $\Upsilon_{j}$

In this section of the appendix, we characterize the behavior of
$\Upsilon_{j}$ to understand the asymptotics as $\delta\rightarrow 1.$ The
notation that we will use for the constants are just for this section of the
appendix alone and so we may occasionally reuse some letters that were used
for other constants elsewhere in the paper for constants here. Our analysis
here is intended to show that, while Polyak-Ruppert averaged TD learning
converges to a Gaussian random vector for every $\delta\in(0.5,1)$ at rate
$\mathcal{O}(1/n^{\delta-0.5})$, the constant hidden in the $\mathcal{O}$
notation may go to $\infty$ as $\delta\rightarrow 1.$ For this purpose, we
examine certain expressions which are shown be either bounded or converging to
zero in [11], and derive more explicit bounds for them as a function of
$\delta$.

From [11, page 846], we know that $\Upsilon_{j}^{t}$ can be written as

$\Upsilon_{j}^{t}=C_{j}^{t}+\bar{A}^{-1}D_{j}^{t},$

where

$||C_{j}^{t}||_{op}\leq\frac{\epsilon_{j}-\epsilon_{j+1}}{\epsilon_{j}}\sum_{i=j}^{t}m_{j}^{i}e^{-\lambda
m_{j}^{i}},\qquad D_{j}^{t}=\prod_{l=j}^{t-1}(I-\epsilon_{l}\bar{A}),\mbox{
and }||D_{j}^{t}||_{op}\leq K\exp(-\alpha m_{j}^{t-1}),$

for $t\geq j+1$ and some $K,\alpha,\lambda>0,$ and
$m_{j}^{i}=\sum_{k=j}^{i}\epsilon_{k}.$ Since

$\epsilon_{k}=1/(k+1)^{\delta}\geq\mu/(k+1)$

for any $\mu>0$ and any $k\geq\lceil\mu^{1/(1-\delta)}=:K_{\mu,\delta}\rceil,$

$m_{j}^{i}=\sum_{k=j}^{i}\epsilon_{k}\geq\mu\sum_{k=j}^{i}\frac{1}{k+1}\geq\mu\log\frac{i+1}{j+1}$

for all $j\geq K_{\mu,\delta}.$ Thus, since $t\geq j+1,$

$||D_{j}^{t}||_{op}\leq K\exp(-\alpha\mu\log(t)/{j+1}))\leq
K(\frac{j+1}{t})^{\alpha\mu}\leq K,\qquad j\geq K_{\mu,\delta}.$

For $j<K_{\mu,\delta},$

$||D_{j}^{t}||_{op}\leq K\exp(-\alpha
m_{j}^{K_{\mu,\delta}-1})+K(\frac{K_{\mu,\delta}+1}{t})^{\alpha\mu}.$

Next, note that by choosing $\mu$ sufficiently large such that $\alpha\mu\geq
2$, $\sum_{j=0}^{t-1}||D_{j}^{t}||_{op}$ is bounded and hence

$\frac{1}{t}\sum_{j=0}^{t-1}||D_{j}^{t}||_{op}=\mathcal{O}\left(\frac{1}{t}\right),$

Note that since $K_{\mu,\delta}$ goes to infinity as $\delta\rightarrow 1,$
the rate at which the average of $||D_{j}^{t}||$ goes to zero has a constant
that blows up to infinity as $\delta\rightarrow 1.$

Next, let us consider $C_{j}^{i}.$ First we note that, using the fact that
$(1+x)^{\delta}\leq 1+\delta x$ for $x\geq 0$ and $\delta\in(0,1),$ we have

$\frac{\epsilon_{j}-\epsilon_{j+1}}{\epsilon_{j}}\leq\frac{\delta}{j+2}.$

Thus,

$\displaystyle||C_{j}^{t}||_{op}\leq\frac{\epsilon_{j}-\epsilon_{j+1}}{\epsilon_{j}}\sum_{i=j}^{t}m_{j}^{i}e^{-\lambda
m_{j}^{i}}\leq\frac{1}{j+2}\sum_{i=j}^{t}m_{j}^{i}e^{-\lambda m_{j}^{i}}$
$\displaystyle\qquad\leq\frac{1}{(j+2)\lambda^{p}}\sum_{i=j}^{t}\frac{1}{(m_{j}^{i})^{p-1}}\leq\frac{1}{(j+2)\lambda^{p}}\sum_{i=1}^{t}\frac{1}{(m_{1}^{i})^{p-1}},$

for any $p>1.$ Noting that $m_{1}^{i}$ is $O(i^{1-\delta})$ and choosing $p$
such that $(p-1)(1-\delta)\geq 2$ shows that $||C_{j}^{t}||_{op}$ is bounded
independent of $\delta$ and

$\frac{1}{t}\sum_{j=0}^{t-1}||C_{j}^{t}||_{op}=\mathcal{O}\left(\frac{\log
t}{t}\right),$

with a constant independent of $\delta.$

Finally, since $||\Upsilon_{j}||_{op}^{2}\leq
2(||C_{j}^{t}||_{op}^{2}+||D_{j}^{t}||_{op}^{2}),$ and $D_{j}^{t}$ and
$C_{j}^{t}$ are uniformly bounded for each fixed $\delta,$ it follows that

$\frac{1}{t}\sum_{j=0}^{t-1}||\Upsilon_{j}^{t}||_{op}^{2}=\mathcal{O}\left(\frac{\log
t}{t}\right),$

with a constant that goes to $\infty$ as $\delta\rightarrow 1.$