\nameJonathan Baxter \emailjbaxter@whizbang.com
\addr4616 Henry Street Pittsburgh, PA 15213
\namePeter L. Bartlett \emailbartlett@barnhilltechnologies.com
###### Abstract

Gradient-based approaches to direct policy search in reinforcement learning have received much recent attention as a means to solve problems of partial observability and to avoid some of the problems associated with policy degradation in value-function methods. In this paper we introduce , a simulation-based algorithm for generating a biased estimate of the gradient of the average reward in Partially Observable Markov Decision Processes (s) controlled by parameterized stochastic policies. A similar algorithm was proposed by \citeAkimura95. The algorithm’s chief advantages are that it requires storage of only twice the number of policy parameters, uses one free parameter (which has a natural interpretation in terms of bias-variance trade-off), and requires no knowledge of the underlying state. We prove convergence of , and show how the correct choice of the parameter is related to the mixing time of the controlled . We briefly describe extensions of  to controlled Markov chains, continuous state, observation and control spaces, multiple-agents, higher-order derivatives, and a version for training stochastic policies with internal states. In a companion paper [6] we show how the gradient estimates generated by  can be used in both a traditional stochastic gradient algorithm and a conjugate-gradient procedure to find local optima of the average reward.

## 1 Introduction

Dynamic Programming is the method of choice for solving problems of decision making under uncertainty [9]. However, the application of Dynamic Programming becomes problematic in large or infinite state-spaces, in situations where the system dynamics are unknown, or when the state is only partially observed. In such cases one looks for approximate techniques that rely on simulation, rather than an explicit model, and parametric representations of either the value-function or the policy, rather than exact representations.

Simulation-based methods that rely on a parametric form of the value function tend to go by the name “Reinforcement Learning,” and have been extensively studied in the Machine Learning literature [8, 25]. This approach has yielded some remarkable empirical successes in a number of different domains, including learning to play checkers [20], backgammon [27, 28], and chess [7], job-shop scheduling [30] and dynamic channel allocation [22].

Despite this success, most algorithms for training approximate value functions suffer from the same theoretical flaw: the performance of the greedy policy derived from the approximate value-function is not guaranteed to improve on each iteration, and in fact can be worse than the old policy by an amount equal to the maximum approximation error over all states. This can happen even when the parametric class contains a value function whose corresponding greedy policy is optimal. We illustrate this with a concrete and very simple example in Appendix A.

An alternative approach that circumvents this problem—the approach we pursue here—is to consider a class of stochastic policies parameterized by , compute the gradient with respect to of the average reward, and then improve the policy by adjusting the parameters in the gradient direction. Note that the policy could be directly parameterized, or it could be generated indirectly from a value function. In the latter case the value-function parameters are the parameters of the policy, but instead of being adjusted to minimize error between the approximate and true value function, the parameters are adjusted to directly improve the performance of the policy generated by the value function.

These “policy-gradient” algorithms have a long history in Operations Research, Statistics, Control Theory, Discrete Event Systems and Machine Learning. Before describing the contribution of the present paper, it seems appropriate to introduce some background material explaining this approach. Readers already familiar with this material may want to skip directly to section 1.2, where the contributions of the present paper are described.

### 1.1 A Brief History of Policy-Gradient Algorithms

For large-scale problems or problems where the system dynamics are unknown, the performance gradient will not be computable in closed form111See equation (17) for a closed-form expression for the performance gradient.. Thus the challenging aspect of the policy-gradient approach is to find an algorithm for estimating the gradient via simulation. Naively, the gradient can be calculated numerically by adjusting each parameter in turn and estimating the effect on performance via simulation (the so-called crude Monte-Carlo technique), but that will be prohibitively inefficient for most problems. Somewhat surprisingly, under mild regularity conditions, it turns out that the full gradient can be estimated from a single simulation of the system. The technique is called the score function or likelihood ratio method and appears to have been first proposed in the sixties [2, 17] for computing performance gradients in i.i.d. (independently and identically distributed) processes.

Specifically, suppose is a performance function that depends on some random variable , and is the probability that , parameterized by . Under mild regularity conditions, the gradient with respect to of the expected performance,

 η(θ)=\boldmathEr(X), (1)

may be written

 ∇η(θ)=\boldmathEr(X)∇q(θ,X)q(θ,X). (2)

To see this, rewrite (1) as a sum

 η(θ)=∑xr(x)q(θ,x),

differentiate (one source of the requirement of “mild regularity conditions”) to obtain

 ∇η(θ)=∑xr(x)∇q(θ,x),

rewrite as

 ∇η(θ)=∑xr(x)∇q(θ,x)q(θ,x)q(θ,x),

and observe that this formula is equivalent to (2).

If a simulator is available to generate samples distributed according to , then any sequence generated i.i.d. according to gives an unbiased estimate,

 ^∇η(θ)=1NN∑i=1r(Xi)∇q(θ,Xi)q(θ,Xi), (3)

of . By the law of large numbers, with probability one. The quantity is known as the likelihood ratio or score function in classical statistics. If the performance function also depends on , then is replaced by in (2).

#### 1.1.1 Unbiased Estimates of the Performance Gradient for Regenerative Processes

Extensions of the likelihood-ratio method to regenerative processes (including Markov Decision Processes or s) were given by \citeAglynn86,glynn90,glynn95 and \citeAreiman86,reiman89, and independently for episodic Partially Observable Markov Decision Processes (s) by \citeAwilliams92, who introduced the  algorithm222A thresholded version of these algorithms for neuron-like elements was described earlier in \citeABarSutAnd83.. Here the i.i.d. samples of the previous section are sequences of states (of random length) encountered between visits to some designated recurrent state , or sequences of states from some start state to a goal state. In this case can be written as a sum

 ∇q(θ,X)q(θ,X)=T−1∑t=0∇pXtXt+1(θ)pXtXt+1(θ), (4)

where is the transition probability from to given parameters . Equation (4) admits a recursive computation over the course of a regenerative cycle of the form , and after each state transition ,

 zt+1=zt+∇pXtXt+1(θ)pXtXt+1(θ), (5)

so that each term in the estimate (3) is of the form333The vector is known in reinforcement learning as an eligibility trace. This terminology is used in \citeABarSutAnd83. . If, in addition, can be recursively computed by

 r(X0,…,Xt+1)=ϕ(r(X0,…,Xt),Xt+1)

for some function , then the estimate for each cycle can be computed using storage of only parameters ( for and parameter to update the performance function ). Hence, the entire estimate (3) can be computed with storage of only real parameters, as follows.

##### Algorithm 1.1: Policy-Gradient Algorithm for Regenerative Processes.
1. Set , , , and ().

2. For each state transition :

• If the episode is finished (that is, , set
,
,
,
.

• Otherwise, set

.

3. If return , otherwise goto 2.

Examples of recursive performance functions include the sum of a scalar reward over a cycle, where is a scalar reward associated with state (this corresponds to being the average reward multiplied by the expected recurrence time ); the negative length of the cycle (which can be implemented by assigning a reward of to each state, and is used when the task is to mimimize time taken to get to a goal state, since in this case is just ); the discounted reward from the start state, , where is the discount factor, and so on.

As \citeAwilliams92 pointed out, a further simplification is possible in the case that is a sum of scalar rewards depending on the state and possibly the time since the starting state (such as , or as above). In that case, the update from a single regenerative cycle may be written as

 Δ=T−1∑t=0∇pXtXt+1(θ)pXtXt+1(θ)[t∑s=0r(Xs,s)+T∑s=t+1r(Xs,s)].

Because changes in have no influence on the rewards associated with earlier states (), we should be able to drop the first term in the parentheses on the right-hand-side and write

 Δ=T−1∑t=0∇pXtXt+1(θ)pXtXt+1(θ)T∑s=t+1r(Xs,s). (6)

Although the proof is not entirely trivial, this intuition can indeed be shown to be correct.

Equation (6) allows an even simpler recursive formula for estimating the performance gradient. Set , and introduce a new variable . As before, set and if , or and otherwise. But now, on each iteration, set . Then is our estimate of . Since is updated on every iteration, this suggests that we can do away with altogether and simply update directly: , where the are suitable step-sizes444The usual requirements on for convergence of a stochastic gradient algorithm are , , and .. Proving convergence of such an algorithm is not as straightforward as normal stochastic gradient algorithms because the updates are not in the gradient direction (in expectation), although the sum of these updates over a regenerative cycle are. \citeAmarbach98 provide the only convergence proof that we know of, albeit for a slightly different update of the form , where is a moving estimate of the expected performance, and is also updated on-line (this update was first suggested in the context of s by \shortciteAjaakola95).

\citeA

marbach98 also considered the case of -dependent rewards (recall the discussion after (3)), as did \citeAbaird98 with their “” algorithm (Value And Policy Search). This last paper contains an interesting insight: through suitable choices of the performance function , one can combine policy-gradient search with approximate value function methods. The resulting algorithms can be viewed as actor-critic techniques in the spirit of \citeABarSutAnd83; the policy is the actor and the value function is the critic. The primary motivation is to reduce variance in the policy-gradient estimates. Experimental evidence for this phenomenon has been presented by a number of authors, including \citeABarSutAnd83, \citeAkimura98a, and \citeAbaird98. More recent work on this subject includes that of \shortciteAsutton99 and \shortciteAkonda99. We discuss the use of -style updates further in Section 6.2.

So far we have not addressed the question of how the parameterized state-transition probabilities arise. Of course, they could simply be generated by parameterizing the matrix of transition probabilities directly. Alternatively, in the case of s or s, state transitions are typically generated by feeding an observation that depends stochastically on the state into a parameterized stochastic policy, which selects a control at random from a set of available controls (approximate value-function based approaches that generate controls stochastically via some form of lookahead also fall into this category). The distribution over successor states is then a fixed function of the control. If we denote the probability of control given parameters and observation by , then all of the above discussion carries through with replaced by . In that case, Algorithm 1.1 is precisely Williams’  algorithm.

Algorithm 1.1 and the variants above have been extended to cover multiple agents \shortcitepeshkin00, policies with internal state \shortcitemeuleau99, and importance sampling methods \shortcitemeuleau00. We also refer the reader to the work of \citeArubinstein93 and \citeArubinstein98 for in-depth analysis of the application of the likelihood-ratio method to Discrete-Event Systems (), in particular networks of queues. Also worth mentioning is the large literature on Infinitesimal Perturbation Analysis (IPA), which seeks a similar goal of estimating performance gradients, but operates under more restrictive assumptions than the likelihood-ratio approach; see, for example, \citeAho91.

#### 1.1.2 Biased Estimates of the Performance Gradient

All the algorithms described in the previous section rely on an identifiable recurrent state , either to update the gradient estimate, or in the case of the on-line algorithm, to zero the eligibility trace . This reliance on a recurrent state can be problematic for two main reasons:

1. The variance of the algorithms is related to the recurrence time between visits to , which will typically grow as the state space grows. Furthermore, the time between visits depends on the parameters of the policy, and states that are frequently visited for the initial value of the parameters may become very rare as performance improves.

2. In situations of partial observability it may be difficult to estimate the underlying states, and therefore to determine when the gradient estimate should be updated, or the eligibility trace zeroed.

If the system is available only through simulation, it seems difficult (if not impossible) to obtain unbiased estimates of the gradient direction without access to a recurrent state. Thus, to solve 1 and 2, we must look to biased estimates. Two principle techniques for introducing bias have been proposed, both of which may be viewed as artificial truncations of the eligibility trace . The first method takes as a starting point the formula555For ease of exposition, we have kept the expression for in terms of the likelihood ratios which rely on the availability of the underlying state . If is not available, should be replaced with . for the eligibility trace at time :

 zt=t−1∑s=0∇pXsXs+1(θ)pXsXs+1(θ)

and simply truncates it at some (fixed, not random) number of terms looking backwards [13, 18, 19, 11]:

 zt(n):=t−1∑s=t−n∇pXsXs+1(θ)pXsXs+1(θ). (7)

The eligibility trace is then updated after each transition by

 zt+1(n)=zt(n)+∇pXtXt+1(θ)pXtXt+1(θ)−∇pXt−nXt−n+1(θ)pXt−nXt−n+1(θ), (8)

and in the case of state-based rewards , the estimated gradient direction after steps is

 ^∇nη(θ):=1T−n+1T∑t=nzt(n)r(Xt). (9)

Unless exceeds the maximum recurrence time (which is infinite in an ergodic Markov chain), is a biased estimate of the gradient direction, although as , the bias approaches zero. However the variance of diverges in the limit of large . This illustrates a natural trade-off in the selection of the parameter : it should be large enough to ensure the bias is acceptable (the expectation of should at least be within of the true gradient direction), but not so large that the variance is prohibitive. Experimental results by \citeAcao98 illustrate nicely this bias/variance trade-off.

One potential difficulty with this method is that the likelihood ratios must be remembered for the previous time steps, requiring storage of parameters. Thus, to obtain small bias, the memory may have to grow without bound. An alternative approach that requires a fixed amount of memory is to discount the eligibility trace, rather than truncating it:

 zt+1(β):=βzt(β)+∇pXtXt+1(θ)pXtXt+1(θ), (10)

where and is a discount factor. In this case the estimated gradient direction after steps is simply

 ^∇βη(θ):=1TT−1∑t=0r(Xt)zt(β). (11)

This is precisely the estimate we analyze in the present paper. A similar estimate with replaced by where is a reward baseline was proposed by \shortciteAkimura95,kimura97 and for continuous control by \citeAkimura98. In fact the use of in place of does not affect the expectation of the estimates of the algorithm (although judicious choice of the reward baseline can reduce the variance of the estimates). While the algorithm presented by \citeAkimura95 provides estimates of the expectation under the stationary distribution of the gradient of the discounted reward, we will show that these are in fact biased estimates of the gradient of the expected discounted reward. This arises because the stationary distribution itself depends on the parameters. A similar estimate to (11) was also proposed by \citeAmarbach98, but this time with replaced by , where is an estimate of the average reward, and with zeroed on visits to an identifiable recurrent state.

As a final note, observe that the eligibility traces and defined by (10) and (8) are simply filtered versions of the sequence , a first-order, infinite impulse response filter in the case of and an -th order, finite impulse response filter in the case of . This raises the question, not addressed in this paper, of whether there is an interesting theory of optimal filtering for policy-gradient estimators.

### 1.2 Our Contribution

We describe , a general algorithm based upon (11) for generating a biased estimate of the performance gradient in general s controlled by parameterized stochastic policies. Here denotes the average reward of the policy with parameters .  does not rely on access to an underlying recurrent state. Writing for the expectation of the estimate produced by , we show that , and more quantitatively that is close to the true gradient provided exceeds the mixing time of the Markov chain induced by the 666The mixing-time result in this paper applies only to Markov chains with distinct eigenvalues. Better estimates of the bias and variance of  may be found in \citeAjcss_01, for more general Markov chains than those treated here, and for more refined notions of the mixing time. Roughly speaking, the variance of  grows with , while the bias decreases as a function of .. As with the truncated estimate above, the trade-off preventing the setting of arbitrarily close to is that the variance of the algorithm’s estimates increase as approaches . We prove convergence with probability 1 of  for both discrete and continuous observation and control spaces. We present algorithms for both general parameterized Markov chains and s controlled by parameterized stochastic policies.

There are several extensions to  that we have investigated since the first version of this paper was written. We outline these developments briefly in Section 7.

In a companion paper we show how the gradient estimates produced by  can be used to perform gradient ascent on the average reward [6]. We describe both traditional stochastic gradient algorithms, and a conjugate-gradient algorithm that utilizes gradient estimates in a novel way to perform line searches. Experimental results are presented illustrating both the theoretical results of the present paper on a toy problem, and practical aspects of the algorithms on a number of more realistic problems.

## 2 The Reinforcement Learning Problem

We model reinforcement learning as a Markov decision process () with a finite state space , and a stochastic matrix777A stochastic matrix has for all and for all . giving the probability of transition from state to state . Each state has an associated reward888All the results in the present paper apply to bounded stochastic rewards, in which case is the expectation of the reward in state . . The matrix belongs to a parameterized class of stochastic matrices, . Denote the Markov chain corresponding to by . We assume that these Markov chains and rewards satisfy the following assumptions:

###### Assumption 1.

Each has a unique stationary distribution satisfying the balance equations

 π′(θ)P(θ)=π′(θ) (12)

(throughout denotes the transpose of ).

###### Assumption 2.

The magnitudes of the rewards, , are uniformly bounded by for all states .

Assumption 1 ensures that the Markov chain forms a single recurrent class for all parameters . Since any finite-state Markov chain always ends up in a recurrent class, and it is the properties of this class that determine the long-term average reward, this assumption is mainly for convenience so that we do not have to include the recurrence class as a quantifier in our theorems. However, when we consider gradient-ascent algorithms \citeAjair_01b, this assumption becomes more restrictive since it guarantees that the recurrence class cannot change as the parameters are adjusted.

Ordinarily, a discussion of s would not be complete without some mention of the actions available in each state and the space of policies available to the learner. In particular, the parameters would usually determine a policy (either directly or indirectly via a value function), which would then determine the transition probabilities . However, for our purposes we do not care how the dependence of on arises, just that it satisfies Assumption 1 (and some differentiability assumptions that we shall meet in the next section). Note also that it is easy to extend this setup to the case where the rewards also depend on the parameters or on the transitions . It is equally straightforward to extend our algorithms and results to these cases. See Section 6.1 for an illustration.

The goal is to find a maximizing the average reward:

 η(θ):=limT→∞\boldmathEθ[1TT−1∑t=0r(Xt)∣∣ ∣∣X0=i],

where denotes the expectation over all sequences with transitions generated according to . Under Assumption 1, is independent of the starting state and is equal to

 η(θ)=n∑i=1π(θ,i)r(i)=π′(θ)r, (13)

where [9].

## 3 Computing the Gradient of the Average Reward

For general s little will be known about the average reward , hence finding its optimum will be problematic. However, in this section we will see that under general assumptions the gradient exists, and so local optimization of is possible.

To ensure the existence of suitable gradients (and the boundedness of certain random variables), we require that the parameterized class of stochastic matrices satisfies the following additional assumption.

###### Assumption 3.

The derivatives,

 ∇P(θ):=[∂pij(θ)∂θk]i,j=1…n;k=1…K

exist for all . The ratios

 ⎡⎢ ⎢ ⎢⎣∣∣∣∂pij(θ)∂θk∣∣∣pij(θ)⎤⎥ ⎥ ⎥⎦i,j=1…n;k=1…K

are uniformly bounded by for all .

The second part of this assumption allows zero-probability transitions only if is also zero, in which case we set . One example is if is a forbidden transition, so that for all . Another example satisfying the assumption is

 pij(θ)=eθij∑nj=1eθij,

where are the parameters of , for then

 ∂pij(θ)/∂θijpij(θ) =1−pij(θ),and ∂pij(θ)/∂θklpij(θ) =−pkl(θ).

Assuming for the moment that exists (this will be justified shortly), then, suppressing dependencies,

 ∇η=∇π′r, (14)

since the reward does not depend on . Note that our convention for in this paper is that it takes precedence over all other operations, so . Equations like (14) should be regarded as shorthand notation for equations of the form

 ∂η(θ)∂θk=[∂π(θ,1)∂θk,…,∂π(θ,n)∂θk][r(1),…,r(n)]′

where . To compute , first differentiate the balance equations (12) to obtain

 ∇π′P+π′∇P=∇π′,

and hence

 ∇π′(I−P)=π′∇P. (15)

The system of equations defined by (15) is under-constrained because is not invertible (the balance equations show that has a left eigenvector with zero eigenvalue). However, let denote the -dimensional column vector consisting of all s, so that is the matrix with the stationary distribution in each row. Since , we can rewrite (15) as

 ∇π′[I−(P−eπ′)]=π′∇P.

To see that the inverse exists, let be any matrix satisfying . Then we can write

 limT→∞[(I−A)T∑t=0At] =limT→∞[T∑t=0At−T+1∑t=1At] =I−limT→∞AT+1 =I.

Thus,

 (I−A)−1=∞∑t=0At.

It is easy to prove by induction that which converges to as by Assumption 1. So exists and is equal to . Hence, we can write

 ∇π′=π′∇P[I−P+eπ′]−1, (16)

and so999The argument leading to (16) coupled with the fact that is the unique solution to (12) can be used to justify the existence of . Specifically, we can run through the same steps computing the value of for small and show that the expression (16) for is the unique matrix satisfying .

 ∇η=π′∇P[I−P+eπ′]−1r. (17)

For s with a sufficiently small number of states, (17) could be solved exactly to yield the precise gradient direction. However, in general, if the state space is small enough that an exact solution of (17) is possible, then it will be small enough to derive the optimal policy using policy iteration and table-lookup, and there would be no point in pursuing a gradient based approach in the first place101010Equation (17) may still be useful for s, since in that case there is no tractable dynamic programming algorithm..

Thus, for problems of practical interest, (17) will be intractable and we will need to find some other way of computing the gradient. One approximate technique for doing this is presented in the next section.

## 4 Approximating the Gradient in Parameterized Markov Chains

In this section, we show that the gradient can be split into two components, one of which becomes negligible as a discount factor approaches .

For all , let denote the vector of expected discounted rewards from each state :

 Jβ(θ,i):=\boldmathEθ[∞∑t=0βtr(Xt)∣∣ ∣∣X0=i]. (18)

Where the dependence is obvious, we just write .

###### Proposition 1.

For all and ,

 ∇η=(1−β)∇π′Jβ+βπ′∇PJβ. (19)
###### Proof.

Observe that satisfies the Bellman equations:

 Jβ=r+βPJβ. (20)

[9]. Hence,

 ∇η =∇π′r =∇π′[Jβ−βPJβ] =∇π′Jβ−β∇π′Jβ+βπ′∇PJβ by (15) =(1−β)∇π′Jβ+βπ′∇PJβ.

We shall see in the next section that the second term in (19) can be estimated from a single sample path of the Markov chain. In fact, Theorem 1 in [14] shows that the gradient estimates of the algorithm presented in that paper converge to . By the Bellman equations (20), this is equal to , which implies . Thus the algorithm of \citeAkimura97 also estimates the second term in the expression for given by (19). It is important to note that —the two quantities disagree by the first term in (19). This arises because the the stationary distribution itself depends on the parameters. Hence, the algorithm of \citeAkimura97 does not estimate the gradient of the expected discounted reward. In fact, the expected discounted reward is simply times the average reward \shortcite[Fact 7]singh94a, so the gradient of the expected discounted reward is proportional to the gradient of the average reward.

The following theorem shows that the first term in (19) becomes negligible as approaches . Notice that this is not immediate from Proposition 1, since can become arbitrarily large in the limit .

###### Theorem 2.

For all ,

 ∇η=limβ→1∇βη, (21)

where

 ∇βη:=π′∇PJβ. (22)
###### Proof.

Recalling equation (17) and the discussion preceeding it, we have111111Since , (23) motivates a different kind of algorithm for estimating based on differential rewards [16].

 ∇η=π′∇P∞∑t=0[Pt−eπ′]r. (23)

But since is a stochastic matrix, so (23) can be rewritten as

 ∇η=π′[∞∑t=0∇PPt]r. (24)

Now let be a discount factor and consider the expression

 f(β):=π′[∞∑t=0∇P(βP)t]r (25)

Clearly . To complete the proof we just need to show that .

Since , we can invoke the observation before (16) to write

 ∞∑t=0(βP)t=[I−βP]−1.

In particular, converges, so we can take back out of the sum in the right-hand-side of (25) and write121212We cannot back out of the sum in the right-hand-side of (24) because diverges (). The reason converges is that becomes orthogonal to in the limit of large . Thus, we can view as a sum of two orthogonal components: an infinite one in the direction and a finite one in the direction . It is the finite component that we need to estimate. Approximating with is a way of rendering the -component finite while hopefully not altering the -component too much. There should be other substitutions that lead to better approximations (in this context, see the final paragraph in Section 1.1).

 f(β)=π′∇P[∞∑t=0βtPt]r. (26)

But . Thus . ∎

Theorem 2 shows that is a good approximation to the gradient as approaches , but it turns out that values of very close to lead to large variance in the estimates of that we describe in the next section. However, the following theorem shows that need not be too small, provided the transition probability matrix has distinct eigenvalues, and the Markov chain has a short mixing time. From any initial state, the distribution over states of a Markov chain converges to the stationary distribution, provided the assumption (Assumption 1) about the existence and uniqueness of the stationary distribution is satisfied [<]see, for example,¿[Theorem 15.8.1, p. 552]lancaster85. The spectral resolution theorem [15, Theorem 9.5.1, p. 314] implies that the distribution converges to stationarity at an exponential rate, and the time constant in this convergence rate (the mixing time) depends on the eigenvalues of the transition probability matrix. The existence of a unique stationary distribution implies that the largest magnitude eigenvalue is and has multiplicity , and the corresponding left eigenvector is the stationary distribution. We sort the eigenvalues in decreasing order of magnitude, so that for some . It turns out that determines the mixing time of the chain.

The following theorem shows that if is small compared to , the gradient approximation described above is accurate. Since we will be using the estimate as a direction in which to update the parameters, the theorem compares the directions of the gradient and its estimate. In this theorem, denotes the spectral condition number of a nonsingular matrix , which is defined as the product of the spectral norms of the matrices and ,

 κ2(A)=∥A∥2∥A−1∥2,

where

 ∥A∥2=maxx:∥x∥=1∥Ax∥,

and denotes the Euclidean norm of the vector .

###### Theorem 3.

Suppose that the transition probability matrix satisfies Assumption 1 with stationary distribution , and has distinct eigenvalues. Let be the matrix of right eigenvectors of corresponding, in order, to the eigenvalues . Then the normalized inner product between and satisfies

 1−∇η⋅β∇βη∥∇η∥2≤κ2(Π1/2S)∥∇(√π1,…,√πn)∥∥∇η∥√r′Πr1−β1−β|λ2|, (27)

where .

Notice that is the expectation under the stationary distribution of .

As well as the mixing time (via ), the bound in the theorem depends on another parameter of the Markov chain: the spectral condition number of . If the Markov chain is reversible (which implies that the eigenvectors are orthogonal), this is equal to the ratio of the maximum to the minimum probability of states under the stationary distribution. However, the eigenvectors do not need to be nearly orthogonal. In fact, the condition that the transition probability matrix have distinct eigenvalues is not necessary; without it, the condition number is replaced by a more complicated expression involving spectral norms of matrices of the form .

###### Proof.

The existence of distinct eigenvalues implies that can be expressed as , where [15, Theorem 4.10.2, p 153]. It follows that for any polynomial , we can write .

Now, Proposition 1 shows that . But

 (1−β)Jβ =(1−β)(r+βPr+β2P2r+⋯) =(1−β)(I+βP+β2P2+⋯)r =(1−β)S(∞∑t=0βtΛt)S−1r =(1−β)n∑j=1xjy′j(∞∑t=0(βλj)t)r,

where .

It is easy to verify that is the left eigenvector corresponding to , and that we can choose and . Thus we can write

 (1−β)Jβ =(1−β)eπ′r+n∑j=2xjy′j(∞∑t=0(1−β)(βλj)t)r =(1−β)eη+n∑j=2xjy′j(1−β1−βλj)r =(1−β)eη+SMS−1r,

where

 M=diag(0,1−β1−βλ2,…,1−β1−βλn).

It follows from this and Proposition 1 that

 1−∇η⋅β∇βη∥∇η∥2 =1−∇η⋅(∇η−∇π′(1−β)Jβ)∥∇η∥2 =∇η⋅∇π′(1−β)Jβ∥∇η∥2 =∇η⋅∇π′((1−β)eη+SMS−1r)∥∇η∥2 =∇η⋅∇π′SMS−1r∥∇η∥2 ≤∥∥∇π′SMS−1r∥∥∥∇η∥,

by the Cauchy-Schwartz inequality. Since , we can apply the Cauchy-Schwartz inequality again to obtain

 1−∇η⋅β∇βη∥∇η∥2≤∥∥∇(√π′)∥∥∥∥Π1/2SMS−1r∥∥∥∇η∥. (28)

We use spectral norms to bound the second factor in the numerator. It is clear from the definition that the spectral norm of a product of nonsingular matrices satisfies , and that the spectral norm of a diagonal matrix is given by . It follows that

 ∥∥Π1/2SMS−1r∥∥ =∥∥Π1/2SMS−1Π−1/2Π1/2r∥∥ ≤∥∥Π1/2S∥∥2∥∥S−1Π−1/2∥∥2∥∥Π1/2r∥∥∥M∥2 ≤κ2(Π1/2S)√r′Πr1−β1−β|λ2|.

Combining with Equation (28) proves (27). ∎

## 5 Estimating the Gradient in Parameterized Markov Chains

Algorithm 1 introduces  (Markov Chain Gradient), an algorithm for estimating the approximate gradient from a single on-line sample path from the Markov chain .  requires only reals to be stored, where is the dimension of the parameter space: parameters for the eligibility trace , and parameters for the gradient estimate . Note that after time steps is the average so far of ,

 ΔT=1TT−1∑t=0ztr(Xt).