Some Limit Properties of Markov Chains Induced by Stochastic Recursive Algorithms

TL;DR

This paper studies the limit properties of Markov chains generated by stochastic recursive algorithms, showing convergence to deterministic contraction trajectories and invariant distributions, with applications to machine learning and dynamic programming.

Contribution

It establishes weak convergence and ergodic properties of Markov chains induced by iterated random operators, extending understanding of stochastic recursive algorithms.

Findings

Random sequences converge weakly to contraction operator trajectories.

Time averages of the sequences converge to invariant distribution means.

Applications include logistic regression and dynamic programming algorithms.

Abstract

Recursive stochastic algorithms have gained significant attention in the recent past due to data driven applications. Examples include stochastic gradient descent for solving large-scale optimization problems and empirical dynamic programming algorithms for solving Markov decision problems. These recursive stochastic algorithms approximate certain contraction operators and can be viewed within the framework of iterated random operators. Accordingly, we consider iterated random operators over a Polish space that simulate iterated contraction operator over that Polish space. Assume that the iterated random operators are indexed by certain batch sizes such that as batch sizes grow to infinity, each realization of the random operator converges (in some sense) to the contraction operator it is simulating. We show that starting from the same initial condition, the distribution of the random sequence generated by the iterated random operators converges weakly to the trajectory generated by the contraction operator. We further show that under certain conditions, the time average of the random sequence converges to the spatial mean of the invariant distribution. We then apply these results to logistic regression, empirical value iteration, and empirical Q value iteration for finite state finite action MDPs to illustrate the general theory develop here.

Figures

Equations

$$\underset{x\in {\mathbb{R}}^n}{\min }L(x)=\frac{1}{N}\sum _{i=1}^N{L}_i(x).$$

$$y_{k+1}=y_k-\beta {\nabla }_{x}L(y_k)=:T(y_k),$$

$${\hat{z}}_{k+1}^n={\hat{z}}_k^n-\frac{\beta }{n}\sum _{i\in {\mathcal{N}}_k}{\nabla }_x{L}_i({\hat{z}}_k^n)=:{\hat{T}}_k^n({\hat{z}}_k^n).$$

$$\displaystyle T(v)(s,a)=\min_{a}\bigg{[}c(s,a)+\alpha\mathds{E}\left[v(s^{\prime})|s,a\right]\bigg{]}.$$

$${\hat{T}}_k^n(v)(s,a)=\underset{a}{\min }\left[c(s,a)+\alpha \frac{1}{n}\sum _{i=1}^{n}v(s_i^{\prime }(s,a))\right].$$

$$\Vert {\hat{z}}_k^n-x^{\ast }\Vert \le \Vert {\hat{z}}_k^n-y_k\Vert +\Vert y_k-x^{\ast }\Vert \le \Vert {\hat{z}}_k^n-y_k\Vert +{\alpha }^k\Vert x_0-x^{\ast }\Vert .$$

$$\rho (T(x_1),T(x_2))\le \alpha \rho (x_1,x_2)\text{ and }T(x^{\ast })=x^{\ast }.$$

$$y_k=T(y_{k-1})\text{ for }k\in \mathbb{N}.$$

$$\rho (y_k,x^{\ast })\le {\alpha }^k\rho (y_0,x^{\ast }).$$

$$Q_n(B|x)=\text{mathds}P\left\{{\hat{z}}_k^n\in B|{\hat{z}}_{k-1}^n=x\right\}$$

$${\mu }_n(B_0\times B_1\times B_2\times \dots )={\int }_{B_0\times B_1\times B_2\times \dots }\text{mathds}{1}_{\{y_0\}}(d{x}_0)Q_n(d{x}_1|x_0)Q_n(d{x}_2|x_1)\dots ,$$

$$\psi (B_0\times B_1\times B_2\times \dots )={\int }_{B_0\times B_1\times B_2\times \dots }\text{mathds}{1}_{\{y_0\}}(d{x}_0)\text{mathds}{1}_{\{y_1\}}(d{x}_1)\text{mathds}{1}_{\{y_2\}}(d{x}_2)\dots .$$

$${\pi }^n(B)={\int }_{\mathcal{X}}{Q}_n(B|x){\pi }^n(dx).$$

$$f(u;x)=\sigma (u^{\mathsf{T}}x),\text{ where }\sigma (t)=\frac{1}{1+e^{-t}}.$$

$$L(x)=\frac{1}{N}\sum _{i=1}^N{L}_i(x),\text{ where }{L}_i(x):=-l^i\log f(u^i;x)-(1-l^i)\log (1-f(u^i;x)),$$

$$\displaystyle\nabla_{x}L_{i}(x)=(f(u^{i};x)-l^{i})u^{i},\quad\nabla^{2}_{xx}L_{i}(x)=\Big{(}f(u^{i};x)(1-f(u^{i};x))\Big{)}u^{i}u^{i\mathsf{T}}.$$

$$x_{k+1}=T(x_k):=x_k-\beta {\nabla }_{x}L(x_k)=x_k-\frac{\beta }{N}\sum _{i=1}^N(f(u^i;x_k)-l^i)u^i,$$

$${\hat{z}}_{k+1}^n={\hat{T}}_k^n({\hat{z}}_k^n)={\hat{z}}_k^n-\beta \frac{1}{n}\sum _{j=1}^n{\nabla }_x(L_j({\hat{z}}_k^n)).$$

$$L_i(x):=-l^i\log f(u^i;x)-(1-l^i)\log (1-f(u^i;x))+\frac{\lambda }{2}\Vert x{\Vert }_2^2\text{ with }\lambda >0,$$

$$\underset{n\to \infty }{\lim }\text{mathds}P\left\{\Vert {\hat{T}}_k^n(x)-T(x){\Vert }_2>ϵ\right\}=0.$$

$$\underset{n\to \infty }{\lim }\underset{x\in \mathcal{K}}{\sup }\text{mathds}P\left\{\Vert {\hat{T}}_k^n(x)-T(x){\Vert }_2>ϵ\right\}=0.$$

$$v^{\gamma }(s):=\text{mathds}E\left[\sum _{k=0}^{\infty }{\alpha }^{k}c(s_k,a_k)|s_0=s,a_k=\gamma (s_k)\right],$$

$$T(v)(s)=\inf_{a\in A}\Big{\{}c(s,a)+\alpha\sum_{s^{\prime}\in\mathcal{S}}v(s^{\prime})p(s^{\prime}|s,a)\Big{\}}.$$

$$\text{Initialize v0 randomly and iterate }{v}_{k+1}(s)=T(v_k)(s).$$

$${\hat{T}}_k^n({\hat{v}}_k^n)(s)=\underset{a\in \mathcal{A}}{\inf }\left\{c(s,a)+\alpha \frac{1}{n}\sum _{i=1}^n{\hat{v}}_k^n(s_{k,i}^{\prime }(s,a))\right\},$$

$$\text{Initialize v0 randomly and let }{\hat{v}}_{k+1}^n(s)={\hat{T}}_k^n({\hat{v}}_k^n)(s).$$

$$\underset{n\to \infty }{\lim }\text{mathds}P\left\{\Vert {\hat{T}}_k^n(v)-T(v){\Vert }_{\infty }>ϵ\right\}=0.$$

$$\underset{v\in \mathcal{K}}{\sup }\text{mathds}P\left\{\Vert {\hat{T}}_k^n(v)-T(v){\Vert }_{\infty }>ϵ\right\}\le 2|\mathcal{S}||\mathcal{A}|\exp \left(-\frac{ϵn}{|\mathcal{S}|k^2}\right),$$

$$\underset{n\to \infty }{\lim }\underset{v\in \mathcal{K}}{\sup }\text{mathds}P\left\{\Vert {\hat{T}}_k^n(v)-T(v){\Vert }_{\infty }>ϵ\right\}=0.$$

$$\underset{x\in K}{\sup }\text{mathds}P\left\{\rho ({\hat{T}}_k^n(x),T(x))>ϵ\right\}<\delta \text{ for all }n\ge N_{ϵ,\delta }(K).$$

Full text

{AMS}

93E35, 60J20, 68Q32

1 Introduction

There has been a surge of interest in using randomization to reduce computational burden in machine learning and reinforcement learning. For instance, in training neural networks with a large amount of data, stochastic gradient descent is frequently employed instead of the usual gradient descent. In data-driven Markov decision problems, empirical dynamic programming and generative models have been employed to determine approximately optimal policies and value functions. In these algorithms, instead of computing the expected value of certain functions at each step of the iteration, one computes the empirical expected value that is rather easy to compute if enough data is available. This simple trick reduces the runtime to determine a reasonably good solution.

It turns out that the outputs of these recursive stochastic algorithms (RSAs) can be viewed as Markov chains. Indeed, if the parameters of the algorithm do not change with iteration, then the RSAs can be thought of as an iterated random operator acting onto certain spaces. Consider, for instance, the case of stochastic gradient descent, where the stepsize remains constant, data samples picked at every iteration are i.i.d., and the number of data samples remain constant. Then, each step of the stochastic gradient descent algorithm can be viewed as a random operator. To see this, let us consider the problem of minimizing a sum of $N$ functions, $L_{i}:\mathbb{R}^{n}\rightarrow\mathbb{R}$, $i=1,\ldots,N$:

[TABLE]

In usual gradient descent, one fixes a stepsize $\beta>0$, and runs the iteration

[TABLE]

where we used $T:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ to denote the exact gradient descent map. We note here that under reasonable assumptions on $L$ and $\beta$, $T$ becomes a contraction operator under some norm on $\mathbb{R}^{n}$ (usually $\ell_{2}$ norm is used).

In contrast, in stochastic gradient descent, the operator applied at every step of the algorithm changes. At time step $k$ of the stochastic gradient descent algorithm, let $\mathcal{N}_{k}:=\{i_{1},\ldots,i_{n}\}$ be the set of $n$ indices that are sampled independently and uniformly from the set of all indices $\{1,\ldots,N\}$. Then, we have

[TABLE]

Since the set of (random) indices $\mathcal{N}_{k}$ is i.i.d., the operator $\hat{T}^{n}_{k}$ is independent of the past maps and is “identically distributed”. This implies that the (random) sequence $(\hat{z}^{n}_{k})_{k\in\mathbb{N}}$ is a Markov chain. It should also be noted that the exact gradient descent operator $T$ and stochastic gradient descent operator $\hat{T}^{n}_{k}$ are related – $\hat{T}^{n}_{k}(x)$ is a consistent estimate of the $T(x)$ for any $x\in\mathbb{R}^{n}$.

A similar setup is considered in empirical dynamic programming using a generative model for dynamic decision process. Consider a controlled Markov process in which $s$ is the state of a system and $a$ is the action of the decision maker. Let $p(s^{\prime}|s,a)$ denote the transition probability of the next state being $s^{\prime}$ given the current state $s$ and action $a$ and $c(s,a)$ be the corresponding cost. We use $\alpha$ to denote the discount parameter. In the value iteration algorithm, one needs to compute $\mathds{E}\left[v(s^{\prime})|s,a\right]$, where $v$ is some real-valued function of the state. This leads to the usual Bellman operator $T$ that acts on the space of value functions and is given by

[TABLE]

It is not difficult to see that $T$ is a contraction operator with the contraction coefficient $\alpha$, when the space of value functions is endowed with the sup norm. If there is enough data, one can replace $\mathds{E}\left[v(s^{\prime})|s,a\right]$ with its “empirical” average, given by $\frac{1}{n}\sum_{i=1}^{n}v(s^{\prime}_{i}(s,a))$, where $\{s^{\prime}_{i}(s,a)\}_{i=1}^{n}$ are $n$ samples of the next state given that the current state-action pair is $(s,a)$. Thus, the random Bellman operator $\hat{T}^{n}_{k}$ acts on the space of value functions and is given by

[TABLE]

In this case as well, $\hat{T}^{n}_{k}(v)$ is a consistent estimate of the $T(v)$ for any value function $v$. When the state or action space is continuous, then the above operator also features an additional function fitting component. The analysis of such algorithms involves understanding the error introduced due to function fitting, as well as the number of samples used at every iteration.

In fact, in both of the examples above, it is readily observed that the original operator $T$ is a contraction map. This is no accident – contraction mapping theorem forms the bedrock of most of the proofs of convergence algorithms used in optimization or MDP problems. Some examples are noted below:

The Bellman operator for a class of stochastic shortest path problems is a contraction operator under an appropriate weighted sup norm [9, 7]. 2. 2.

The Bellman operator for a class of average cost MDPs is a contraction operator under the span seminorm over the quotient space [63, 2]. 3. 3.

Some variational inequality problems involve contraction operators under the usual 2 norm [29, p. 143]. 4. 4.

The Bellman operator for continuous state MDP is a contraction operator over the space of continuous and bounded functions over the state space of the MDP (it requires a variety of assumptions as elucidated in [40, 41, 28]). 5. 5.

The resolvent of a strongly monotone operator is a contraction operator [65]. Many other contraction operators used in the context of optimization algorithms are discussed in Section 5 of [65].

In data driven applications, computing the exact (contraction) operator is either computationally expensive or impossible. Thus, one has to use random mappings, that are drawn independent of the past operators, and that simulate the effect of contraction mapping. The two examples explained above are merely instances of this more general methodology. The primary purpose of the paper is to devise sufficient conditions on the relationship between random maps $\hat{T}^{n}_{k}$ and the deterministic contraction operator $T$ so that (a) the convergence properties of new algorithms can be readily established, and (b) find common threads between the convergence.

1.1 Our Contribution

The primary contribution of this paper is to conceptually unify the convergence analysis of certain RSAs in optimization, machine learning, and reinforcement learning using the tools from the Markov chain theory. This is achieved by leveraging several results available for convergence and stability for Feller Markov chains established in [19], [56, Sec. 18.5], [17, Sec. 8], and [49]. Our key contributions are summarized below.

Consider the deterministic sequence $y_{k}=T\circ\ldots\circ T(y_{0})$ ($k$ compositions of the exact operator $T$) and the Markov chain $\hat{z}^{n}_{k}=\hat{T}^{n}_{k-1}\circ\ldots\circ\hat{T}^{n}_{0}(\hat{z}^{n}_{0})$ with $\hat{z}^{n}_{0}=y_{0}$. It is natural to assume that using the deterministic operator yields the best convergence guarantee. We use independent samples of data to approximate $T$ by $\hat{T}^{n}_{k}$, then it introduces error at every time step. How does this error accumulate? Most authors bound $\|\hat{z}^{n}_{k}-x^{*}\|$, where $x^{*}$ is the fixed point of the operator $T$. Instead, we are interested in the error $\|\hat{z}^{n}_{k}-y_{k}\|$, noting that triangle inequality yields

[TABLE]

We derive a sufficient condition on the random operators $\hat{T}^{n}_{k}$ and its relationship with the exact operator $T$ so that the sequence of distributions over the sequence $(\hat{z}^{n}_{k})_{k\in\mathbb{N}}$ output by the RSA converges in weak* topology to the unit mass over the trajectory $(y_{k})_{k\in\mathbb{N}}$ output by the exact algorithm as the parameter $n\rightarrow\infty$. This further implies that $\hat{z}^{n}_{k}$ is close to $y_{k}$ with high probability. Using the inequality in (1) above, we conclude that the error term $\|\hat{z}^{n}_{k}-x^{*}\|$ is less than $\alpha^{k}\|x_{0}-x^{*}\|$ plus a small error with high probability.

We further show that the sufficient condition is satisfied in a sufficiently general class of problems encountered in stochastic gradient descent for strongly convex loss functions and synchronous empirical dynamic programming for MDPs with discounted cost criterion. These examples serve to illustrate how to apply the results to derive this property of an RSA. 2. 2.

Existence of invariant distributions of RSAs is an important property, as it implies some form of stability of the algorithm. A Markov chain does not admit invariant distribution if it features a cyclic behavior or diverges to infinity (there could be other reasons for the non-existence of an invariant distribution, but these two are more common). In an RSA, we do not usually expect a cyclic behavior due to randomization. Thus, if an RSA admits an invariant distribution, then it implies that the iterates will not diverge with high probability. Thus, establishing the existence of an invariant distribution is an important problem, which we address here.

We show that the Markov chains generated by many RSAs satisfy the weak Feller property, that is, if $f$ is a continuous and bounded function, then $\hat{z}^{n}_{0}\mapsto\mathds{E}\left[f(\hat{z}^{n}_{1})|\hat{z}^{n}_{0}\right]$ is also a continuous and bounded function. The existence (and in some case, the uniqueness) of an invariant measure of Feller Markov chains has been presented in [17]. We apply these results to conclude that under some reasonable assumptions, the chains generated by stochastic gradient descent and empirical dynamic programming algorithms admit invariant distributions. In certain cases, we can show that this invariant distribution is unique. 3. 3.

There has been a sustained interest in using time-averages in stochastic gradient descent and deep Q learning. Particularly, references [67, 55, 45, 54, 78] propose that fixing the stepsize in stochastic gradient descent algorithms and using the average of the tail of the random sequence leads to a better performance of the trained algorithm. Within the context of reinforcement learning, [3] and [82] propose averaging the deep Q function iterates to arrive at a solution with lower variance.

Indeed, we show here that under some conditions on the random operators, the variance reduction property of time-average (or in these cases, tail average) is largely due to the fact that the Markov chain output by RSA may be admitting a unique invariant distribution. This part leverages the law of large numbers for Markov chains, presented in [19, 56]. 4. 4.

We complement the theoretical contributions with numerical simulations of two RSAs – stochastic gradient descent for logistic regression, empirical value iteration for discounted MDP with a generative model, and synchronous batch Q value iteration for discounted MDP.

While we present complete proofs of two of our main results stated here, we admit that our proofs require minor tweaks of existing results in the literature. The need for presenting the complete proofs are twofold: Our hypotheses differ in some ways from the hypotheses presented in the standard texts, particularly in [49], [56, Section 18.5], and [19]. Moreover, to construct the complete proofs using these texts under our hypotheses require substantial effort on the part of the reader. To ease this burden, we chose to furnish the complete proofs using the notation adopted here.

1.2 Previous Work

Convergence proofs of randomized optimization and learning algorithms are usually obtained from specifically tailored arguments, which are not generalizable to other settings. For instance, the convergence of stochastic gradient descent, stochastic variance reduction gradient descent (SVRG), and stochastic average gradient (SAG) descent follow a completely different, and often involved, sequence of arguments [18, 14, 47, 36, 67, 24]. The argument usually starts with identifying some conditions on the functions, such that for every iteration $k$, one can upper bound $L(\hat{z}^{n}_{k})-L(x^{*})$ (where we used the notation introduced above) by a function that decays as $k$ grows. These tailored methods usually also yield the convergence rates specific to those algorithms.

It would be conceptually elegant to determine a set of more general conditions which can be readily applied to these algorithms and many of its variants to establish the asymptotic convergence to the fixed point of the map. The stochastic approximation theory is one such elegant theory [16, 43, 50, 15]. There are two types of stochastic approximation algorithms – one with decreasing (also called tapering or diminishing) stepsize and other one with constant stepsize. In decreasing stepsize algorithms, the stepsize has to converge to 0 as the number of iteration goes to infinity (the stepsize is not summable, but is square summable). This leads to the almost sure convergence guarantee to the fixed point in the limit. For constant stepsize, the theory says that the iterates will eventually enter a neighborhood of the fixed point and do a random walk within that neighborhood.

We now present a sample of decreasing stepsize RSAs whose convergence is ascertained using the stochastic approximation theory. Under reasonable assumptions on the loss function, stochastic gradient descent and distributed asynchronous gradient descent methods converge almost surely to the optimal solution [79, 31, 32]. The convergence of reinforcement learning algorithms usually invoke some version of the stochastic approximation theorem. Reference [44, 80, 1, 10] studies the convergence of various types of Q learning algorithms developed for discounted cost or average cost MDPs with finite state and action spaces. The convergence of on-policy reinforcement learning algorithm SARSA is established in [74]. More recently, the stochastic approximation theory has been used to establish the convergence of policy gradient, temporal difference, and other related methods in [76, 90, 94, 91, 88, 89]. For more information on various reinforcement learning algorithms and their convergence proofs, we refer the reader to books [2, 77, 9] and recent survey papers [92, 93].

Decreasing stepsize RSAs do not yield approximate solutions in a reasonable time frame. As a result, constant stepsize algorithms are gaining traction as a way to speed up the computation at the cost of tolerating a small error in the final result; see, for example, [30, 58, 4, 27, 26], where constant stepsizes are used in the context of the stochastic gradient descent-type algorithms and [6, 75, 59, 68] for their usage in the reinforcement learning algorithms.

Constant stepsize stochastic approximation over finite dimensional state space has been studied in [16, 43], where the authors derive the asymptotic concentration results. It is well-understood that in stochastic gradient descent with constant stepsize, the sequence generated by the algorithm gets closer to the optimal solution, but then does a random walk around the optimal solution [26]. The closeness of the random walk to the optimal solution depends on the number of random samples one uses at each iteration of the algorithm. Similarly, in Q learning algorithm, constant stepsize Q learning has been studied in [6] (both synchronous and asynchronous version are studied). Convergence of constant stepsize temporal difference methods with linear function approximation is studied in [51, 23, 11, 35, 75]. Empirical value iteration and their variants, studied in a number of works under varying assumptions [86, 87, 21, 57, 22, 46, 48, 37, 73], are also examples of constant stepsize stochastic approximation algorithms.

When the stepsize is taken as a constant in a RSA, then the output of the RSA forms a Markov chain. The goal of this paper is to study the limit properties of such a Markov chain. There are two ways the limit can be taken. Either the sample size $n$ used at every time step can grow to infinity or the number of iterations $k$ can escape to infinity. We study both the limits in this paper. We note here that the generality of the model and minimal assumptions do not allow us to derive a finite time guarantee, which has significant importance in the machine learning and the reinforcement learning communities. Further, our proof approach is not algorithm-dependent. We leave these important problems for a future work.

Our work is largely motivated by the analysis of empirical dynamic programming in [37]. This work viewed empirical dynamic programs within the framework of iterated random operators. It used stochastic dominance based arguments to derive bounds on the asymptotic probability of error (between the random outputs of the algorithm and the optimal solution) being large. Inspired by this work, we extended the arguments to empirical relative value iteration in [34]. We further relaxed some conditions on random operators assumed in [37] in our follow up work [33]. The aim of this paper is to further expand the analysis and present conditions on random operators and its relationship to the exact operator to arrive at insightful conclusions about the random sequences generated by these RSAs.

1.3 Outline of the Paper

The paper is organized as follows. Section 2 presents a common mathematical framework to study the problem of convergence and stability of Markov chains induced by RSAs. We also state the three main questions we address. Section 3 presents some motivating examples where the mathematical framework we develop can be applied. Through these examples, we also illustrate certain desirable properties that the random operators enjoy. In Section 4, we show that the distributions over the trajectories generated by RSA converges to the Dirac mass over the trajectory generated by the exact algorithm. This constitutes the first main result of the paper. In Section 5, we study some sufficient conditions on the operators $\hat{T}^{n}_{k}$ such that the resulting Markov chain admits an invariant distribution. We also study conditions under which the invariant distribution is unique. Section 6 then introduces the assumptions and establishes the weak law of large numbers for Markov chains. This constitutes the second main result of the paper. The proofs of the two main results are presented in Sections 8 and 9. We finally conclude our discussion in Section 10.

1.4 Notations

Let $(\mathcal{A},\rho)$ be a Polish space, which is defined as a complete separable metric space with metric $\rho$. We use $\wp(\mathcal{A})$ to denote the set of all probability distributions over $\mathcal{A}$. We use $\mathds{1}_{\{a\}}\in\wp(\mathcal{A})$ to denote the Dirac mass over $a\in\mathcal{A}$. The notations $C_{b}(\mathcal{A})$ and $U_{b}(\mathcal{A})$ denote, respectively, the set of all continuous and bounded functions and uniformly continuous and bounded functions over the set $\mathcal{A}$. We use $C(\mathcal{A})$ to denote the set of (possibly unbounded) continuous functions over the set $\mathcal{A}$. We say that a sequence of measures $\{\mu_{n}\}_{n\in\mathbb{N}}\subset\wp(\mathcal{A})$ converges to $\mu$ in weak* sense iff for every $f\in C_{b}(\mathcal{A})$, $\int fd\mu_{n}\rightarrow\int fd\mu$ as $n\rightarrow\infty$. This is usually also referred to as weak convergence in probability theory literature.

We use $1_{\{\cdot\}}$ to denote an indicator function, which takes the value of 1 if $\{\cdot\}$ is true and 0 otherwise. By a slight abuse of notation, we also use $1_{F}$ to be the indicator function over a measurable set $F\subset\mathcal{A}$.

2 Problem Formulation

Let $(\mathcal{X},\rho)$ be a Polish space with metric $\rho$. Consider a contraction operator $T:\mathcal{X}\rightarrow\mathcal{X}$ with contraction coefficient $\alpha\in(0,1)$ and the unique fixed point denoted by $x^{*}$. This means

[TABLE]

Starting from any initial point $y_{0}\in\mathcal{X}$, define the iterates

[TABLE]

By the Banach contraction mapping theorem, this iteration converges to $x^{*}$. In fact, we have

[TABLE]

As discussed previously, in many instances, it is beneficial or required in many iterative algorithms to use randomization to evaluate an approximation of $T(x)$ using a random operator. We now formulate a framework to analyze the output of this RSA rigorously.

Let $(\Omega,\mathscr{F},\mathbb{P})$ be a standard probability space, where $\Omega$ is the set of uncertainties, $\mathcal{F}$ is the Borel $\sigma$-algebra over $\Omega$ and $\mathbb{P}$ be the probability distribution function over $\Omega$. Let $\hat{T}^{n}_{k}:\mathcal{X}\times\Omega\rightarrow\mathcal{X}$ be a random operator that is used at the $k^{th}$ iteration and is indexed by a natural number $n$. The index $n$ would capture, for instance, the stepsize, the number of random samples used to approximate the operator $T$, etc. Although $\hat{T}^{n}_{k}$ is a function of $\omega\in\Omega$, we will suppress this dependence for ease of exposition. Thus, $\hat{T}^{n}_{k}(x):=\hat{T}^{n}_{k}(x;\omega)$. We make the following assumption on the independence of the sequence of operators $\hat{T}^{n}_{k}$.

{assumption}

For every $x\in\mathcal{X}$, $\hat{T}^{n}_{k}(x)$ and $\hat{T}^{n}_{k^{\prime}}(x)$ are statistically independent and identically distributed for $k\neq k^{\prime}$.

2.1 Key Questions

Consider the stochastic process that starts from $\hat{z}^{n}_{0}=y_{0}$ and define $\hat{z}^{n}_{k}=\hat{T}^{n}_{k-1}(\hat{z}^{n}_{k-1})$ for all $k\in\mathbb{N}$. Due to Assumption 2 and the fact that $n$ does not change with time, the stochastic process $\{\hat{z}^{n}_{k}\}_{k\in\mathbb{N}}$ is a time-homogeneous Markov chain. One can view $\hat{z}^{n}_{k}$ as an $\mathcal{X}$-valued Markov chain with the Markov transition kernel given by

[TABLE]

for any Borel set $B\subset\mathcal{X}$. Note that $Q_{n}$ does not depend on the time index $k$, since what we have here is a time-homogeneous Markov process.

2.1.1 Convergence of Distribution of Trajectories

We are interested in deriving conditions on the random maps $\hat{T}^{n}_{k}$ under which the random sequence generated by RSA is close to the deterministic sequence generated by exact algorithm with high probability. Let us formulate the precise mathematical problem. We let $\mu_{n}\in\wp(\mathcal{X}^{\mathbb{N}})$ denote the joint distribution of the sequence $(\hat{z}^{n}_{0},\hat{z}^{n}_{1},\hat{z}^{n}_{2},\ldots)$. Endow $\mathcal{X}^{\mathbb{N}}$ with the product topology so that it becomes a Polish space. Then, $\mu_{n}$ is defined by

[TABLE]

where $B_{0},B_{1},B_{2},\ldots$ are Borel sets in $\mathcal{X}$.

In the similar vein, one can also view the iterates $(y_{k})_{k\in\mathbb{N}}$ defined in (2) as a Markov chain on the same probability space $(\Omega,\mathscr{F},\mathbb{P})$, with the distribution over this sequence defined by

[TABLE]

This is a Dirac mass on the sequence $(y_{0},y_{1},\ldots)$. Our first result, stated in Section 4, proves that under a mild assumption on the random operators $\hat{T}^{n}_{k}$, the sequence of measures $\mu_{n}$ converges in the weak* sense to $\psi$.

A similar setup was considered by Karr [49]. It studies the convergence properties of a class of Feller Markov chains parametrized by $n$ such that the transition probability $Q_{n}$ converges in some sense to a transition probability $Q$ as $n\rightarrow\infty$. Although our assumptions are slightly different, the proof essentially imitates the one in [49] except for a couple of key steps. We also discuss numerical implication of this result in Section 4.

2.1.2 Existence of Invariant Measures for Fixed $n$

For the Markov chain $(\hat{z}^{n}_{k})_{k\in\mathbb{N}}$, one of the key questions is the existence of an invariant distribution. An invariant probability distribution of the Markov chain $(\hat{z}^{n}_{k})_{k\in\mathbb{N}}$ is a probability measure $\pi^{n}$ such that for any Borel set $B\subset\mathcal{X}$, we have

[TABLE]

A desirable property of a Markov chain is to have an invariant distribution, since it implies that the RSA satisfies a form of stability. More importantly, it implies that the Markov chain will not escape to infinity with probability 1. There is a large body of literature that studies the problem of the existence of invariant measures for Harris recurrent Markov chains that take values in continuous state spaces [17, 56]. However, the Markov chain generated by the RSAs seldom satisfy the strong recurrence structure required for Harris recurrent chains.

Instead, these chains satisfy the weak Feller conditions, for which there are limited results in the literature. Nonetheless, we show that many RSAs satisfy certain desirable properties, which can be leveraged to not only guarantee the existence of an invariant distribution, but also establish the uniqueness of the invariant distribution. These properties of the random operators are discussed in Section 5. This further leads to strong conclusions about the weak law of large numbers, as we discuss next.

2.1.3 Convergence of Time Average of Iterates

The weak law of large numbers for independent and identically distributed (i.i.d.) random variables states that the time average of i.i.d. random variables converges to the mean of the distribution in probability under fairly mild conditions. In fact, such a version of the weak law of large numbers is also available for Feller Markov chains. This is established for Feller chains in [19] for chains residing in a compact Hausdorff space with a unique invariant measure, and in [56, Section 18.5] for the non-compact case under certain technical conditions, which include existence of an invariant measure. It turns out that this result can be proved simply under the uniqueness of the invariant measure if the starting point $\hat{z}^{n}_{0}$ is chosen according to certain specific distribution (in fact, we do not need other technical conditions of [56, Section 18.5]). We prove this result in Section 6, the proof of which is adapted from the results from [19] and [56].

We illustrate the theoretical results using numerical simulations for batch gradient descent, empirical value iteration, and synchronous batch Q learning in Section 7.

3 Motivating Examples

We present here two examples where we illustrate how the random operator framework can be applied.

3.1 Stochastic Gradient Descent in Logistic Regression

Logistic regression has been widely used in many binary or multi-class classification problems. For simplicity, we consider the logistic regression with binary classification. Let $\mathcal{U}\subset\mathbb{R}^{m}$ be the set of feature vectors. Let $(u^{i},l^{i})_{i=1}^{N}\subset\mathcal{U}\times\{0,1\}$ denote the labeled dataset with $N$ data points and their labels. Our task is to model conditional probability distribution of label $l$ given the feature vector $u\in\mathcal{U}$. In logistic regression, we model $\mathds{P}\left\{l=1|\ u^{i}\right\}$ as $f(u^{i};x)$ where $x\in\mathcal{X}:=\mathbb{R}^{m}$ are the parameters of $f$ to be learned from the data, where $f$ is defined below:

[TABLE]

Our goal is to compute the parameter $x$ that maximizes the log likelihood (or equivalently, minimizes the negative log likelihood) given the labeled data. The log likelihood $L:\mathcal{X}\longrightarrow\mathbb{R}$ of i.i.d data under conditional distribution $f$ is given by

[TABLE]

It can be shown that the derivatives of $L_{i}$ are given by

[TABLE]

Consequently, each $L_{i}$ is a convex function, and thus, $L$ is a convex function over the space $\mathcal{X}$. If the matrix $[u^{1}|\ldots|u^{N}]$ is full rank and $N>m$, then it immediately follows that $\nabla^{2}_{xx}L(x)$ is a full rank matrix with positive eigenvalues. Consequently, $L$ is strongly convex, and therefore has a unique minimum $x^{*}$. This minimum can be computed using the usual gradient descent algorithm. The algorithm starts at $x_{0}$, picked arbitrarily, and proceeds in the direction of $-\nabla_{x}L(x_{k})$ in small steps of size $\beta$:

[TABLE]

where $T:\mathcal{X}\rightarrow\mathcal{X}$ is the gradient descent map (dependent on the parameter $\beta$). It can be further shown that if $\beta$ is sufficiently small, then the operator $T$ is a contraction on $\mathcal{X}$, endowed with the Euclidean norm. We note here that the above arguments would be true if $\{L_{i}\}_{i=1}^{N}$ is a collection of strongly convex and smooth loss functions.

In practice, the exact gradient computation of loss function $L$ is computationally expensive as it requires evaluating $N$ gradients at every time step. Therefore, to ease the computational burden, a mini-batch Stochastic Gradient Descent (SGD) is employed. In the mini-batch SGD, at every step $k$, the gradient is approximated by a small, randomly sampled, subset (of size $n$) of the data set. To introduce this algorithm, let $\mathcal{N}_{k}\subset\{1,\ldots,N\}$ be the randomly sampled subset of size $n$. The state is updated as

[TABLE]

Note that $\hat{T}^{n}_{k}$ is now a random operator.

Remark 3.1**.**

Any realization of this random operator with small values of $n$ need not be a contraction since $\nabla_{x}L_{j}$ is a rank 1 positive semidefinite matrix. One could add a regularizer to the loss function to make it strongly convex. In particular, if the loss function is chosen as

[TABLE]

then $\hat{T}^{n}_{k}$ is a contraction operator for a sufficiently small $\beta$ irrespective of the $n$ used. $\Box$

Some obvious properties of these random operators are:

$\hat{T}^{n}_{k}$ is continuous in $\hat{z}^{n}_{k}$. 2. 2.

For every $\epsilon>0$ and $x\in\mathcal{X}$, we have

[TABLE] 3. 3.

Suppose that for any compact set $\mathcal{K}\subset\mathcal{X}$, $\|\nabla_{x}L_{i}(x)\|_{2}$ is uniformly bounded, that is, there exists $M_{\mathcal{K}}\in\mathbb{R}$ such that $\sup_{x\in\mathcal{K}}\sup_{1\leq i\leq N}\|\nabla_{x}L_{i}(x)\|_{2}\leq M_{\mathcal{K}}$. Then, for any $\epsilon>0$,

[TABLE]

This statement can readily be proved using the Hoeffding inequality and the union bound.

Let us depart from the specific case of logistic regression, and consider the case where $L_{i}$ can be any strictly concave function for all $i\in\{1,\ldots,N\}$. Then, $\hat{T}^{n}_{k}$ satisfies the following property.

If all the eigenvalues of the Hessian of $L_{i}(x)$ satisfy $0<m\leq\lambda(\nabla^{2}_{xx}L_{i}(x))\leq M<\infty$ for all $i\in\{1,\ldots,N\}$, then every realization of the random operator $\hat{T}^{n}_{k}$ is a contraction map with contraction coefficient $\hat{\alpha}^{n}_{k}\equiv\max\{|1-\beta M|,|1-\beta m|\}<1$ for an appropriately picked $\beta>0$.

We now introduce the empirical dynamic programming algorithm in the context of value iteration for MDP with discounted cost criteria.

3.2 Empirical Value Iteration for Discounted Cost MDP

Consider a Markov Decision Problem (MDP) problem described by 4 tuple $(\mathcal{S},\mathcal{A},c,p)$, where $\mathcal{S}$ is the finite state space, $\mathcal{A}$ is the finite action space, and $c:\mathcal{S}\times\mathcal{A}\longrightarrow\mathcal{R}$ is the cost function. The state transitions according to $s_{t+1}\sim p(\cdot|s_{t},a_{t})$. Let $\Gamma$ denote the set consisting of all possible deterministic policies $\gamma:\mathcal{S}\longrightarrow\mathcal{A}$. The infinite horizon discounted cost $v^{\gamma}:\mathcal{S}\longrightarrow\mathcal{R}$ starting from state $s$ and following policy $\gamma$ is given by

[TABLE]

where $\alpha\in(0,1)$ is the discount factor. The goal is to compute the optimal value $v^{*}(s)=\inf_{\gamma\in\Gamma}v^{\gamma}(s)$. Let $\mathcal{V}$ be the set of all $v:\mathcal{S}\longrightarrow\mathbb{R}$; this space is isomorphic to the Euclidean space $\mathbb{R}^{|\mathcal{S}|}$.

It can be shown that the optimal infinite horizon discounted cost is a fixed point of a contraction map $T:\mathcal{V}\rightarrow\mathcal{V}$, where $T$ is the Bellman operator given by

[TABLE]

Due to the Banach contraction mapping theorem, $T:\mathcal{V}\rightarrow\mathcal{V}$ admits a unique fixed point, which is equal to $v^{*}$. The iterative process of finding this unique fixed point is called the Value Iteration algorithm:

[TABLE]

In data driven applications, it is often the case that for all possible state-action pairs, multiple realizations of the next states are available. In this situation, we can replace the expectation in (8) to the empirical average. This algorithm is referred to as empirical dynamic programming, and is written as

[TABLE]

where $s_{k,i}^{\prime}(s,a)$ are $n$ independent and identically distributed samples of the next state given the current state-action pair $(s.a)$, redrawn at every $k$ independently from the past realizations. The above intuition can be turned into an algorithm to determine an approximately optimal value function, and is known as the empirical value iteration algorithm:

[TABLE]

Note that $\hat{T}^{n}_{k}$ is a random operator, and its realization is dependent on the samples generated $(s_{k,i}^{\prime}(s,a))_{i=1}^{n}$. The following properties of $\hat{T}^{n}_{k}$ are obvious:

$\hat{T}^{n}_{k}$ is a contraction with contraction coefficient $\alpha$. Therefore, $\hat{T}^{n}_{k}$ is continuous. 2. 2.

For every $\epsilon>0$ and $v\in\mathcal{V}$, we have

[TABLE] 3. 3.

In fact, we have a stronger property here. For any compact set $\mathcal{K}\subset\mathcal{V}$, we have for every $\epsilon>0$

[TABLE]

where $k=\max_{v\in\mathcal{K}}\|v\|_{\infty}$. This immediately yields

[TABLE] 4. 4.

Let $\mathcal{V}$ be endowed with the partial order $\preceq$ such that $v_{1}\preceq v_{2}$ implies $v_{1}(s)\leq v_{2}(s)$ for all $s\in\mathcal{S}$. If $c\geq 0$, then $\hat{T}^{n}_{k}$ satisfies

1. (a)

If $v_{0}=0$, then $v_{0}\preceq\hat{T}^{n}_{k}(v_{0})$. 2. (b)

If $v_{1}\preceq v_{2}$, then $\hat{T}^{n}_{k}(v_{1})\preceq\hat{T}^{n}_{k}(v_{2})$. 3. (c)

If $v_{l}\uparrow v$, then $\hat{T}^{n}_{k}(v_{l})\rightarrow\hat{T}^{n}_{k}(v)$ as $l\rightarrow\infty$.

3.3 Observations

Through the two examples above, we observed that the approximate operator $\hat{T}^{n}_{k}$ corresponding to the contraction operator $T$ is context dependent. In the case of stochastic gradient descent, it is constructed by picking certain loss functions randomly and then averaging their gradients. In the case of empirical dynamic programming, the approximate operator involves computing the empirical average of the future expected value. Nonetheless, there are some fundamental properties that the empirical operator satisfies in both situations. For instance, the property stated in (7) in the context of logistic regression is (mathematically) the same as the property stated in (12) in the context of empirical value iteration. Similarly, every realization of the random operator $\hat{T}^{n}_{k}$ is a contraction map under certain reasonable assumptions. We will consider more examples in Section 7, where we show that these properties (or some minor variant of these properties) are enjoyed by other empirical dynamic programming algorithms as well.

The other important observation is that every realization of random operator $\hat{T}^{n}_{k}$ may also satisfy some other desirable properties. For instance, in the empirical value iteration example, if we endow $\mathcal{V}$ with a partial order $\preceq$ and the cost is nonnegative, then every realization of $\hat{T}^{n}_{k}$ satisfies certain monotonicity properties. This is very useful in establishing the existence of unique invariant measure, as we show in Section 5. This property is, unfortunately, not satisfied by the logistic regression problem. This property is also not satisfied by the empirical relative value iteration for the average cost MDP. However, we will show that the realizations of the random operators in these cases have some other desirable properties that lead to the existence and uniqueness of the invariant measure.

We now turn our attention to introducing our first main result in the next section.

4 Weak* Convergence of the Distribution of Trajectories

We now study the convergence property of the sequence of distributions $\mu_{n}$, which is defined in (3). Before we study that, we need to ensure that the random operator $\hat{T}^{n}_{k}$ is “close to” the operator $T$ in some sense. Accordingly, we make the following assumption.

{assumption}

For every compact set $K\subset\mathcal{X}$, $\epsilon>0$, and $\delta>0$, there exists $N_{\epsilon,\delta}(K)>0$ such that

[TABLE]

We recall here that this assumption is satisfied by the logistic regression and empirical value iteration for discounted cost MDP (see (7) and (12) within the discussion at the end of Subsections 3.1 and 3.2). We are now in a position to introduce our first main result.

Theorem 4.1**.**

If Assumptions 2 and 4 hold, then $\mu_{n}$ converges in weak topology to $\psi$ as $n\rightarrow\infty$, where $\psi$ is defined in (4). *

Proof 4.2**.**

*The proof is based on the proof by [49], except that our hypotheses are slightly different from those in [49]. For completeness, we present a proof in Section 8. *

Levy-Prohorov’s metric over the space of probability measures over Polish spaces metrizes the weak* topology [60]. Generally speaking, if distributions of two random variables are close to each other in Levy-Prohorov’s metric, then it does not imply that the random variables will be close to each other. As an instance, the Levy-Prohorov’s metric between the measures of two independent and identically distributed random variables is zero, but the difference between the random variables themselves is not zero. If one of the random variables is deterministic (that is, its distribution is a Dirac mass), then the random variable must be close to the deterministic variable with high probability. This is established in the lemma below.

Let $\mathcal{A}$ be a Polish space with metric $\rho_{\mathcal{A}}$. Let $d_{P}$ be the Levy-Prohorov’s metric on the space of probability measures $\wp(\mathcal{A})$ over $\mathcal{A}$. This metric is defined as follows. For a Borel set $A\subset\mathcal{A}$, let $A^{\epsilon}$ be defined as

[TABLE]

Let $\mu,\nu\in\wp(\mathcal{A})$. Then, $d_{P}(\mu,\nu)$ is defined by

[TABLE]

We are now in a position to introduce our next result. We believe that this result may not be new, but we could not locate a reference where this result is proved.

Lemma 4.3**.**

Let $\mu\in\wp(\mathcal{A})$ and $\mathds{1}_{\{a^{*}\}}$ be a unit mass at point $a^{*}\in\mathcal{A}$. If $d_{P}(\mu,\mathds{1}_{\{a^{*}\}})<\epsilon$, then for any random variable $W$ distributed according to the law $\mu$, we have

[TABLE]

Proof 4.4**.**

Let $B_{\epsilon}$ be an open $\epsilon$ ball around $a^{*}$. Let $(B_{\epsilon}^{\complement})^{\epsilon}$ be defined as

[TABLE]

Note that $(B_{\epsilon}^{\complement})^{\epsilon}=\mathcal{A}\setminus a^{*}$, which implies $\mathds{1}_{\{a^{*}\}}((B_{\epsilon}^{\complement})^{\epsilon})=0$. Let $W$ be a random variable distributed according to the law $\mu$. Then, from the definition of Levy-Prohorov’s metric, we know that

[TABLE]

The proof then follows from noting that

[TABLE]

*The proof is established. *

As a consequence of the lemma above, we conclude that since the distribution $\mu_{n}$ converges to the Dirac delta function $\psi$, it implies that for $n$ sufficiently large, the random sequence generated by RSA lies within a small tube around the trajectory induced by the deterministic contraction operator with high probability. This is illustrated in the Figure 1.

To see this, let $\hat{\boldsymbol{z}}^{n}:=(\hat{z}^{n}_{0},\hat{z}^{n}_{1},\ldots)$ and $\boldsymbol{y}=(y_{0},y_{1},\ldots)$. Endow the space $\mathcal{X}^{\mathbb{N}}$ with the following metric:

[TABLE]

It can be readily established that $\rho_{\mathcal{X}^{\mathbb{N}}}$ defined above is a metric on $\mathcal{X}^{\mathbb{N}}$. Then, due to Lemma 4.3, we conclude that

[TABLE]

Next, note that if $\hat{\boldsymbol{z}}^{n}$ satisfies $\rho_{\mathcal{X}^{\mathbb{N}}}(\hat{\boldsymbol{z}}^{n},\boldsymbol{y})<\epsilon$, then $\hat{z}^{n}_{k}$ is within $\epsilon$ neighborhood of $y_{k}$ for most of the $k\in\mathbb{N}$.

5 Existence of Invariant Measures for Fixed $n$

In this section, we identify conditions on the operators $(\hat{T}^{n}_{k})_{k\in\mathbb{N}}$ so that the Markov chain $\hat{z}^{n}_{k}$ admits a unique invariant distribution, denoted by $\pi^{n}$, as $k\rightarrow\infty$. The existence of an invariant measure yields insight about “stability” of an RSA. In particular, if there is no invariant distribution, then it is likely the case that the sequence generated by RSA can blow up with positive probability. Thus, by tweaking the RSA (for instance, by changing the stepsize or increasing the number of samples), one can ensure that the sufficient conditions noted below are satisfied, thereby establishing that the RSA is stable and yields finite values with probability 1.

In the case where there exists unique invariant measure under certain assumptions on the initial condition, then it means that any element in the tail of the random sequence generated by the RSA will have its law as the invariant distribution. This is a crucial step in proving that the time average of $f(\hat{z}^{n}_{k})$ for any $f\in C_{b}(\mathcal{X})$ converges in probability to the spatial average $\int fd\pi^{n}$ of the function with respect to the invariant measure $\pi^{n}$. This important result is established in the next section.

To state the assumptions, we drop the subscript $k$ wherever possible since the statistical properties of $\hat{T}^{n}_{k}$ and $\hat{T}^{n}_{k^{\prime}}$ are independent and identical to each-other as long as $k\neq k^{\prime}$. Below, we list three assumptions under which we can show that $\hat{z}^{n}_{k}$ admits an invariant distribution.

To introduce our first assumption, we need to assume a partial order $\preceq$ exists on the normed space $\mathcal{X}$ (we will drop the completeness requirement on the space $\mathcal{X}$ for this assumption). For instance, $\mathcal{X}$ could be the Euclidean space with the natural partial order, wherein $x,y\in\mathbb{R}^{n}$, $x\preceq y$ implies $x_{i}\leq y_{i}$ for all $i\in\{1,\ldots,n\}$. The limit of an increasing sequence $(x_{k})_{k\in\mathbb{N}}$ satisfying $x_{1}\preceq x_{2}\preceq x_{3}\preceq\ldots$ under this partial order may escape to infinity. Thus, it is required to bound the space $\mathcal{X}$ to ensure that the sequences do not diverge along any of the coordinates. Another example of a normed space with partial order is the space of measurable functions from a Euclidean space $\mathbb{R}^{n}$ to $[-M,M]$, denoted by $\mathcal{L}_{\infty}(\mathbb{R}^{n},[-M,M])$ (where $M$ is a fixed positive real number). This space features a natural partial order, wherein $f_{1}\preceq f_{2}$ iff $f_{1}(s)\leq f_{2}(s)$ for every $s\in\mathbb{R}^{n}$.

{assumption}

The following conditions are satisfied:

$\mathcal{X}$ is a bounded normed space (not necessarily Polish) with partial order $\preceq$. This ordering satisfies the following property: For any sequence $(x_{k})_{k\in\mathbb{N}}$ satisfying $x_{1}\preceq x_{2}\preceq x_{3}\preceq\ldots$, there exists a minimal element $\bar{x}\in\mathcal{X}$ such that $x_{k}\preceq\bar{x}$ for all $k\in\mathbb{N}$. This is denoted as $x_{k}\uparrow\bar{x}$. 2. 2.

The operator $\hat{T}^{n}$ satisfy

1. (a)

Monotonicity 1: exists $x_{0}\in\mathcal{X}$ such that $x_{0}\preceq\hat{T}^{n}(x_{0})$ almost surely. 2. (b)

Monotonicity 2: If $x_{1}\preceq x_{2}$, then $\hat{T}^{n}(x_{1})\preceq\hat{T}^{n}(x_{2})$ almost surely. 3. (c)

Continuity: If $x_{k}\uparrow x$, then $\hat{T}^{n}(x_{k})\rightarrow\hat{T}^{n}(x)$ as $k\rightarrow\infty$ almost surely.

Assumption 5 is satisfied in Markov decision processes with non-negative cost functions. This has been noted in Subsection 3.2 in the context of empirical value iteration for discounted cost criterion. However, this is also satisfied in MDP for total cost criterion with an absorbing state and having a proper policy (that is, there is a stationary policy under which the MDP terminates in an absorbing state with probability 1). The proof of the last claim follows from arguments similar to the one made in Subsection 3.2; see, for example, Chapter 1 of [9], and standard texts on MDPs with total cost criterion [40, 41, 8].

For RSAs involving MDPs, the point $x_{0}$ can be chosen easily. Assuming that the costs are non-negative, one could take the zero function as the initial value function or Q function – and the Monotonicity 1 condition is satisfied almost surely.

{assumption}

For $n\in\mathbb{N}$ and $m\in\mathbb{N}$, let $\hat{\alpha}_{m}$ denote the Lipschitz coefficient of $\hat{T}^{n}_{m}\circ\hat{T}^{n}_{m-1}\circ\ldots\circ\hat{T}^{n}_{1}$. The following conditions are satisfied:

$\mathcal{X}$ is a compact Polish space. 2. 2.

For $n\in\mathbb{N}$, there exists $m\in\mathbb{N}$, such that for any $\epsilon>0$, there exists a $\delta_{\epsilon}>0$ such that

[TABLE]

The notable point in Assumption 5 is that the assumption requires $\mathcal{X}$ to be compact. It is satisfied in empirical value iteration for MDP with average cost criterion, as long as we project the value functions outside a sufficiently large compact set back to that compact set. We adopted this approach earlier in [34] to ensure that the value functions obtained through repeated use of empirical operators do not blow up. During simulations, however, we never needed to use projection, as the value functions were bounded.

{assumption}

The following conditions are satisfied:

$\mathcal{X}$ is a Polish space. 2. 2.

There exists $a,b>0$ such that the operator $\hat{T}^{n}$ satisfy

[TABLE]

where $x^{*}$ is the fixed point of $T$. 3. 3.

Let $\hat{\alpha}^{n}$ denote the Lipschitz coefficient of $(\hat{T}^{n})_{k\in\mathbb{N}}$. Then,

[TABLE]

Assumption 5 is satisfied in stochastic gradient descent of strongly convex and smooth functions as noted in Subsection 3.1. This is also trivially satisfied in the empirical value iteration for an MDP with a discounted reward criterion, as we have noted in Subsection 3.2.

Remark 5.1**.**

*Assumption 5(ii) is typically proved using concentration of measures results, like Hoeffding inequality [42, 52, 64] or the theory of empirical processes [61, 81]. Some references where such inequalities have been used in empirical dynamic programming for continuous state MDPs are [57, 38, 33]. *

Our next theorem summarizes the main result of this section. Let $\mu^{n}_{k}$ denote the distribution of $\hat{z}^{n}_{k}$.

Theorem 5.2**.**

Suppose that Assumption 2 holds. Additionally, if either one of three assumptions, Assumptions 5, 5, 5, holds, then there exists an invariant measure $\pi^{n}$ such that $\mu^{n}_{k}$ converges to $\pi^{n}$ in weak topology. Further, the invariant measure is unique under either of the following circumstances:*

Assumption 5 holds, and the RSA is always initialized with $x_{0}$, where $x_{0}$ is defined in Assumption 5. 2. 2.

Assumption 5 holds. 3. 3.

Assumption 5 holds.

Proof 5.3**.**

*Under Assumption 5, the existence of invariant measures is proved in [17, Theorem 8.1, p. 79-81]. Under Assumption 5, the existence of invariant measures is proved in [12] and [17, Theorem 8.2, p. 82-83]. Under Assumption 5, the existence and uniqueness result is established in [25, Theorem 1.1, p. 87]. *

Remark 5.4**.**

We can replace the assumption of $\mathcal{X}$ being a compact Polish space in Assumption 5 by making the following assumption. There exists $x_{0}\in\mathcal{X}$ such that if $\hat{z}^{n}_{0}=x_{0}$, then for any $\delta>0$, there exists $N_{\delta}\in\mathbb{N}$ such that for any $k\geq 1$, we have

[TABLE]

*If the above condition and Assumption 5 (2) holds, then one can show that for any initial condition $\hat{z}^{n}_{0}\in\mathcal{X}$, a unique invariant distribution exists. For a proof, we refer the reader to [17, p. 179]. However, proving (13) is satisfied in usual RSAs appears to be difficult in our experience. $\Box$ *

We now turn our attention to establishing the law of large numbers for time averages of the outputs from a RSA.

6 Averaging of Iterates and The Weak Law of Large Numbers

RSAs with constant stepsize and averaging of iterates have received significant attention recently. Nemirovski et. al. in [58] develops an algorithm for stochastic gradient descent with averaging of the last few iterates and shows that the algorithm is robust to stepsize selection and yields better results when properties of the loss functions (such as strong convexity parameter, Lipschitz coefficient of the derivative, etc.) are unknown. Along similar lines, [5, 4, 27] study stochastic gradient descent based algorithms with constant stepsize and averaging and derive the finite time guarantees on the loss achieved with averaged output. We note here that the constant stepsize algorithms with averaging are fundamentally different from Polyak-Ruppert averaging used in stochastic approximation [62, 50], where the stepsizes decreases at a rate slower than $1/k$.

Motivated by the above references, in this section, we consider the problem of convergence of the sequence of averages of the random sequence $(\hat{z}^{n}_{k})_{k\in\mathbb{N}}$ (or of $(f(\hat{z}^{n}_{k}))_{k\in\mathbb{N}}$ for some $f\in C_{b}(\mathcal{X})$). As stated in Subsection 1.1, averaging of the last few outputs of an RSA is used to reduce the variance in the final output of the algorithm. Specifically, if we terminate the iteration of RSA at $k$, and output $\hat{z}^{n}_{k}$, there is a small chance that this output is far from $x^{*}$. It is generally believed that if we average the last few outputs (assuming none of them are $\infty$), then it reduces the variance in the output. We establish this result rigorously here, under the assumption that the Markov chain has a unique invariant distribution.

We will assume throughout that the Markov chain $(\hat{z}^{n}_{k})_{k\in\mathbb{N}}$ admits an invariant distribution (the precise assumption is stated in the sequel). The time-average of the Markov chain is precisely the law of large numbers for Markov chains. It has been studied within the context of Markov chains over compact spaces in [19] and over general locally compact spaces in [56, Sec 18.5]. Let us formulate the problem precisely.

Consider a continuous function $f\in C_{b}(\mathcal{X})$. Let $\pi^{n}$ denote the invariant measure of the Markov chain $(\hat{z}^{n}_{k})_{k\in\mathbb{N}}$. We have already studied the conditions under which such an invariant measure would exist in Theorem 5.2. In what follows, we show that under relatively mild assumptions, we have

[TABLE]

Note that in the expression above, the term $\int f(x)\pi^{n}(dx)$ is the spatial average of the function $f$ at the invariant distribution. Thus, what we show is that the time average converges in probability to the spatial average – which implies that $(f(\hat{z}^{n}_{k}))_{k\in\mathbb{N}}$ is a discrete time ergodic process.

Let us define the operator $F:C_{b}(\mathcal{X})\rightarrow C_{b}(\mathcal{X})$ and its adjoint $F^{*}:\wp(\mathcal{X})\rightarrow\wp(\mathcal{X})$, where $\wp(\mathcal{X})$ is endowed with the weak* topology, as follows:

[TABLE]

{assumption}

The distribution $\mu$ of the initial condition $\hat{z}^{n}_{0}$ is picked from a set $\mathcal{M}\subset\wp(\mathcal{X})$. Either of the following two conditions holds:

There exists a unique invariant measure $\pi^{n}$ such that for any $\mu\in\mathcal{M}$, $(F^{*})^{k}(\mu)$ converges in weak* topology to $\pi^{n}$. 2. 2.

There exists a unique invariant measure $\pi^{n}$ such that for any $\mu\in\mathcal{M}$, the averaged operator satisfies

[TABLE]

where $(F^{*})^{k}(\mu)$ denotes $k$ compositions of $F^{*}$ applied on $\mu$ and the convergence is in weak* sense.

From the Stolz-Cesaro theorem [20, Theorem 2.7.1, p.59], it is easy to show that if Assumption 6(1) holds, then Assumption 6(2) holds as well. However, the converse may not be true. To prove our next result, we only need Assumption 6(2) to hold. It should be noted that Assumption 6(1) may be rather easy to prove using Theorem 5.2.

Theorem 6.1**.**

*If Assumptions 2 and 6 hold, then (14) holds. *

Proof 6.2**.**

*The proof essentially follows the steps in [19] and [56, Theorem 18.5.1, p. 478], except that we relax the assumption on compactness of the state space as assumed in [19] and replace the hypotheses in [56] with Assumption 6. For completeness, a detailed proof is presented in Section 9. *

One way the presentation of Assumption 6 departs from the traditional Markov chain literature is as follows. It is generally assumed that $\mathcal{M}=\wp(\mathcal{X})$, that is, for every $\mu\in\wp(\mathcal{X})$, $(F^{*})^{k}(\mu)$ converges in weak* topology to $\pi^{n}$. This is a very strong assumption from the applicability viewpoint in RSAs. In particular, it is possible to pick the most suitable initialization for the RSAs, which implies that $\mathcal{M}$ can be picked appropriately. For example, in empirical dynamic programming for MDPs, one can initialize the value function to be 0. Then, we can utilize Assumption 5 (with $x_{0}=0$) to establish the existence of a unique invariant distribution using Theorem 5.2. Incidentally, for the law of large numbers to hold, we do not need the stronger condition of $\mathcal{M}=\wp(\mathcal{X})$.

We now have a corollary of the result above, which we capture as a theorem below due to its importance and applicability to RSAs.

Theorem 6.3**.**

Suppose that $\mathcal{X}\subset\mathbb{R}^{n}$ is a compact set. Suppose further that either of the following conditions hold:

Assumption 5 holds, and the RSA is always initialized with $x_{0}$, where $x_{0}$ is defined in Assumption 5. 2. 2.

Assumption 5 holds. 3. 3.

Assumption 5 holds.

Consider the average $\hat{a}^{n}_{k}$ of the Markov chain

[TABLE]

*If Assumptions 2 holds, then for any initialization $\hat{z}^{n}_{0}\in\mathcal{X}$, a unique invariant distribution $\pi^{n}$ exists and $\hat{a}^{n}_{k}$ converges almost surely to the mean $\bar{a}^{n}$ of the distribution $\pi^{n}$. *

Proof 6.4**.**

*The existence of a unique invariant measure is due to Theorem 5.2. Since $\mathcal{X}$ is a compact set, we can take $f(x)=x_{i}$ to conclude that $(\hat{a}^{n}_{k})_{i}$ converges almost surely to $\bar{a}^{n}_{i}$ by [19]. *

The result in (14) requires $f$ to be a bounded function over $\mathcal{X}$. We now consider the case where $f$ is a potentially unbounded continuous function. To establish essentially the same result, we need to make the following additional assumption. {assumption} There exists a random variable $C(\omega)$ such that

[TABLE]

and $C(\omega)<\infty$ for $\mathbb{P}$-almost all $\omega\in\Omega$.

Theorem 6.5**.**

*Suppose that $f:\mathcal{X}\rightarrow\mathbb{R}$ is a continuous function (potentially unbounded). If Assumptions 2, 6, and 6 hold, then (14) holds. *

Proof 6.6**.**

*The proof is presented in Subsection 9.1. *

It is easy to see in practice if Assumption 6 holds or not, particularly when $\mathcal{X}=\mathbb{R}^{n}$. If the random iterates are uniformly bounded during a single run of the RSA, then Assumption 6 holds along that trajectory. If almost all independent runs of the RSA is not expected to “blow up”, then the time average of the iterates is likely going to have the variance reduction property–the time average converges in probability to the spatial average.

We now turn our attention to illustrating the application of Theorems 4.1, 6.1, and 6.5 in various optimization and empirical dynamic programming algorithms.

7 Numerical Simulations

In this section, we complement the theoretical results proved above with extensive numerical simulations. We conduct simulations of minibatch stochastic gradient descent for logistic regression on MNIST dataset, empirical value iteration for discounted and average cost MDPs, and empirical Q-value iteration.

7.1 Stochastic Gradient Descent

Consider the task of classifying a subset of MNIST handwritten digits, where we consider only the images corresponding to the numbers 0 and 1. Each data point is a $28\times 28$ pixel image with the corresponding label (either [math] or $1$). We use logistic regression and Poisson regression with an $\ell_{2}$ regularizer for the classification task. We refer the reader to Subsection 3.1 for details of this problem for the logistic regression. The loss function for the logistic regression with the regularizer is

[TABLE]

The loss function for Poisson regression with regularizer is

[TABLE]

The first and the second derivatives of this loss function is

[TABLE]

We use here the notations introduced in Subsection 3.1. We transform each image into a vector and append $1$ at the beginning of the vector. Thus, the space $\mathcal{U}=\{1\}\times[0,1]^{784}$. Thus, the space $\mathcal{X}=\mathbb{R}^{785}$. As mentioned previously, the variable $n$ represents the batch size picked at every SGD iteration step. We pick $y_{0}$ arbitrarily in $\mathcal{X}$ and set $\hat{z}^{n}_{0}=y_{0}$. Then, we run the exact gradient descent and the minibatch SGD as follows:

[TABLE]

As discussed in Subsection 3.1, it is clear that the exact gradient descent is a contraction map for stepsize $\beta$ small enough. Therefore, $y_{k}$ converges to the optimal solution $x^{*}$. We can make the following claim about the operator $\hat{T}^{n}_{k}$:

Theorem 7.1**.**

The random operator $\hat{T}^{n}_{k}$ satisfies Assumptions 2, 4. If $\beta$ is sufficiently small, then $\hat{T}^{n}_{k}$ also satisfies Assumption 5. Let $\mu_{n}$ denote the distribution of $(\hat{z}^{n}_{0},\hat{z}^{n}_{1},\ldots)$ and $\psi$ be the Dirac mass on $(y_{0},y_{1},\ldots)$. We have $\{\mu_{n}\}_{n\in\mathbb{N}}$ converges to $\psi$ in weak topology as $n\rightarrow\infty$ and there exists a unique invariant distribution of the Markov chain $(\hat{z}^{n}_{k})_{k\in\mathbb{N}}$ for any $n\in\mathbb{N}$. *

Proof 7.2**.**

*The first conclusion follows from the discussion leading to (7) in Subsection 3.1. For $\beta$ sufficiently small, every realization of $\hat{T}^{n}_{k}$ is a contraction as discussed in Remark 3.1. We only need to show that Part 2 of Assumption 5 is satisfied. Since $n$ is finite, $\hat{T}^{n}_{k}(x^{*})$ can take only finitely many values, and therefore Part 2 of Assumption 5 is trivially satisfied. The existence of invariant distribution then follows from Theorem 5.2. *

As a result of the theorem above, we conclude that $\|y_{k}-\hat{z}^{n}_{k}\|_{2}$ is generally small for most $k$, and its variance converges to 0 as $n$ increases. This is evident in Figures 2 and 4, where we plot the distance $\|y_{k}-\hat{z}^{n}_{k}\|_{2}$ for various values of $n$ (batch sizes) and $k$ ranging from 0 to 5000. The black curve is the mean $\mathds{E}\left[\|y_{k}-\hat{z}^{n}_{k}\|_{2}\right]$ and the red region shows one standard deviation of the distance $\|y_{k}-\hat{z}^{n}_{k}\|_{2}$ over 1000 sample paths of the minibatch SGD iterations. We picked the same initial condition $\hat{z}^{n}_{0}=y_{0}=0$ for all the sample paths to generate the figure.

Figures 3 and 5 plots $\|\hat{a}^{n}_{k}-x^{*}\|_{2}$, and we observe that the variance of time averages reduces as $k$ grows large for any batch size $n$. This variance reduction of time averages can be attributed to the contraction property of the random maps $\hat{T}^{n}_{k}$, which in turn is due to the addition of a strongly convex regularizer to the empirical loss function.

7.2 Empirical Value Iteration for Discounted Cost MDP

We consider here the empirical value iteration for discounted Markov decision processes (MDP) as described in Subsection 3.2. Consider the value iteration algorithm applied to an MDP in which there are 20 states and 5 actions. We generate the state transition probability matrix for this MDP randomly at the beginning of the simulation.

We use here the notation introduced in Subsection 3.2. We initialize $v_{0}$ arbitrarily and set $\hat{z}^{n}_{0}=v_{0}$, and define the iterates of exact value iteration and empirical value iteration as

[TABLE]

We can prove the following result.

Theorem 7.3**.**

The random operator $\hat{T}^{n}_{k}$ satisfies Assumptions 2, 4, and 5. As a result, we conclude:

Let $\mu_{n}$ denote the distribution of $(\hat{z}^{n}_{0},\hat{z}^{n}_{1},\ldots)$ and $\psi$ be the Dirac mass on $(y_{0},y_{1},\ldots)$. We have $\{\mu_{n}\}_{n\in\mathbb{N}}$ converges to $\psi$ in weak topology as $n\rightarrow\infty$.* 2. 2.

There exists a unique invariant distribution $\pi^{n}$ of the Markov chain $(\hat{z}^{n}_{0},\hat{z}^{n}_{1},\ldots)$. Further, the sequence of distribution of $\hat{z}^{n}_{k}$ converges to this invariant distribution $\pi^{n}$. 3. 3.

The time average $\hat{a}^{n}_{k}$ converges in probability to the mean of $\pi^{n}$ as $k\rightarrow\infty$.

Proof 7.4**.**

The proof follows from the discussions in Subsection 3.2, where we show that $\hat{T}^{n}_{k}$ satisfies Assumptions 4 and 5.

*Moreover, since every realization of $\hat{T}^{n}_{k}$ is a contraction operator with coefficient $\alpha$, if $\|v\|_{\infty}\leq\frac{\|c\|_{\infty}}{1-\alpha}$, then $\|\hat{T}^{n}_{k}(v)\|_{\infty}\leq\frac{\|c\|_{\infty}}{1-\alpha}$. Thus, there is no loss of generality in restricting $\mathcal{X}$ to be in the compact set $\{v\in\mathbb{R}^{|\mathcal{S}|}:\|v\|_{\infty}\leq\frac{\|c\|_{\infty}}{1-\alpha}\}$. The fact that the time averages converge in probability to the mean of the invariant distribution is a direct application of Theorem 6.3. *

Figure 6 shows the mean (black line) and the variance (red region) of $\|v_{k}-\hat{z}^{n}_{k}\|_{\infty}$ for every time step $k=1,\ldots,10000$ for $n=1,25,400$. To compute the mean and the variance, 1000 sample paths were taken. As proved in Theorem 7.3 above, the mean of $\|v_{k}-\hat{z}^{n}_{k}\|_{\infty}$ becomes smaller as $n$ grows. Figure 7 shows the mean and the variance of $\|v^{*}-\hat{a}^{n}_{k}\|_{\infty}$ for $n=1,25,400$ for $k=1,\ldots,10000$. For every $n$, since the sequence of distribution of $\hat{z}^{n}_{k}$ converges to a unique invariant distribution, $\hat{a}^{n}_{k}$ converges in probability to the mean of the invariant distribution. We observe in the figure that the variance of $\|v^{*}-\hat{a}^{n}_{k}\|_{\infty}$ reduces as $k$ grows implying that the average $\hat{a}^{n}_{k}$ is close but not equal to $v^{*}$. We further observe that the variance of $\|v^{*}-\hat{a}^{n}_{k}\|_{\infty}$ is vanishingly small for large values of $n$ – a result that we have not formally proved here, and can be addressed in a future work.

7.3 Synchronous Batch Q-Value Iteration for Discounted Cost MDP

Q-value iteration is another algorithm that, like value iteration, computes the optimal value function in MDPs. Let $\mathcal{Q}$ denote the set of all Q-value functions $Q:\mathcal{S}\times\mathcal{A}\rightarrow\mathbb{R}$. Similar to Bellman operator of value iteration, we define an operator $T:\mathcal{Q}\rightarrow\mathcal{Q}$ as

[TABLE]

Similar to Bellman operator, $T$ is a contraction on $(\mathcal{Q},\|\cdot\|_{\infty})$. Further, it can be shown that the fixed point of $T$ is $q^{*}$, which is defined as $q^{*}(s,a)=c(s,a)+\alpha\mathds{E}\left[v^{*}(s^{\prime})|s,a\right]$, where $v^{*}(\cdot)=\min_{a\in\mathcal{A}}q(\cdot,a)$. The Q-value iteration starts with an arbitrary $q_{0}\in\mathcal{Q}$ and generates the sequence according to $q_{k+1}=T(q_{k})$, which converges to $q^{*}$ as $k\rightarrow\infty$.

The exact operator, as in other cases considered in the paper, can be approximated by the empirical operator $\hat{T}^{n}_{k}$:

[TABLE]

where $\{s^{\prime}_{k,i}(s,a)\}_{i=1}^{n}$ are $n$ i.i.d. samples of the next state given the current state-action pair $(s,a)$.

Let us define $\hat{z}^{n}_{k+1}=\hat{T}^{n}_{k}(\hat{z}^{n}_{k})$, where $\hat{z}^{n}_{0}=q_{0}$. Let $\hat{a}^{n}_{k}$ be the time averaged version of $\hat{z}^{n}_{k}$. The properties of the random operator $\hat{T}^{n}_{k}$ for empirical Q-value iteration has the same properties as listed in Theorem 7.3 for the case of empirical value iteration. Thus, we omit repetition of essentially the same result here. The simulation results are plotted in Figure 8 and 9.

Figure 8 shows the mean (black line) and the variance (red region around the mean) of $\|q_{k}-\hat{z}^{n}_{k}\|_{\infty}$ for various values of $n$ for $k$ ranging from $1$ to $10000$. This mean and variance is computed using 1000 sample paths of the algorithm starting from the same initial $\hat{z}^{n}_{0}=q_{0}$ for all sample paths. As is evident, the mean reduces as we increase the sample size $n$.

Figure 9 plots the mean and the variance of $\|q^{*}-\hat{a}^{n}_{k}\|$. Once again, we observe that the variance reduces as $k$ increases since the Markov chain admits a unique invariant measure.

8 Proof of Theorem 4.1

To prove Theorem 4.1, we need to introduce some further notation. Define ${\mathbb{Z}_{+}}=\{0,1,2,\ldots\}$ and we use $\mu_{n}\Rightarrow\psi$ to denote the convergence in weak topology. Let $\Pi_{k}:\mathcal{X}^{\mathbb{Z}_{+}}\rightarrow\mathcal{X}^{k+1}$ denote the projection operator that projects a sequence to its first $(k+1)$ components, that is,

[TABLE]

For a measure $\mu\in\wp(\mathcal{X}^{\mathbb{Z}_{+}})$, let $\mu\circ\Pi_{k}^{-1}\in\wp(\mathcal{X}^{k+1})$ denote the pullback of the measure to the first $k+1$ components. We note the following fact from probability theory.

Theorem 8.1** ([13], Theorem 2.8).**

*Let $(\nu_{n})_{n\in\mathbb{N}}\subset\wp(\mathcal{X}^{\mathbb{Z}_{+}})$ be a sequence of measures. Then, $(\nu_{n})_{n\in\mathbb{N}}\Rightarrow\nu_{\infty}$ if and only if $\nu_{n}\circ\Pi_{k}^{-1}\Rightarrow\nu_{\infty}\circ\Pi_{k}^{-1}$ for every $k\in{\mathbb{Z}_{+}}$. *

Recall that $\mu_{n}$ is a measure defined on a stochastic sequence and we would like to show that $\mu_{n}\Rightarrow\psi$. Due to the theorem above, all we need to establish is that $\mu_{n}\circ\Pi_{k}^{-1}\Rightarrow\psi\circ\Pi_{k}^{-1}$ for every $k$. We will establish this result via an induction argument following [49].

Fix $n\in\mathbb{N}$. For $k=0$, let $\mu_{n}\circ\Pi_{0}^{-1}(F)=\mathds{1}_{\{y_{0}\}}(F)=\psi\circ\Pi_{0}^{-1}(F)$ for all measurable $F\subset\mathcal{X}$. Thus, the statement is true for $k=0$. In the induction step, we prove that if $\mu_{n}\circ\Pi_{k-1}^{-1}\Rightarrow\psi\circ\Pi_{k-1}^{-1}$ and for every closed set $F_{1}\subset\mathcal{X}^{k}$ and $F_{2}\subset\mathcal{X}$, we have

[TABLE]

By the Portmanteau theorem [13, Theorem 2.1], we then conclude that $\mu_{n}\circ\Pi_{k}^{-1}\Rightarrow\psi\circ\Pi_{k}^{-1}$ for every $k\in{\mathbb{Z}_{+}}$. We divide the proof into three steps.

Step 1: We show that for any uniformly continuous function $g$, the map $g\circ T$ is also uniformly continuous since $T$ is Lipschitz operator.

Lemma 8.2**.**

*If $g\in U_{b}(\mathcal{X})$, then $g\circ T\in U_{b}(\mathcal{X})$. *

Proof 8.3**.**

Fix $\epsilon>0$. Since $g\in U_{b}(\mathcal{X})$, we can pick a $\delta_{\epsilon}>0$ such that for any $x,x^{\prime}\in\mathcal{X}$ with $\rho(x,x^{\prime})<\delta_{\epsilon}$, $|g(x)-g(x^{\prime})|<\epsilon$. Now, for any $y,y^{\prime}\in\mathcal{X}$ with $\rho(y,y^{\prime})<\delta_{\epsilon}$, we know that

[TABLE]

Taking $x=T(y)$ and $x^{\prime}=T(y^{\prime})$, we conclude that

[TABLE]

*This implies that $g\circ T\in U_{b}(\mathcal{X})$. *

Step 2: Next, we establish a consequence of Assumption 4 and Lemma 8.2. Let $x_{0:l}=(x_{0},\ldots,x_{l})$. We assume that $\mu_{n}\circ\Pi_{k}^{-1}\Rightarrow\psi\circ\Pi_{k}^{-1}$ for all $k=0,\ldots,l$, which serves as the induction hypothesis.

Lemma 8.4**.**

Suppose that Assumption 4 holds. Then, for every $\epsilon>0$, there exists an $N_{\epsilon}$ such that

[TABLE]

and

[TABLE]

Proof 8.5**.**

*Lemma 8.2 implies that $g\circ T$ is uniformly continuous. The first result is a direct consequence of the induction hypothesis. The proof of the second result is presented in Appendix A. *

Step 3: We now complete the proof by establishing the induction step using the result from Lemma 8.4. We again follow [49]. Let $F_{1}\subset\mathcal{X}^{l+1}$ be any closed set. Then, using Lemma 8.4, for any $g\in U_{b}(\mathcal{X})$, there exists $N_{\epsilon}$ such that

[TABLE]

Now, for any $F_{2}\subset\mathcal{X}$ closed, we can construct a sequence of $(g_{m})_{m\in\mathbb{N}}\subset U_{b}(\mathcal{X})$ such that $g_{m+1}\leq g_{m}$ for all $m\in\mathbb{N}$ and $g_{m}\rightarrow 1_{F_{2}}$. This leads us to the following inequality for every $m\in\mathbb{N}$:

[TABLE]

Taking the limsup on both sides and using (15), we arrive at the following inequality

[TABLE]

Now, since the right side holds for every $m\in\mathbb{N}$, we take the limit $m\rightarrow\infty$ and use the bounded convergence theorem to conclude

[TABLE]

Collecting the two inequalities in (16) and (17), we conclude that

[TABLE]

The proof of the lemma is complete.

9 Proof of Theorem 6.1

We first introduce some notations. For $f\in C_{b}(\mathcal{X})$, we define $g_{k}\in C_{b}(\mathcal{X})$, $k=0,1,\ldots$ as

[TABLE]

The following equation follows immediately from the above definitions:

[TABLE]

Define $\bar{F}^{*}_{k}=\frac{1}{k}\sum_{i=1}^{k}(F^{*})^{i}$. Define the constant function $c_{f}\in C_{b}(\mathcal{X})$ as

[TABLE]

The average of the functions $g_{k}$, denoted by $\bar{g}_{k}$, is

[TABLE]

For a function $f\in C_{b}(\mathcal{X})$ and a set $C\subset\mathcal{X}$, we use $f|_{C}$ to denote the restriction of the function on the set $C$. We now prove three lemmas that lead to the result. For the next result, let us define the occupation measure $\eta_{k}$ over the set $C\subset\mathcal{X}$ as

[TABLE]

We claim the following.

Lemma 9.1**.**

If Assumptions 2 and 6 holds, then for every $\epsilon>0$, there exists a compact set $C_{\epsilon}\subset\mathcal{X}$ such that

[TABLE]

Proof 9.2**.**

Note that since $\bar{F}^{*}_{k}(\mu)\rightarrow\pi^{n}$ for any $\mu\in\mathcal{M}$ by Assumption 6, we conclude from Prohorov’s theorem that the sequence $\{\bar{F}^{*}_{k}(\mu)\}_{k\in\mathbb{N}}$ is tight. Thus, for $\epsilon>0$, let $C_{\epsilon}$ be the compact set such that

[TABLE]

Further, we note that for any set $C\subset\mathcal{X}$, we have

[TABLE]

Using the above identity and using Markov’s inequality, we conclude that

[TABLE]

*The proof of the theorem is complete. *

Lemma 9.3**.**

*Let $C\subset\mathcal{X}$ be a compact set. Then, the sequence of functions $(\bar{g}_{k}|_{C})_{k\in\mathbb{N}}$ is uniformly bounded and equicontinuous and converges uniformly to $c_{f}|_{C}$. *

Proof 9.4**.**

First, we note that $\|g_{i}\|_{\infty}\leq\|f\|_{\infty}$ for all $i\in\mathbb{N}$, which implies that $\bar{g}_{k}|_{C}$ is uniformly bounded.

The proof of equicontinuity follows directly from Assumption 6(2) and Ascoli theorem. Note that as $k\rightarrow\infty$, we get

[TABLE]

*Thus, by Ascoli’s theorem, $\{\bar{g}_{k}|_{C}\}_{k\in\mathbb{N}}$ is an equicontinuous sequence of functions. The result then follows using Assumption 6(2). *

Next, we establish the following result.

Lemma 9.5**.**

For every $M\in\mathbb{N}$, we have

[TABLE]

Proof 9.6**.**

*See Appendix B. *

The proof can now be completed easily. Fix $\epsilon>0$ and recall the definition of the set $C_{\epsilon}$ from Lemma 9.1. We now note that for every $K\in\mathbb{N}$ and $M\in\mathbb{N}$, we have

[TABLE]

Since $\bar{g}_{M}(x)\rightarrow\int fd\pi^{n}$ uniformly on the compact set $C_{\epsilon}$ due to Lemma 9.3, we can pick $M_{\epsilon}$ sufficiently large such that for all $K\in\mathbb{N}$ and $M\geq M_{\epsilon}$, we get

[TABLE]

For such $M_{\epsilon}$, as $K\rightarrow\infty$, the first summand goes to 0 by Lemma 9.5. In the third summand, we know that $\eta_{K}(C_{\epsilon}^{\complement})$ is less than $\epsilon$ with probability at least $1-\epsilon$ due to Lemma 9.1 for sufficiently large $K$. Collecting all these results, we conclude that

[TABLE]

This completes the proof of the theorem.

9.1 Proof of Theorem 6.5

In order to prove the result, let us first consider the function $f:\mathcal{X}\rightarrow[0,\infty)$, which always takes non-negative values. Let $f_{m}(x):=\min\{f(x),m\}$ be the clipped function, in which case, $f_{m}\in C_{b}(\mathcal{X})$. Define $\Omega_{0}=\{\omega:C(\omega)<\infty\}$, and note that by assumption, $\mathds{P}\left\{\Omega_{0}\right\}=1$. It is obvious that for any $x\in\mathcal{X}$, if we pick $m\geq f(x)$, then $f_{m}(x)=f(x)$.

Note that for this case, for any $m\in\mathbb{N}$, the following inequality holds:

[TABLE]

In the next lemma, we show that each of the summand on the right is small with high probability.

Lemma 9.7**.**

*Let $f:\mathcal{X}\rightarrow[0,\infty)$. If Assumptions 2, 6, and 6 hold, then (14) holds. *

Proof 9.8**.**

For any $\epsilon>0$, there exists (a random natural number) $M_{\epsilon}(\omega)$ such that

[TABLE]

Indeed, one can take $M_{\epsilon}(\omega)=\lceil C(\omega)\rceil$, in which case,

[TABLE]

Further, due to the monotone convergence theorem, there exists $\bar{M}$ such that

[TABLE]

Now, pick $\epsilon,\delta>0$. Let us define $\Omega_{m,K}\subset\Omega$ as the set such that

[TABLE]

For every $m$, pick $K_{m}$ such that $\mathds{P}\left\{\Omega_{m,K_{m}}\right\}<\frac{\delta}{2^{m}}$; such a $K_{m}$ always exists due to Theorem 6.1. Define $\check{\Omega}=\Omega_{0}\bigcap\bigcap_{m=1}^{\infty}\Omega_{m,K_{m}}$. This immediately implies that

[TABLE]

For every $\omega\in\check{\Omega}$, pick $m$ sufficiently large such that $m\geq M_{\epsilon}(\omega)$, $m\geq\bar{M}$, and pick any $K=K_{m}$. We use (20) to conclude that

[TABLE]

*Since $\mathds{P}\left\{\omega\in\check{\Omega}\right\}>1-\delta$, the proof of the lemma is complete. *

We can now complete the proof of Theorem 6.5 using the result above. Pick $f\in C(\mathcal{X})$ and define $f_{+},f_{-}:\mathcal{X}\rightarrow[0,\infty)$ such that $f_{+}=\max\{f,0\}$ and $f_{-}=\max\{-f,0\}$. This immediately yields $f=f_{+}-f_{-}$. Further, we have

[TABLE]

Now, it is easy to conclude from Lemma 9.7 that the convergence in probability result in (14) holds for both $f_{+}$ and $f_{-}$, due to which the convergence in probability result holds for the function $f$ itself.

10 Conclusion

In this paper, we studied the convergence of random sequences generated from certain RSAs used in machine and reinforcement learning problems. If the randomization device used within the algorithm is independent at every iteration, and the maps do not change (for instance, the stepsize is taken as constant), then the random sequence generated can be viewed using the lens of Markov chains. We leveraged the theory of Feller Markov chains to deduce many interesting characteristics of the random sequence and their distributions. Specifically, under reasonable conditions, we showed that the entire random sequence is close to the sequence generated by the exact algorithm with high probability for $n$ sufficiently large. We further showed that the average of the random sequence converges to the mean of the invariant distribution if the sequence if there exists a unique invariant measure of the Markov chain.

We expect that the results presented here can be applied to MDPs over continuous state and action spaces (referred to as continuous MDPs). Indeed, finite time guarantees of empirical value iteration for continuous MDPs with function approximator have been presented in [57, 72, 71, 38, 70, 39, 33] under a variety of assumptions on the MDPs and performance criteria. Convergence of asynchronous algorithm for continuous MDPs with non-parametric function approximation is presented in [69]. It will be interesting to investigate if the output of these algorithms satisfy the sufficient conditions for Theorem 4.1 and 6.1. It will also be interesting to apply the results presented here to variance reduced algorithms [73, 83, 66, 85, 84] that have been developed recently.

Another problem left for future research is to determine bounds on $\mathds{P}\left\{\rho(\hat{z}^{n}_{k},y_{k})\geq\epsilon\right\}$ for any $\epsilon>0$ for every $n$ and $k$. Such bounds would unify our understanding of finite time guarantees for RSAs and allow us to improve the existing algorithms. We hope that the unified framework developed in this paper will be useful for analyzing many other learning algorithms in the future, particularly for analyzing MDPs over compact uncountable state and action spaces.

Appendix A Proof of Lemma 8.4

Since $\mu_{n}\circ\Pi_{l}^{-1}\rightarrow\psi\circ\Pi_{l}^{-1}$ in weak topology as $n\rightarrow\infty$, we conclude that \Big{(}\mu_{n}\circ\Pi_{l}^{-1}\Big{)}_{n\in\mathbb{N}} is tight. For a fixed $\epsilon>0$, let $F_{2}\subset\mathcal{X}$ be the compact set such that

[TABLE]

We now need the following result.

Lemma A.1**.**

If Assumption 4 holds, then for any $g\in U_{b}(\mathcal{X})$, compact set $K\subset\mathcal{X}$ and $\epsilon>0$, there exists $N\in\mathbb{N}$ such that

[TABLE]

Proof A.2**.**

Since $g$ is uniformly continuous, for every $\epsilon>0$, there exists a $\delta_{\epsilon}>0$ such that for any $x,x^{\prime}\in\mathcal{X}$ with $\rho(x,x^{\prime})<\delta_{\epsilon}$, we have $|g(x)-g(x^{\prime})|<\epsilon$. Since Assumption 4 holds, there exists $N_{\epsilon}(g,K)$ such that

[TABLE]

This implies

[TABLE]

*The proof of the lemma is complete. *

We are now in a position to prove the result. Consider the following expressions

[TABLE]

The proof of the lemma is complete.

Appendix B Proof of Lemma 9.5

Define $g_{p}(\hat{z}^{n}_{-i}):=0$ for all $i\in\mathbb{N}$ and $p\in\{0,1,\ldots\}$. Let us define random variables $u_{p,l}$ for $l,p=0,1,\ldots$ as

[TABLE]

Thus, $u_{p,0}=g_{p}(\hat{z}^{n}_{0})$. By definition of $u_{p,l}$, we have for any $l,k\geq 0$

[TABLE]

We use the above expression to establish the following identity.

Lemma B.1**.**

For $M\leq N$, we have

[TABLE]

Proof B.2**.**

For any $l\in\{0,1,\ldots\}$, we have

[TABLE]

This yields

[TABLE]

For any $l,k\geq 0$, (22) yields

[TABLE]

Consider the second summand in (24). For any $k\in\mathbb{N}$, we have

[TABLE]

These expressions immediately establish the equality in (23).

Lemma B.3**.**

For a fixed $p\in\{0,1,\ldots\}$, we have

[TABLE]

Proof B.4**.**

For any $p\in\{0,1,\ldots\}$ and $l\in\mathbb{N}$, let $\mathcal{F}_{p,l-1}$ denote the $\sigma$-algebra generated by $(\hat{z}^{n}_{0},\ldots,\hat{z}^{n}_{l-1},u_{p,0},\ldots,u_{p,l-1})$. We show that $\{u_{p,l}\}_{l=1}^{\infty}$ forms a martingale with respect to the $\sigma$-algebra $\mathcal{F}_{p,l-1}$. By the definition of $u_{p,l}$, we immediately have for any $l\geq 1$,

[TABLE]

This implies

[TABLE]

*The proof then is an immediate consequence of the strong law of large numbers for martingales in [53, p. 66]. *

Now note that for any $M\in\mathbb{N}$, as $N\rightarrow\infty$, the right side of (23) converges to 0 almost surely by Lemma B.3. This yields the result.