Batch Stationary Distribution Estimation

Introduction

Markov chains are a pervasive modeling tool in applied mathematics of particular importance in stochastic modeling and machine learning. A key property of an ergodic Markov chain is the existence of a unique stationary distribution; i.e., the long-run distribution of states that remains invariant under the transition kernel. In this paper, we consider a less well studied but still important version of the stationary distribution estimation problem, where one has access to a set of sampled transitions from a given Markov chain, but does not know the mechanism by which the probe points were chosen, nor is able to gather additional data from the underlying process. Nevertheless, one would still like to estimate target properties of the stationary distribution, such as the expected value of a random variable of interest.

This setting is inspired by many practical scenarios where sampling from the Markov process is costly or unavailable, but data has already been collected and available for analysis. A simple example is a queueing system consisting of a service desk that serves customers in a queue. Queue length changes stochastically as customers arrive or leave after being served. The long-term distribution of queue length (i.e., the stationary distribution of the underlying Markov chain) is the object of central interest for managing such a service (Haviv, 2009; Serfozo, 2009). In practice, however, queue lengths are physical quantities that can only be measured for moderate periods, perhaps on separate occasions, but rarely for sufficient time to ensure the (stochastic) queue length has reached the stationary distribution. Since the measurement process itself is expensive, it is essential to make reasonable inferences about the stationary distribution from the collected data alone.

We investigate methods for estimating properties of the stationary distribution solely from a batch of previously collected data. The key idea is to first estimate a correction ratio function over the given data, which can then be used to estimate expectations of interest with respect to the stationary distribution. To illustrate, consider an ergodic Markov chain with state space $\mathcal{X}$, transition kernel ${\mathcal{T}}$, and a unique stationary distribution $\mu$ that satisfies

Assume we are given a fixed sample of state transitions, $\mathcal{D}=\left\{\left(x,x^{\prime}\right)_{i=1}^{n}\right\}\sim{\mathcal{T}}\left(x^{\prime}|x\right)p\left(x\right)$, such that each $x$ has been sampled according to an unknown probe distribution $p$, but each $x^{\prime}$ has been sampled according to the true underlying transition kernel, $x^{\prime}|x\sim{\mathcal{T}}\left(x^{\prime}|x\right)$. Below we investigate procedures for estimating the point-wise ratios, $\widehat{\tau}\left(x_{i}\right)\approx\frac{\mu\left(x_{i}\right)}{p\left(x_{i}\right)}$, such that the weighted empirical distribution

can be used to approximate $\mu$ directly, or further used to estimate the expected value of some target function(s) of $x$ with respect to $\mu$. Crucially, the approach we propose does not require knowledge of the probe distribution $p$, nor does it require additional access to samples drawn from the transition kernel ${\mathcal{T}}$, yet we will be able to establish consistency of the estimation strategy under general conditions.

In addition to developing the fundamental approach, we demonstrate its applicability and efficacy in a range of important scenarios beyond queueing, including:

Stochastic differential equations (SDEs) SDEs are an essential modeling tool in many fields like statistical physics (Kadanoff, 2000), finance (Oksendal, 2013) and molecular dynamcis (Liu, 2001). An autonomous SDE describes the instantaneous change of a random variable $X$ by

where $f\left(X\right)$ is a drift term, $\sigma\left(X\right)$ a diffusion term, and $W$ the Wiener process. Given data \mathcal{D}=\big{\{}\left(x,x^{\prime}\right)_{i=1}^{n}\big{\}} such that $x\sim p\left(x\right)$ is drawn from an unknown probe distribution and $x^{\prime}$ is the next state after a small time step according to (2), we consider the problem of estimating quantities of the stationary distribution $\mu$ when one exists.

Off-policy evaluation (OPE) Another important application is behavior-agnostic off-policy evaluation (Nachum et al., 2019) in reinforcement learning (RL). Consider a Markov decision process (MDP) specified by $M=\langle\mathcal{S},\mathcal{A},P,R\rangle$, such that $\mathcal{S}$ and $\mathcal{A}$ are the state and action spaces, $P$ is the transition function, and $R$ is the reward function (Puterman, 2014). Given a policy $\pi$ that maps $s\in\mathcal{S}$ to a distribution over $\mathcal{A}$, a random trajectory can be generated starting from an initial state $s_{0}$: $(s_{0},a_{0},r_{0},s_{1},a_{1},r_{1},\ldots)$, where $a_{t}\sim\pi(\cdot|s_{t})$, $s_{t+1}\sim P(\cdot|s_{t},a_{t})$ and $r_{t}\sim R\left(s_{t},a_{t}\right)$. The value of a policy $\pi$ is defined to be its long-term average per-step reward:

where $d_{\pi}$ denotes the limiting distribution over states $\mathcal{S}$ of the Markov process induced by $\pi$. In behavior-agnostic off-policy evaluation, one is given a target policy $\pi$ and a set of transitions $\mathcal{D}=\left\{(s,a,r,s^{\prime})_{i=1}^{n}\right\}\sim P\left(s^{\prime}|s,a\right)p\left(s,a\right)$, potentially generated by multiple behavior policies. From such data, an estimate for $\rho\left(\pi\right)$ can be formed in terms of a stationary ratio estimator:

We refer the interested readers to Section 5.4 and Appendix C for further discussion.

For the remainder of the paper, we will outline four main contributions. First, we generalize the classical power iteration method to obtain an algorithm, the Variational Power Method (VPM), that can work with arbitrary parametrizations in a functional space, allowing for a flexible yet practical approach. Second, we prove the consistency and convergence of VPM. Third, we illustrate how a diverse set of stationary distribution estimation problems, including those above, can be addressed by VPM in a unified manner. Finally, we demonstrate empirically that VPM significantly improves estimation quality in a range of applications, including queueing, sampling, SDEs and OPE.

Variational Power Method

To develop our approach, first recall the definition of ${\mathcal{T}}$ and $\mu$ in (1). We make the following assumption about ${\mathcal{T}}$ and $\mu$ throughout the paper.

The transition operator ${\mathcal{T}}$ has a unique stationary distribution, denoted $\mu$.

Conditions under which this assumption holds are mild, and have been extensively discussed in standard textbooks (Meyn et al., 2009; Levin and Peres, 2017).

Next, to understand the role of the probe distribution $p$, note that we can always rewrite the stationary distribution as $\mu=p\circ\tau$ (i.e., $\mu\left(x\right)\!=\!p\left(x\right)\tau\left(x\right)$, hence $\tau\left(x\right)\!=\!\frac{\mu\left(x\right)}{p\left(x\right)}$), provided the following assumption holds.

The stationary distribution $\mu$ is absolutely continuous w.r.t. $p$. That is, there exists $C<\infty$ such that $\left\|\tau\right\|_{\infty}\leqslant C$.

Assumption 2 follows previous work (Liu and Lee, 2017; Nachum et al., 2019), and is common in density ratio estimation (Sugiyama et al., 2008; Gretton et al., 2009) and off-policy evaluation (Wang et al., 2017; Xie et al., 2019).

Combining these two assumptions, definition (1) yields

This development reveals how, under the two stated assumptions, there is sufficient information to determine the unique ratio function $\tau$ that ensures $p\circ\tau=\mu$ in principle. Given such a function $\tau$, we can then base inferences about $\mu$ solely on data sampled from $p$ and $\tau$.

To develop a practical algorithm for recovering $\tau$ from the constraint (4), in function space, we first consider the classical power method for recovering the $\mu$ that satisfies (1). From (1) it can be seen that the stationary distribution $\mu$ is an eigenfunction of ${\mathcal{T}}$. Moreover, it is the principal eigenfunction, corresponding to the largest eigenvalue $\lambda_{1}=1$. In the simpler case of finite $\mathcal{X}$, the vector $\mu$ is the principal (right) eigenvector of the transposed transition matrix. A standard approach to computing $\mu$ is then the power method:

whose iterates converge to $\mu$ at a rate linear in $\left|\lambda_{2}\right|$, where $\lambda_{2}$ is the second largest eigenvalue of ${\mathcal{T}}$. For ergodic Markov chains, one has $|\lambda_{2}|<1$ (Meyn et al., 2009, Chap 20).

Our initial aim is to extend this power iteration approach to the constraint (4) without restricting the domain $\mathcal{X}$ to be finite. This can be naturally achieved by the update

where the division is element-wise. Clearly the fixed point of (6) corresponds to the solution of (4) under the two assumptions stated above. Furthermore, just as for $\mu_{t}$ in (5), $\tau_{t}$ in (6) also converges to $\tau$ at a linear rate for finite $\mathcal{X}$. Unfortunately, the update (6) cannot be used directly in a practical algorithm for two important reasons. First, we do not have a point-wise evaluator for ${\mathcal{T}}_{p}$, but only samples from ${\mathcal{T}}_{p}$. Second, the operator ${\mathcal{T}}_{p}$ is applied to a function $\tau_{t}$, which typically involves an intractable integral over $\mathcal{X}$ in general. To overcome these issues, we propose a variational method that considers a series of reformulated problems whose optimal solutions correspond to the updates (6).

To begin to develop a practical variational approach, first note that (6) operates directly on the density ratio, which implies the density ratio estimation techniques of Nguyen et al. (2008) and Sugiyama et al. (2012) can be applied. Let $\phi$ be a lower semicontinuous, convex function satisfying $\phi\left(1\right)=0$, and consider the induced $f$-divergence,

To apply this construction to our setting, first consider solving a problem of the following form in the dual space:

where to achieve (9) we have applied the inductive assumption that $\tau_{t}=\partial\phi^{*}(\nu_{t})$. Then, by the optimality property of $\nu_{t+1}$, we know that the solution $\nu_{t+1}$ must satisfy

hence the updated ratio $\tau_{t+1}$ in (6) can be directly recovered from the dual solution $\nu_{t+1}$, while also retaining the inductive property that $\tau_{t+1}=\partial\phi^{*}(\nu_{t+1})$ for the next iteration.

These developments can be further simplified by considering the specific choice $\phi^{*}\left(x\right)=x^{2}/2$, which satisfies $\phi^{*}=\phi$ and simplifies the overall update to

Crucially, this variational update (11) determines the same update as (6), but overcomes the two aforementioned difficulties. First, it bypasses the direct evaluation of ${\mathcal{T}}_{p}$ and $p$, and allows these to be replaced by unbiased estimates of expectations extracted from the data. Second, it similarly bypasses the intractability of the operator application ${\mathcal{T}}_{p}\tau_{t}$ in the functional space, replacing this with an expectation of $\tau_{t}\circ\tau$ that can also be directly estimated from the data.

We now discuss some practical refinements of the approach.

2 Maintaining Normalization

Although it might appear that solving (12) requires one to solve a sequence of regularized problems

with increasing $\lambda\rightarrow\infty$, to ensure the constraint in (12) is satisfied exactly, we note that this additional expense can be entirely avoided for the specific problem we are considering.

3 Avoiding Double Sampling

Crucially, the dual variable $v$ is a scalar, making this problem much simpler than dual embedding (Dai et al., 2017), where the dual variables form a parameterized function that introduces approximation error. The problem (14) is a straightforward convex-concave objective with respect to $\left(\tau,v\right)$ that can be optimized by stochastic gradient descent.

4 Damped Iteration

To control the error due to sampling, we introduce a damped version of the update (Ryu and Boyd, 2016), where instead of performing a stochastic update $\tau_{t+1}={\textstyle\frac{\widehat{{\mathcal{T}}}_{p}}{p}}\tau_{t}$, we instead perform a damped update given by

where $\alpha_{t}\in(0,1)$ is a stepsize parameter. Intuitively, the update error introduced by the stochasticity of $\widehat{{\mathcal{T}}}_{p}$ is now controlled by the stepsize $\alpha_{t}$. The choice of stepsize and convergence of the algorithm is discussed in Section 3.

The damped iteration can be conveniently implemented with minor modifications to the previous objective. We only need to change the sample from ${\mathcal{T}}_{p}$ in (14) by a weighted sample:

5 A Practical Algorithm

After convergence of $\tau_{\theta}$ in each iteration, the reference network is updated by setting $\tau_{t+1}=\tau_{\theta}$. Note that one may apply other gradient-based optimizers instead of SGD.

Convergence Analysis

We now demonstrate that the final algorithm obtains sufficient control over error accumulation to achieve consistency. For notation brevity, we discuss the result for the simpler form (5) instead of the ratio form (6). The argument easily extends to the ratio form.

Starting from the plain stochastic update $\mu_{t}=\widehat{{\mathcal{T}}}\mu_{t-1}$, the damped update can be expressed by

where $\epsilon$ is the error due to stochasticity in $\widehat{{\mathcal{T}}}$. The following theorem establishes the convergence properties of the damped iteration.

Under mild conditions, after $t$ iteration with step-size $\alpha_{t}=1/\sqrt{t}$, we have

The precise version of the theorem statement, together with a complete proof, is given in Appendix B.

Note that the optimization quality depends on the number of samples, the approximation error of the parametric family, and the optimization algorithm. There is a complex trade-off between these factors (Bottou and Bousquet, 2008). On one hand, with more data, the statistical error is reduced, but the computational cost of the optimization increases. On the other hand, with a more flexible parametrization, such as neural networks, reduces the approximation error, but adds to the difficulty of optimization as the problem might no longer be convex. Alternatively, if the complexity of the parameterized family is increased, the consequences of statistical error also increases.

Representing $\tau$ in a reproducing kernel Hilbert space (RKHS) is a particularly interesting case, because the problem (14) becomes convex, hence the optimization error of the empirical surrogate is reduced to zero. Nguyen et al. (2008, Theorem 2) show that, under mild conditions, the statistical error can be bounded in rate $\mathcal{O}\left(n^{-\frac{1}{2+\beta}}\right)$ in terms of Hellinger distance ($\beta$ denotes the exponent in the bracket entropy of the RKHS), while the approximation error will depend on the RKHS (Bach, 2014).

Related Work

The algorithm we have developed reduces distribution estimation to density ratio estimation, which has been extensively studied in numerous contexts. One example is learning under covariate shift (Shimodaira, 2000), where the ratio $\tau$ can be estimated by different techniques (Gretton et al., 2009; Nguyen et al., 2008; Sugiyama et al., 2008; Sugiyama and Kawanabe, 2012). These previous works differ from the current setting in that they require data to be sampled from both the target and proposal distributions. By contrast, we consider a substantially more challenging problem, where only data sampled from the proposal is available, and the target distribution is given only implicitly by (1) through the transition kernel ${\mathcal{T}}$. A more relevant approach is Stein importance sampling (Liu and Lee, 2017), where the ratio is estimated by minimizing the kernelized Stein discrepancy (Liu et al., 2016). However, it requires additional gradient information about the target potential, whereas our method only requires sampled transitions. Moreover, the method of Liu and Lee (2017) is computationally expensive and does not extrapolate to new examples.

The algorithm we develop in this paper is inspired by the classic power method for finding principal eigenvectors. Many existing works have focused on the finite-dimension setting (Balsubramani et al., 2013; Hardt and Price, 2014; Yang et al., 2017), while Kim et al. (2005) and Xie et al. (2015) have extended the power method to the infinite-dimension case using RKHS. Not only do these algorithms require access to the transition kernel ${\mathcal{T}}$, but they also require tractable operator multiplications. In contrast, our method avoids direct interaction with the operator ${\mathcal{T}}$, and can use flexible parametrizations (such as neural networks) to learn the density ratio without per-step renormalization.

Another important class of methods for estimating or sampling from stationary distributions are based on simulations. A prominent example is Markov chain Monte Carlo (MCMC), which is widely used in many statistical inference scenarios (Andrieu et al., 2003; Koller and Friedman, 2009; Welling and Teh, 2011). Existing MCMC methods (e.g., Neal et al., 2011; Hoffman and Gelman, 2014) require repeated, and often many, interactions with the transition operator ${\mathcal{T}}$ to acquire a single sample from the stationary distribution. Instead, VPM can be applied when only a fixed sample is available. Interestingly, this suggests that VPM can be used to “post-process” samples generated from typical MCMC methods to possibly make more effective use of the data. We demonstrated this possibility empirically in Section 5. Unlike VPM, other post-processing methods (Oates et al., 2017) require additional information about the target distribution (Robert and Casella, 2004). Recent advances have also shown that learning parametric samplers can be beneficial (Song et al., 2017; Li et al., 2019), but require the potential function. In contrast, VPM directly learns the stationary density ratio solely from transition data.

One important application of VPM is off-policy RL (Precup et al., 2001). In particular, in off-policy evaluation (OPE), one aims to evaluate a target policy’s performance, given data collected from a different behavior policy. This problem matches our proposed framework as the collected data naturally consists of transitions from a Markov chain, and one is interested in estimating quantities computed from the stationary distribution of a different policy. (See Appendix C for a detailed description of how the VPM algorithm can be applied to OPE, even when $\gamma=1$.) Standard importance weighting is known to have high variance, and various techniques have been proposed to reduce variance (Precup et al., 2001; Jiang and Li, 2016; Rubinstein and Kroese, 2016; Thomas and Brunskill, 2016; Guo et al., 2017). However, these methods still exhibit exponential variance in the trajectory length (Li et al., 2015; Jiang and Li, 2016).

More related to the present paper is the recent work on off-policy RL that avoids the exponential blowup of variance. It is sufficient to adjust observed rewards according to the ratio between the target and behavior stationary distributions (Hallak and Mannor, 2017; Liu et al., 2018; Gelada and Bellemare, 2019). Unfortunately, these methods require knowledge of the behavior policy, $p(a|s)$, in addition to the transition data, which is not always available in practice. In this paper, we focus on the behavior-agnostic scenario where $p(a|s)$ is unknown. Although the recent work of Nachum et al. (2019) considers the same scenario, their approach is only applicable when the discount factor $\gamma<1$, whereas the method in this paper can handle any $\gamma\in$.

Experimental Evaluation

In this section, we demonstrate the advantages of VPM in four representative applications. Due to space limit, experiment details are provided in Appendix D.

In this subsection, we use VPM to estimate the stationary distribution of queue length. Following the standard Kendall’s notation in queueing theory (Haviv, 2009; Serfozo, 2009), we analyze the discrete-time Geo/Geo/1 queue, which is commonly used in the literature (Atencia and Moreno, 2004; Li and Tian, 2008; Wang et al., 2014). Here the customer inter-arrival time and service time are geometrically distributed with one service desk. The probe distribution $p(x)$ is a uniform distribution over the states in a predefined range $[0,B)$. The observed transition $(x,x^{\prime})$ is the length change in one time step. The queue has a closed-form stationary distribution that we can compare to (Serfozo, 2009, Sec.1.11).

Fig. 1 provides the log KL divergence between the estimated and true stationary distributions. We compare VPM to a model-based approach, which estimates the transition matrix $\widehat{{\mathcal{T}}}(x^{\prime}|x)$ from the same set of data, then simulates a long trajectory using $\widehat{{\mathcal{T}}}$. It can be seen that our method can be more effective across different sample sizes and queue configurations.

2 Solving SDEs

We next apply VPM to solve a class of SDEs known as the Ornstein-Uhlenbeck process (OUP), which finds many applications in biology (Butler and King, 2004), financial mathematics and physical sciences (Oksendal, 2013). The process is described by the equation:

where $\mu$ is the asymptotic mean, $\sigma>0$ is the deviation, $\theta>0$ determines the strength, and $W$ is the Wiener process. The OUP has a closed-form solution, which converges to the stationary distribution, a normal distribution $\mathcal{N}(\mu,\sigma^{2}/2\theta)$, as $t\to\infty$. This allows us to conveniently calculate the Maximum Mean Discrepancy (MMD) between the adjusted sample to a true sample. We compare our method with the Euler-Maruyama (EM) method (Gardiner, 2009), which is a standard simulation-based method for solving SDEs. VPM uses samples from the EM steps to train the ratio network and the learned ratio is used to compute weighted MMD.

The results are shown in Fig. 2, with different configurations of parameters $(\mu,\sigma,\theta)$. It can be seen that VPM consistently improves over the EM method in terms of the log MMD to a true sample from the normal distribution. The EM method only uses the most recent data, which can be wasteful since the past data can carry additional information about the system dynamics.

In addition, we perform experiment on real-world phylogeny studies. OUP is widely used to model the evolution of various organism traits. The results of two configurations (Beaulieu et al., 2012; Santana et al., 2012, Tab.3&1 resp.) are shown in Fig. 2(d). Notably VPM can improve over the EM method by correcting the sample with learned ratio.

3 Post-processing MCMC

In this experiment, we demonstrate how VPM can post-process MCMC to use transition data more effectively in order to learn the target distributions. We use four common potential functions as shown in the first column of Fig. 3 (Neal, 2003; Rezende and Mohamed, 2015; Li et al., 2018). A point is sampled from the uniform distribution $p(x)=\text{Unif}(x;^{2})$, then transitioned through an HMC operator (Neal et al., 2011). The transitioned pairs are used as training set $\mathcal{D}$.

We compare VPM to a model-based method that explicitly learns a transition model $\widehat{{\mathcal{T}}}(x^{\prime}|x)$, parametrized as a neural network to produce Gaussian mean (with fixed standard deviation of $0.1$). Then, we apply $\widehat{{\mathcal{T}}}$ to a hold-out set drawn from $p(x)$ sufficiently many times, and use the final instances as limiting samples (second column of Fig. 3). As for VPM, since $p$ is uniform, the estimated $\widehat{\tau}$ is proportional to the true stationary distribution. To obtain limiting samples (third column of Fig. 3), we resample from a hold-out set drawn from $p(x)$ with probability proportional to $\widehat{\tau}$.

The results are shown in Fig. 3. Note that the model-based method quickly collapses all training data into high-probability regions as stationary distributions, which is an inevitable tendency of restricted parametrized $\widehat{{\mathcal{T}}}$. Our learned ratio faithfully reconstructs the target density as shown in the right-most column of Fig. 3. The resampled data of VPM are much more accurate and diverse than that of the model-based method. These experiments show that VPM can indeed effectively use a fixed set of data to recover the stationary distribution without additional information.

To compare the results quantitatively, Fig. 4 shows the MMD of the estimated sample to a “true” sample. Since there is no easy way to sample from the potential function, the “true” sample consists of data after $2k$ HMC steps with rejection sampler. After each MCMC step, VPM takes the transition pairs as input and adjusts the sample importance according to the learned ratio. As we can see, after each MCMC step, VPM is able to post-process the data and further reduce MMD by applying the ratio. The improvement is consistent along different MCMC steps across different datasets.

4 Off-Policy Evaluation

Finally, we apply our method to behavior-agnostic off-policy evaluation outlined in Section 1, in which only the transition data and the target policy are given, while the behavior policy is unknown. Concretely, given a sample $\mathcal{D}=\left\{\left(s,a,r,s^{\prime}\right)_{i=1}^{n}\right\}$ from the behavior policy, we compose each transition in $\mathcal{D}$ with a target action $a^{\prime}\sim\pi\left(\cdot|s^{\prime}\right)$. Denoting $x=\left(s,a\right)$, the data set can be expressed as $\mathcal{D}=\left\{\left(x,x^{\prime}\right)_{i=1}^{n}\right\}$. Applying the proposed VPM with ${\mathcal{T}}(x^{\prime}|x)$, we can estimate $\frac{\mu\left(s,a\right)}{p\left(s,a\right)}$, hence the average accumulated reward can be obtained via (3). Additional derivation and discussion can be found in Appendix C.

We conduct experiments on the (discrete) Taxi environment as in Liu et al. (2018), and the challenging (continuous) environments including the Reacher, HalfCheetah and Ant.

Taxi is a gridworld environment in which the agent navigates to pick up and drop off passengers in specific locations. The target and behavior policies are set as in Liu et al. (2018). For the continuous environments, the Reacher agent tries to reach a specified location by swinging an robotic arm, while the HalfCheetah/Ant agents are complex robots that try to move forward as much as possible. The target policy is a pre-trained PPO or A2C neural network, which produces a Gaussian action distribution $\mathcal{N}(m_{t},\Sigma_{t})$. The behavior policy is the same as target policy but using a larger action variance $\Sigma_{b}=(1-\alpha)\Sigma_{t}+2\alpha\Sigma_{t},\alpha\in(0,1]$. We collect $T$ trajectories of $n$ steps each, using the behavior policy.

We compare VPM to a model-based method that estimates both the transition ${\mathcal{T}}$ and reward $R$ functions. Using behavior cloning, we also compare to the trajectory-wise and step-wise weighted importance sampling (WIST,WISS) (Precup et al., 2001), as well as Liu et al. (2018) with their public code for the Taxi environment.

The results are shown in Fig. 5. The $x$-axes are different configurations and the $y$-axes are the log Mean Square Error (MSE) to the true average target policy reward, estimated from abundant on-policy data collected from the target policy. As we can see, VPM outperforms all baselines significantly across different settings, including number of trajectories, trajectory length and behavior policies. The method by Liu et al. (2018) can suffer from not knowing the behavior policy, as seen in the Taxi environment. Weighted importance sampling methods (WIST,WISS) also require access to the behavior policy.

Conclusion

We have formally considered the problem of estimating stationary distribution of an ergodic Markov chain using a fixed set of transition data. We extended a classical power iteration approach to the batch setting, using an equivalent variational reformulation of the update rule to bypass the agnosticity of transition operator and the intractable operations in a functional space, yielding a new algorithm Variational Power Method (VPM). We characterized the convergence of VPM theoretically, and demonstrated its empirical advantages for improving existing methods on several important problems such as queueing, solving SDEs, post-processing MCMC and behavior-agnostic off-policy evaluation.

References

Appendix A Consistency of the Objectives

We just need to show $\frac{{\mathcal{T}}_{p}\tau_{t}}{p}$ is also the solution to (13). Specifically, we have

which can be attained by plugging in $\tau=\frac{{\mathcal{T}}_{p}\tau_{t}}{p}$. Finally, we conclude the proof by noticing that (13) is strictly convex so the optimal solution is unique.

Appendix B Convergence Analysis

with suitable step-sizes $\alpha_{t}\in(0,1)$, where $\epsilon\in\mathcal{L}^{2}(X)$ is a random field due to stochacity in $\widehat{{\mathcal{T}}}$. To this end, we will use the following lemma.

This can be proved by expanding both sides. Now we state our main convergence result.

Suppose $\mu_{0}\in\mathcal{L}^{2}(X)$, the step size is $\alpha_{t}=1/\sqrt{t}$, $\epsilon\in\mathcal{L}^{2}(X)$ is a random field and ${\mathcal{T}}$ has a unique stationary distribution $\mu$. After $t$ iterations, define the probability distribution over the iterations as

Then there exist some constants $C_{1},C_{2}>0$ such that

where the expectation is taken over $R$. Consequently, $\mu_{R}$ converges to $\mu$ for ergodic ${\mathcal{T}}$.

Proof Using Lemma 3 and the fact that ${\mathcal{T}}$ is non-expansive, we have

Divide both sides by $\sum_{k=1}^{t}\alpha_{k}(1-\alpha_{k})$ (taking expectation over iterations) gives

So for big enough $t$, there exists $C_{0}>0$ such that

Appendix C Application to Off-policy Stationary Ratio Estimation

We provide additional details describing how the variational power method we have developed in the main body of the paper can be applied to the behavior-agnostic off-policy estimation problem (OPE). The general framework has been introduced in Section 1 and the implementation for the undiscounted case ($\gamma=1$) is demonstrated in Section 5.4. Specifically, given a sample $\mathcal{D}=\left\{\left(s,a,r,s^{\prime}\right)_{i=1}^{n}\right\}$ from the behavior policy, we compose each transition in $\mathcal{D}$ with a target action $a^{\prime}\sim\pi\left(\cdot|s^{\prime}\right)$. Denoting $x=\left(s,a\right)$, the data set can be expressed as $\mathcal{D}=\left\{\left(x,x^{\prime}\right)_{i=1}^{n}\right\}$. Applying the proposed VPM with ${\mathcal{T}}(x^{\prime}|x)$, we can estimate $\frac{\mu\left(s,a\right)}{p\left(s,a\right)}$. Here the $\mu(s,a)=d_{\pi}(s)\pi(a|s)$ consists of the stationary state occupancy $d_{\pi}$ and the target policy $\pi$, while $p(s,a)$ is the data-collecting distribution. Then the average accumulated reward can be obtained via (3).

Here we elaborate on how the discounted case (i.e., $\gamma\in(0,1)$) can be handled by our method. We first introduce essential quantities similar to the undiscounted setting. For a trajectory generated stochastically using policy $\pi$ from an initial state $s_{0}$: $(s_{0},a_{0},r_{0},s_{1},a_{1},r_{1},\ldots)$, where $a_{t}\sim\pi(\cdot|s_{t})$, $s_{t+1}\sim P(\cdot|s_{t},a_{t})$ and $r_{t}\sim R\left(s_{t},a_{t}\right)$, the the policy value is

where $\mu_{0}$ is the initial-state distribution. Denote

Then, we can re-express the discounted accumulated reward via $\mu_{\gamma}$ and the stationary density ratio,

The proposed VPM is applicable to estimating the density ratio in this discounted case. Denoting $x=\left(s,a\right)$, $x^{\prime}=\left(s^{\prime},a^{\prime}\right)$ respectively for notational consistency, we expand $\mu_{\gamma}$ and use the definition of $d_{t}^{\pi}$:

where $\mu_{0}\pi\left(x^{\prime}\right)=\mu_{0}\left(s^{\prime}\right)\pi\left(a^{\prime}|s^{\prime}\right)$ and ${\mathcal{T}}_{p}\left(x,x^{\prime}\right)=\pi\left(a^{\prime}|s^{\prime}\right)P\left(s^{\prime}|s,a\right)p\left(s,a\right)$.

It has been shown that the RHS of (C) is contractive [Sutton and Barto, 1998, Mohri et al., 2012], therefore, the fix-point iteration,

converges to the true $\tau$ as $t\rightarrow\infty$, provided the update above is carried out exactly. Compared to (6), we can see that the RHS of (24) is now a mixture of $\mu_{0}\pi$ and ${\mathcal{T}}_{p}$, with respective coefficients $(1-\gamma)$ and $\gamma$.

Similarly, we construct the $(t+1)$-step variational update as

Compared to (11), we see that the main difference is the third term of (25) involves the initial distribution. As $\gamma\rightarrow 1$, (25) reduces to (11).

Appendix D Experiment Details

Here we provide additional details about the experiments. In all experiments, the regularization $\lambda=1$ and the optimizer is Adam with $\beta_{1}=0.5$.

For Geo/Geo/1 queue, when the arrival and finish probabilities are $q_{a},q_{f}\in(0,1)$ respectively with $q_{f}>q_{a}$, the stationary distribution is $P(X=i)=(1-\rho)\rho^{i}$ where $\rho=q_{a}(1-q_{f})/[q_{f}(1-q_{a})]$ [Serfozo, 2009, Sec.1.11]. The defaults are $(n,q_{a},q_{f})=(100,0.8,0.9)$ for the figures. $\rho$ is called traffic intensity in the queueing literature and we set $B=\lceil 40\rho\rceil$ in the experiment. The mean and standard error of the log KL divergence is computed based on 10 runs. We conduct closed-form update for $1000$ steps. As for the model-based method, we simulate the transition chain for $200$ steps to attain the estimated stationary distribution.

D.2 Solving SDEs

Using initial samples are uniformly spaced in $$, we run the Euler-Maruyama (EM) method and evaluate the MMD along the path. The$\tau$model is a neural network with 2 hidden layers of 64 units each with ReLU and Softplus for the final layer. Numbers of outer and inner steps are$T=50,M=10$. The learning rate is 0.0005. At each evaluation time step$t$, we use the most recent$1\%$of evolution data to train our model$\tau$. The plots are reporting the mean and standard deviation over 10 runs. For the phylogeny studies, the number of particles is$1k$and$dt=0.0005$for the EM simulation, while the rest settings using$dt=0.001$.

D.3 Post-processing MCMC

The potential functions are collected from several open-source projectshttps://github.com/kamenbliznashki/normalizing_flowshttps://github.com/kevin-w-li/deep-kexpfam. $50k$ examples are sampled from the uniform distribution $p(x)=\text{Unif}(x;^{2})$, then transition each $x$ through an HMC operator (one leapfrog step of size $0.5$). The $\tau$ model is a neural network with 4 hidden layers of 128 units each with ReLU activation and softplus activation for the output. The model-based $\widehat{{\mathcal{T}}}$ has a similar structure except the final layer has 2D output without activation to estimate the Gaussian mean. The mini-batch size is $B=1k$, the maximum number of power iterations $T=150$ and the number of inner optimization steps is $M=10$. The model-based ${\mathcal{T}}$ is given the same number of iterations ($MT=1500$). The learning rate is 0.001 for $\tau$ and 0.0005 for $\widehat{{\mathcal{T}}}$. To compute the model-based sample, we apply the estimated transition $100$ time steps. The MMD plot is based on a “true sample” of size $2k$ from the stationary distribution (estimated by 2k HMC transition steps). The numbers are mean and standard deviation over 10 runs. The MMD is computed by the Gaussian kernel with the median pairwise distance as kernel width.

The quality of the transition kernel and the generated data is critical. Since $x$ and $x^{\prime}$ are supposed to be related, we use an HMC kernel with one leap-frog step. The initial $x$ is effectively forgotten if using too many leap-frog steps. The main point is to show that our method can utilize the intermediate samples from the chain other than the final point. Moreover, to conform with Assumption 2, the potential functions are numerically truncated.

To verify the convergent behavior of our method, Fig. 6 shows how the ratio network improves as we train the model. It can be seen that the our method quickly concentrates its mass to the region with high potentials.

D.4 Off-policy Evaluation

Taxi is a $5\times 5$ gridworld in which the taxi agent navigates to pick up and drop off passengers in specific locations. It has a total of $2000$ states and $6$ actions. Each step incurs a $-1$ reward unless the agent picks up or drops off a passenger in the correct locations. The behavior policy is set to be the policy after $950$ Q-learning iterations and the target policy is the policy after $1000$ iterations. In the Taxi experiment, given a transition $(s,a,s^{\prime})$, instead of sampling one single action from the target policy $\pi(a^{\prime}|s^{\prime})$, we use the whole distribution $\pi(\cdot|s^{\prime})$ for estimation. We conduct closed-form update in the power method and the number of steps is $T=100$.

The results in the plots are mean and standard deviation from 10 runs.