Fast Bidirectional Probability Estimation in Markov
Models

arXiv:1507.05998v2 [cs.DS] 4 Nov 2015

Siddhartha Banerjee ∗
sbanerjee@cornell.edu

Peter Lofgren†
plofgren@cs.stanford.edu

Abstract
We develop a new bidirectional algorithm for estimating Markov chain multi-step
transition probabilities: given a Markov chain, we want to estimate the probability of hitting a given target state in ` steps after starting from a given source
distribution. Given the target state t, we use a (reverse) local power iteration to
construct an ‘expanded target distribution’, which has the same mean as the quantity we want to estimate, but a smaller variance – this can then be sampled efficiently by a Monte Carlo algorithm. Our method extends to any Markov chain on
a discrete (finite or countable) state-space, and can be extended to compute functions of multi-step transition probabilities such as PageRank, graph diffusions, hitting/return times, etc. Our main result is that in ‘sparse’ Markov Chains – wherein
the number of transitions between states is comparable to the number of states –
the running time of our algorithm for a uniform-random target node is order-wise
smaller than Monte Carlo and power iteration based algorithms;
in particular, our
√
method can estimate a probability p using only O(1/ p) running time.

1

Introduction

Markov chains are one of the workhorses of stochastic modeling, finding use across a variety of
applications – MCMC algorithms for simulation and statistical inference; to compute network centrality metrics for data mining applications; statistical physics; operations management models for
reliability, inventory and supply chains, etc. In this paper, we consider a fundamental problem associated with Markov chains, which we refer to as the multi-step transition probability estimation
(or MSTP-estimation) problem: given a Markov Chain on state space S with transition matrix P ,
an initial source distribution σ over S, a target state t ∈ S and a fixed length `, we are interested in
computing the `-step transition probability from σ to t. Formally, we want to estimate:
p`σ [t] := hσP ` , et i = σP ` eTt ,

(1)

where et is the indicator vector of state t. A natural parametrization for the complexity of MSTPestimation is in terms of the minimum transition probabilities we want to detect: given a desired
minimum detection threshold δ, we want algorithms that give estimates which guarantee small relative error for any (σ, t, `) such that p`σ [t] > δ.
Parametrizing in terms of the minimum detection threshold δ can be thought of as benchmarking
against a standard Monte Carlo algorithm, which estimates p`σ [t] by sampling independent `-step
paths starting from states sampled from σ. An alternate technique for MSTP-estimation is based on
linear algebraic iterations, in particular, the (local) power iteration. We discuss these in more detail
in Section 1.2. Crucially, however, both these techniques have a running time of Ω(1/δ) for testing
if p`σ [t] > δ (cf. Section 1.2).
∗

Siddhartha Banerjee is an assistant professor at the School of Operations Research and Information Engineering at Cornell (http://people.orie.cornell.edu/sbanerjee).
†
Peter Lofgren is a graduate student in the Computer Science Department at Stanford (http://cs.
stanford.edu/people/plofgren/).

1

1.1

Our Results

To the best of our knowledge, our work gives the first bidirectional algorithm for MSTP-estimation
which works for general discrete state-space Markov chains1 . The algorithm we develop is very
simple, both in terms of implementation and analysis. Moreover, we prove that in many settings, it
is order-wise faster than existing techniques.
Our algorithm consists of two distinct forward and reverse components, which are executed sequentially. In brief, the two components proceed as follows:
• Reverse-work: Starting from the target node t, we perform a sequence of reverse local power
iterations – in particular, we use the REVERSE-PUSH operation defined in Algorithm 1.
• Forward-work: We next sample a number of random walks of length `, starting from σ and
transitioning according to P , and return the sum of residues on the walk as an estimate of p`σ [t].
This full algorithm, which we refer to as the Bidirectional-MSTP estimator, is formalized in
Algorithm 2. It works for all countable-state Markov chains, giving the following accuracy result:
Theorem 1 (For details, refer Section 2.3). Given any Markov chain P , source distribution σ,
terminal state t, length `, threshold δ and relative error , Bidirectional-MSTP (Algorithm 2)
b `σ [t] for p`σ [t], which, with high probability, satisfies:
returns an unbiased estimate p

b `σ [t] − p`σ [t] < max p`σ [t], δ .
p
Since we dynamically adjust the number of REVERSE-PUSH operations to ensure that all residues
are small, the proof of the above theorem follows from straightforward concentration bounds.
Since Bidirectional-MSTP combines local power iteration and Monte Carlo techniques, a natural question is when the algorithm is faster than both. It is easy to to construct scenarios where the
runtime of Bidirectional-MSTP is comparable to its two constituent algorithms – for example,
if t has more than 1/δ in-neighbors. Surprisingly, however, we show that in sparse Markov chains
and for typical target states, Bidirectional-MSTP is order-wise faster:
Theorem 2 (For details, refer Section 2.3). Given any Markov chain P , source distribution σ,
length `, threshold δ and desired accuracy ; then for a uniform q
random choice of t ∈ S, the
e 3/2 d/δ), where d is the average
Bidirectional-MSTP algorithm has a running time of O(`
number of neighbors of nodes in S.
√
Thus, for typical targets, we can estimate transition probabilities of order δ in time only O(1/ δ).
Note that we do not need for every state that the number of neighboring states is small, but rather,
that they are small on average – for example, this is true in ‘power-law’ networks, where some
nodes have very high degree, but the average degree is small. The proof of this result is based on a
modification of an argument in [2] – refer Section 2.3 for details.
Estimating transition probabilities to a target state is one of the fundamental primitives in Markov
chain models – hence, we believe that our algorithm can prove useful in a variety of application
domains. In Section 3, we briefly describe how to adapt our method for some of these applications – estimating hitting/return times and stationary probabilities, extensions to non-homogenous
Markov chains (in particular, for estimating graph diffusions and heat kernels), connections to local algorithms and expansion testing. In addition, our MSTP-estimator could be useful in several
other applications – estimating ruin probabilities in reliability models, buffer overflows in queueing
systems, in statistical physics simulations, etc.
1.2

Existing Approaches for MSTP-Estimation

There are two main techniques used for MSTP-estimation. The first is a natural Monte Carlo algorithm: we estimate p`σ [t] by sampling independent `-step paths, each starting from a random
state sampled from σ. A simple concentration argument shows that for a given value of δ, we need
e
Θ(1/δ)
samples to get an accurate estimate of p`σ [t], irrespective of the choice of t, and the structure
1
Bidirectional estimators have been developed before for reversible Markov chains [1]; our method however
is not only more general, but conceptually and operationally simpler than these techniques (cf. Section 1.2).

2

of P . Note that this algorithm is agnostic of the terminal state t; it gives an accurate estimate for any
t such that p`σ [t] > δ.
On the other hand, the problem also admits a natural linear algebraic solution, using the standard
power iteration starting with σ, or the reverse power iteration starting with et (which is obtained
by re-writing Equation (1) as p`σ [t] := σ(et (P T )` )T ). When the state space is large, performing a
direct power iteration is infeasible – however, there are localized versions of the power iteration that
are still efficient. Such algorithms have been developed, among other applications, for PageRank
estimation [3, 4] and for heat kernel estimation [5]. Although slow in the worst case 2 , such local
update algorithms are often fast in practice, as unlike Monte Carlo methods they exploit the local
structure of the chain. However even in sparse Markov chains and for a large fraction of target
states, their running time can be Ω(1/δ). For example, consider a random walk on a random dregular graph and let δ = o(1/n) – then for ` ∼ logd (1/δ), verifying p`es [t] > δ is equivalent to
uncovering the entire logd (1/δ) neighborhood of s. Since a large random d-regular graph is (whp)
an expander, this neighborhood has Ω(1/δ) distinct nodes. Finally, note that as with Monte Carlo,
power iterations can be adapted to either the source or terminal state, but not both.
For reversible Markov chains, one can get a bidirectional algorithms for estimating p`es [t] based on
colliding random walks. For example, consider the problem of estimating length-2` random walk
transition probabilities in a regular undirected graph G(V, E) on n vertices [1, 6]. The main idea
is that to test if a random walk goes from s to t in 2` steps with probability ≥ δ, we can generate
two independent random walks of length `, starting from s and t respectively, and detect if they
terminate at the same intermediate node. Suppose pw , qw are the probabilities that a length-` walk
from s andPt respectively terminate at node w – then from the reversibility of the chain, we have that
p2`
σ [t] = p w∈V pw qw ; this is also the collision probability. The critical observation is that if we
generate 1/δ walks from s and t, then we get 1/δ potential collisions, which is sufficient to detect
if p2`
σ [t] > δ. This argument forms the basis of the birthday-paradox, and similar techniques used
in a variety of estimation problems (eg., see [7]). Showing concentration for this estimator is tricky
as the samples are not independent; moreover, to control the variance of the samples, the algorithms
often need to separately deal with ‘heavy’ intermediate nodes, where pw or qw are much larger than
O(1/n). Our proposed approach is much simpler both in terms of algorithm and analysis, and more
significantly, it extends beyond reversible chains to any general discrete state-space Markov chain.
The most similar approach to ours is the recent FAST-PPR algorithm of Lofgren et al. [2] for PageRank estimation; our algorithm borrows several ideas and techniques from that work. However, the
FAST-PPR algorithm relies heavily on the structure of PageRank – in particular, the fact that the
PageRank walk has Geometric(α) length (and hence can be stopped and restarted due to the memoryless property). Our work provides an elegant and powerful generalization of the FAST-PPR
algorithm, extending the approach to general Markov chains.

2
2.1

The Bidirectional MSTP-estimation Algorithm
Algorithm

As described in Section 1.1, given a target state t, our bidirectional MSTP algorithm keeps track of
a pair of vectors – the estimate vector qkt ∈ Rn and the residual vector rkt ∈ Rn – for each length
k ∈ {0, 1, 2, . . . , `}. The vectors are initially all set to 0 (i.e., the all-0 vector), except rt0 which is
initialized as et . Moreover, they are updated using a reverse push operation defined as:
Algorithm 1 REVERSE-PUSH(v, i)
Inputs: Transition matrix P , estimate vector qit , residual vectors rit , ri+1
t
1: return New estimate vectors {e
qit } and residual-vectors {reit } computed as:
eit ← qit + hrit , ev iev ;
q

e
rit ← rit − hrit , ev iev ;

e
ri+1
← ri+1
+ hrit , ev i ev P T
t
t



2
In particular, local power iterations are slow if a state has a very large out-neighborhood (for the forward
iteration) or in-neighborhood (for the reverse update).

3

The main observation behind our algorithm is that we can re-write p`σ [t] in terms of {qkt , rkt } as an
expectation over random sample-paths of the Markov chain as follows (cf. Equation (3)):
p`σ [t] = hσ, q`t i +

`
X



EVk ∼σP k r`−k
(Vk )
t

(2)

k=0

In other words, given vectors {qkt , rkt }, we can get an unbiased estimator for p`σ [t] by sampling a
length-` random trajectory {V0 , V1 , . . . , V` } of the Markov chain P starting at a random state V0
sampled from the source distribution σ, and then adding the residuals along the trajectory as in
Equation (2). We formalize this bidirectional MSTP algorithm in Algorithm 2.
Algorithm 2 Bidirectional-MSTP(P, σ, t, `max , δ)
Inputs: Transition matrix P , source distribution σ, target state t, maximum steps `max , minimum
probability threshold δ, relative error bound , failure probability pf
1: Set accuracy parameter c based on  and pf and set reverse threshold δr (cf. Theorems 1 and 2)
p
(in our experiments we use c = 7 and δr = δ/c)
2: Initialize: Estimate vectors qkt = 0 , ∀ k ∈ {0, 1, 2, . . . , `},
Residual vectors r0t = et and rkt = 0 , ∀ k ∈ {1, 2, 3, . . . , `}
3: for i ∈ {0, 1, . . . , `max } do
4:
while ∃ v ∈ S s.t. rit [v] > δr do
5:
Execute REVERSE-PUSH(v, i)
6:
end while
7: end for
8: Set number of sample paths nf = c`max δr /δ (See Theorem 1 for details)
9: for index i ∈ {1, 2, . . . , nf } do
10:
Sample starting node Vi0 ∼ σ
11:
Generate sample path Ti = {Vi0 , Vi1 , . . . , Vi`max } of length `max starting from Vi0
`
12:
For ` ∈ {1, 2, . . . , `max }: sample k ∼ U nif orm[0, `] and compute St,i
= `rt`−k [Vik ]
(We reinterpret the sum over k in Equation 2 as an expectation and sample k rather sum over
k ≤ ` for computational speed.)
13: end for
Pnf `
b `σ [t] = hσ, q`t i + (1/nf ) i=1
14: return {b
p`σ [t]}`∈[`max ] , where p
St,i
2.2

Some Intuition Behind our Approach

Before formally analyzing the performance of our MSTP-estimation algorithm, we first build some
intuition as to why it works. In particular, it is useful to interpret the estimates and residues in
probabilistic/combinatorial terms. In Figure 1, we have considered a simple Markov chain on three
states – Solid, Hollow and Checkered (henceforth (S, H, C)). On the right side, we have illustrated
an intermediate stage of reverse work using S as the target, after performing the REVERSE-PUSH
operations (S, 0), (H, 1), (C, 1) and (S, 2) in that order. Each push at level i uncovers a collection

Figure 1: Visualizing a sequence of REVERSE-PUSH operations: Given the Markov chain on the
left with S as the target, we perform REVERSE-PUSH operations (S, 0), (H, 1), (C, 1),(S, 2).
4

of length-(i + 1) paths terminating at S – for example, in the figure, we have uncovered all length
2 and 3 paths, and several length 4 paths. The crucial observation is that each uncovered path of
length i starting from a node v is accounted for in either qiv or riv . In particular, in Figure 1, all paths
starting at solid nodes are stored in the estimates of the corresponding states, while those starting at
blurred nodes are stored in the residue. Now we can use this set of pre-discovered paths to boost the
estimate returned by Monte Carlo trajectories generated starting from the source distribution. The
dotted line in the figure represents the current reverse-work frontier – it separates the fully uncovered
neighborhood of (S, 0) from the remaining states (v, i).
In a sense, what the REVERSE-PUSH operation does is construct a sequence of importancesampling weights, which can then be used for Monte Carlo. An important novelty here is that
the importance-sampling weights are: (i) adapted to the target state, and (ii) dynamically adjusted
to ensure the Monte Carlo estimates have low variance. Viewed in this light, it is easy to see how
the algorithm can be modified to applications beyond basic MSTP-estimation: for example, to nonhomogenous Markov chains, or for estimating the probability of hitting a target state t for the first
time in ` steps (cf. Section 3). Essentially, we only need an appropriate reverse-push/dynamic
programming update for the quantity of interest (with associated invariant, as in Equation (2)).
2.3

Performance Analysis

We first formalize the critical invariant introduced in Equation (2):
Lemma 1. Given a terminal state t, suppose we initialize q0t = 0, r0t = et and qkt , rkt = 0 ∀ k ≥
0. Then for any source distribution σ and length `, after any arbitrary sequence of REVERSEPUSH(v, k) operations, the vectors {qkt , rkt } satisfy the invariant:
p`σ [t]

=

hσ, q`t i

+

`
X

hσP k , rt`−k i

(3)

k=0

Proof. The proof follows the outline of a similar result in AndersenDet al. [4] for PageRank estiE
P`
mation. First, note that Equation (3) can be re-written as p`σ [t] = σ, q`t + k=0 rt`−k (P k )T .
Note that under
conditions specified
in Algorithm 2, for any σ, `, the invariant reduces
D the initial
E
P`
`−k
`
k T
`
= hσ, et (P T )` i = hσP ` , et i, which is true by definition.
to: pσ [t] = σ, qt + k=0 rt (P )
We now prove the invariant by induction. Suppose at some stage, vectors {qkt , rkt } satisfy Equation (3); to complete the proof, we need to show that for any pair (v, i), the REVERSE-PUSH(v, i)
rkt } denote the estimate and residual vectors after the
operation preserves the invariant. Let {e
qkt , e
REVERSE-PUSH(v, i) operation is applied, and define:
!
!
`
`
X
X
`−k
`−k
i
`
k T
`
k T
et +
e
∆v = q
rt (P )
− qt +
rt (P )
k=0

k=0

Now, to prove the invariant, it is sufficient to show ∆iv = 0 for any choice of (v, i). Clearly this is
true if ` < i; this leaves us with two cases:

 

• If ` = i we have: ∆iv = qet i − qit + e
rit − rit = hrkt , ev iev − hrkt , ev iev = 0
• If ` > i, we have:

 `−i−1
i+1
i
i
T
∆iv = e
ri+1
−
r
+
(e
r
−
r
)P
PT
t
t
t
t

 `−i−1
= hrkt , ev i(et P T ) − (hrkt , ev iev )P T P T
=0
This completes the proof of the lemma.
Using this result, we can now characterize the accuracy of the Bidirectional-MSTP algorithm:
Theorem 1. We are given any Markov chain P , source distribution σ, terminal state t, maximum
length `max and also parameters δ, pf and  (i.e., the desired threshold, failure probability and
5

relative error). Suppose we choose any reverse threshold δr > δ, and set the number of samplepaths nf = cδr /δ, where c = max 6e/2 , 1/ ln 2 ln (2`max /pf ). Then for any length ` ≤ `max
with probability at least 1 − pf , the estimate returned by Bidirectional-MSTP satisfies:

b `σ [t] − p`σ [t] < max p`σ [t], δ .
p
Proof. Given any Markov chain P and terminal state t, note first that for a given length ` ≤ `max ,
b `σ [t] is an unbiased estimator. Now, for any random-trajectory
Equation (2) shows that the estimate p
`
`
`
Tk , we have that the score St,k obeys: (i) E[St,k
] ≤ p`σ [t] and (ii) St,k
∈ [0, `δr ]; the first inequality again follows from Equation (2), while the second follows from the fact that we executed
REVERSE-PUSH operations until all residual values were less than δr .
P
`
Now consider the rescaled random variable Xk = St,k
/(`δr ) and X =
k∈[nf ] Xk ; then we
have that Xk ∈ [0, 1], E[X] ≤ (nf /`δr )p`σ [t] and also (X − E[X]) = (nf /`δr )(b
p`σ [t] − p`σ [t]).
Moreover, using standard Chernoff bounds (cf. Theorem 1.1 in [8]), we have that:
 2

 E[X]
P [|X − E[X]| > E[X]] < 2 exp −
and P[X > b] ≤ 2−b for any b > 2eE[X]
3
Now we consider two cases:
`
1. E[St,k
] > δ/2e (i.e., E[X] > nf δ/2e`δr = c/2e): Here, we can use the first concentration
bound to get:



 `
nf `
`
`
b σ [t] − pσ [t] ≥ pσ [t] = P |X − E[X]| ≥
P p
p [t] ≤ P [|X − E[X]| ≥ E[X]]
`δr σ

 2 
 2
 E[X]
 c
≤ 2 exp −
≤ 2 exp −
,
3
6e

where we use that nf = c`max δr /δ (cf. Algorithm 2). Moreover, by the union bound, we have:


 2 
[ 
 c
`
`
`


b σ [t] − pσ [t] ≥ pσ [t]
≤ 2`max exp −
p
P
,
32e
`≤`max


Now as long as c ≥ 6e/2 ln (2`max /pf ), we get the desired failure probability.
`
2. E[St,k
] < δ/2e (i.e., E[X] < c/2e): In this case, note first that since X > 0, we have that
`
b `σ [t] ≤ (nf /`δr )E[X] ≤ δ/2e < δ. On the other hand, we also have:
pσ [t] − p


 `

nf δ
`
b
P pσ [t] − pσ [t] ≥ δ = P X − E[X] ≥
≤ P [X ≥ c] ≤ 2−c ,
`δr
where the last inequality follows from our second concentration bound, which holds since we
have c > 2eE[X]. Now as before, we can use the union bound to show that the failure probability
is bounded by pf as long as c ≥ log2 (`max /pf ).

Combining
cases, we see that as long as c ≥ max
6e/2 , 1/ ln 2 ln (2`max /pf ), then we
hS the two
i
 `
b σ [t] − p`σ [t] ≥ max{δ, p`σ [t]} ≤ pf .
have P `≤`max p
One aspect that is not obvious from the intuition in Section 2.2 or the accuracy analysis is if using a
bidirectional method actually improves the running time of MSTP-estimation. This is addressed by
the following result, which shows that for typical targets, our algorithm achieves significant speedup:
Theorem 2. Let any Markov chain P , source
σ, maximum length `max and parameters
q distribution
2 δ
δ, pf and  be given. Suppose we set δr = `max log(`max /pf ) . Then for a uniform random choice of


q
3/2
e
t ∈ S, the Bidirectional-MSTP algorithm has a running time of O `max d/δ .

6

Proof. The runtime of Algorithm 2 consists of two parts:
Forward-work (i.e., for generating trajectories): we generate nf = c`max δr /δ sample trajectories,
each of length `max – hence the running time is O cδ`2max /δ for any Markov chain P , source
distribution σ and target
 2node t. Substituting
 for c from Theorem 1, we get that the forward-work
`max δr log(`max /pf )
running time Tf = O
. Reverse-work (i.e., for REVERSE-PUSH operations):
2 δ
Let Tr denote the reverse-work runtime for a uniform random choice of t ∈ S. Then we have:
`

E[Tr ] =

max X
1 XX
(din (v) + 1)1{REVERSE-PUSH(v,k) is executed}
|S|

t∈S k=0 v∈S

Now for a given t ∈ S and k ∈ {0, 1, . . . , `max }, note that the only states v ∈ S on which we
execute REVERSE-PUSH(v, k) are those with residual rkt (v) > δr – consequently, for these states,
we have that qkt (v) > δr , and hence, by Equation (3), we have that pkev [t] ≥ δr (by setting σ = ev ,
i.e., starting from state v). Moreover, aPREVERSE-PUSH(v, k) operation involves updating the
residuals for din (v) + 1 states. Note that t∈S pkev [t] = 1 and hence, via a straightforward counting
P
argument, we have that for any v ∈ S, t∈S 1{pke [t]≥δr } ≤ 1/δr . Thus, we have:
v

`max X
1 XX
(din (v) + 1)1{pke [t]≥δr }
[t]≥δr }
v
v
|S|
t∈S k=0 v∈S
v∈S k=0 t∈S
P




in
1 X
`max d
1
`max
v∈S d (v)
in
≤
=O
(`max + 1) · (d (v) + 1) = O
·
|S|
δr
δr
|S|
δr
v∈S
q
2 δ
to balance Tf and Tr and get the result.
Finally, we choose δr = `max log(`
max /pf )

E[Tr ] ≤

3

1
|S|

max
X `X

X

(din (v) + 1)1{pke

=

Applications of MSTP estimation

• Estimating the Stationary Distribution and Hitting Probabilities: MSTP-estimation can be
used in two ways to estimate stationary probabilities π[t]. First, if we know the mixing time τmix
of the chain P , we can directly use Algorithm 2 to approximate π[t] by setting `max = τmix and
using any source distribution σ. Theorem
q 2 then guarantees that we can estimate a stationary
3/2
probability of order δ in time O(τmix d/δ). In comparison, Monte Carlo has O(τmix /δ) runtime. We note that in practice, we usually do not know the mixing time – in such a setting, our
algorithm can be used to compute an estimate of p`σ [t] for all values of ` ≤ `max .
[t](i.e., the
An alternative is to modify Algorithm 2 to estimate the truncated hitting time pb`,hit
σ
probability of hitting t starting from σ for the first time in ` steps). ByPsetting σ = et , we get
an estimate for the expected truncated return time E[Tt 1{Tt ≤`max } ] = `≤`max `b
p`,hit
et [t]. Now,
using that fact that π[t] = 1/E[Tt ], we can get a lower bound for π[t] which converges to π[t]
as `max → ∞. We note also that the truncated hitting time has been shown to be useful in other
applications such as identifying similar documents on a document-word-author graph [9].
To estimate the truncated hitting time, we modify Algorithm 2 as follows: at each stage i ∈
eit [t] = qit [t] +
{1, 2, . . . , `max } (note: not i = 0), instead of REVERSE-PUSH(t, i), we update q
rit [t], set e
rit [t] = 0 and do not push back rit [t] to the in-neighbors of t in the (i + 1)th stage. The
remaining algorithm remains the same. It is easy to see from the discussion in Section 2.2 that the
resulting quantity pb`,hit
[t] is an unbiased estimate of P[Hitting time of t = `|X0 ∼ σ] – we omit
σ
a formal proof due to lack of space.
• Exact Stationary Probabilities in Strong Doeblin chains: A strong Doeblin chain [10] is obtained by mixing a Markov chain P and a distribution σ as follows: at each transition, the process proceeds according to P with probability α, else samples a state from σ. Doeblin chains are
widely used in ML applications – special cases include the celebrated PageRank metric [11], variants such as HITS and SALSA [12], and other algorithms for applications such as ranking [13] and
structured prediction [14]. An important property of these chains is that if we sample a starting
node V0 from σ and sample a trajectory of length Geometric(α) starting from V0 , then the terminal node is an unbiased sample from the stationary distribution [15]. There are two ways in which
our algorithm can be used for this purpose: one is to replace the REVERSE-PUSH algorithm with
7

Table 1: Datasets used in experiments
Dataset
Type
# Nodes # Edges
Pokec
directed
1.6M
30.6M
LiveJournal undirected
4.8M
69M
Orkut
undirected
3.1M
117M
Twitter-2010
directed
42M
1.5B
a corresponding local update algorithm for the strong Doeblin chain (similar to the one in Andersen et al. [4] for PageRank), and then sample random trajectories of length Geometric(α). A
more direct technique is to choose some `max >> 1/α, estimate {p`σ [t]} ∀ ` ∈ [`max ] and then
P`max `−1
α (1 − α)p`σ [t].
directly compute the stationary distribution as p[t] = `=1
• Graph Diffusions: If we assign a weight αi to random
P∞ walks of length i on a (weighted)
graph, the resulting scoring functions f (P, σ)[t] := i=0 αi σ T P i [t] are known as a graph
diffusions [16] and are used in a variety of applications. The case where αi = αi−1 (1 − α)
corresponds to PageRank. If instead the length is drawn according to a Poisson distribution
(i.e., αi = e−α αi /i!), then the resulting function is called the heat-kernel h(G, α) – this too
has several applications, including finding communities (clusters) in large networks [5].
 Note
P`max
T i
P
obeys
that for any function f as defined above, the truncated sum f `max =
α
p
i
σ
i=0
P∞
||f − f `max ||∞ ≤ `max +1 αi . Thus a guarantee on an estimate for the truncated sum directly
translates to a guarantee on the estimate for the diffusion. We can use MSTP-estimation to efficiently estimate these truncated sums. We perform numerical experiments on heat kernel estimation in the next section.
• Conductance Testing in Graphs: MSTP-estimation is an essential primitive for conductance
testing in large Markov chains [1]. In particular, in regular undirected graphs, Kale et al [6]
develop a sublinear bidirectional estimator based on counting collisions between walks in order
to identify ‘weak’ nodes – those which belong to sets with small conductance. Our algorithm can
be used to extend this process to any graph, including weighted and directed graphs.
• Local Algorithms: There is a lot of interest recently on local algorithms – those which perform
computations given only a small neighborhood of a source node [17]. In this regard, we note that
Bidirectional-MSTP gives a natural local algorithm for MSTP estimation, and thus for the
applications mentioned above – given a k-hop neighborhood around the source and target, we can
perform Bidirectional-MSTP with `max set to k. The proof of this follows from the fact
that the invariant in Equation (2) holds after any sequence of REVERSE-PUSH operations.

4

Experiments

To demonstrate the efficiency of our algorithm on large Markov chains, we use heat kernel estimation (cf. Section 3) as an example application. The heat kernel is a non-homogenous Markov chain,
defined as the probability of stopping at the target on a random walk from the source, where the walk
length is sampled from a P oisson(`) distribution. It is important in practice as it has empirically
been shown to detect communities well, in the sense that the heat kernel from s to t tends to be
larger for nodes t sharing a community with s than for other nodes t. Note that for finding complete
communities around s, previous algorithms which work from a single s to many targets are more
suitable. Our algorithm is relevant for estimating the heat kernel from s to a small number of targets.
For example, in the application of personalized search on a social network, if user s is searching for
nodes t satisfying some query, we might find the top-k results using conventional information retrieval methods for some small k, then compute the heat kernel to each of these k targets, to allow
us to place results t sharing an inferred community with s above less relevant results. A related
example is if a social network user s is viewing a list of nodes t, say the set of users attending some
event, our estimator can be used to rank the users t in that list in decreasing order of heat kernel.
This can help s quickly find the users t closest to them in the network.
We use our estimator to compute heat kernels on diverse, real-world graphs that range from millions
to billions of edges as described in Table 1. For each graph, for random (source, target) pairs, we
compute the heat kernel. We set average walk-length ` = 5 since for larger values of `, the walk
mixes into the stationary distribution, “forgetting” where it started. We set a maximum length of 10
8

Figure 2: Estimating heat kernels: Bidirectional MSTP-estimation vs. Monte Carlo, Forward Push.
For a fair comparison of running time, we choose parameters such that the mean relative error of all
three algorithms is around 10%. Notice that our algorithm is 100 times faster per (source, target)
pair than state-of-the-art algorithms.
√
standard deviations above the expectation ` + 10 ` ≈ 27; the probability of a walk being longer
than this is 10−12 , so we can ignore that event.
We compare Algorithm 2 to two state-of-the-art algorithms: The natural Monte Carlo algorithm, and
the push forward algorithm introduced in [5]. All three have parameters which allow them to trade
off speed and accuracy. For a fair comparison, we choose parameters such that the mean relative
error of all three algorithms is around 10%, and for those parameters we measure the mean running
time of all three algorithms. We implemented all three algorithms in Scala (for the push forward
algorithm, our implementation follows the code linked from [5]).
Our results are presented in Figure 4. We find that on these four graphs, our algorithm is 100x
faster (per (s, t) pair) than these state-of-the-art algorithms. For example, on the Twitter graph, our
algorithm can estimate a heat kernel score is 0.1 seconds, while the state-of-the-art algorithms both
take more than 4 minutes. We note though that the state-of-the-art algorithms were designed for
community detection, and return scores from the source to all targets, rather than just one target.
Our algorithm’s advantage applies in applications where we want the score from a source to a single
target or small set of targets.
The Pokec [18], Live Journal [19], and Orkut [19] datasets were downloaded from the Stanford
SNAP project [20]. The Twitter-2010 [21] was downloaded from the Laboratory for Web Algorithmics [22]. For reproducibility, the source code of our experiments are available on our website3 .
Acknowledgments
Research supported by the DARPA GRAPHS program via grant FA9550-12-1-0411, and by NSF
grant 1447697. Peter Lofgren was supported by an NPSC fellowship. Thanks to Ashish Goel and
other members of the Social Algorithms Lab at Stanford for many helpful discussions.

3

http://cs.stanford.edu/˜plofgren

9

References
[1] Oded Goldreich and Dana Ron. On testing expansion in bounded-degree graphs. In Studies in Complexity
and Cryptography. Miscellanea on the Interplay between Randomness and Computation. Springer, 2011.
[2] Peter Lofgren, Siddhartha Banerjee, Ashish Goel, and C Seshadhri. FAST-PPR: Scaling personalized
PageRank estimation for large graphs. In ACM SIGKDD’14, 2014.
[3] Reid Andersen, Fan Chung, and Kevin Lang. Local graph partitioning using pagerank vectors. In IEEE
FOCS’06, 2006.
[4] Reid Andersen, Christian Borgs, Jennifer Chayes, John Hopcraft, Vahab S Mirrokni, and Shang-Hua
Teng. Local computation of pagerank contributions. In Algorithms and Models for the Web-Graph.
Springer, 2007.
[5] Kyle Kloster and David F Gleich. Heat kernel based community detection. In ACM SIGKDD’14, 2014.
[6] Satyen Kale, Yuval Peres, and C Seshadhri. Noise tolerance of expanders and sublinear expander reconstruction. In IEEE FOCS’08, 2008.
[7] Rajeev Motwani, Rina Panigrahy, and Ying Xu. Estimating sum by weighted sampling. In Automata,
Languages and Programming, pages 53–64. Springer, 2007.
[8] Devdatt P Dubhashi and Alessandro Panconesi. Concentration of measure for the analysis of randomized
algorithms. Cambridge University Press, 2009.
[9] Purnamrita Sarkar, Andrew W Moore, and Amit Prakash. Fast incremental proximity search in large
graphs. In Proceedings of the 25th international conference on Machine learning, pages 896–903. ACM,
2008.
[10] Wolfgang Doeblin. Elements d’une theorie generale des chaines simples constantes de markoff. In
Annales Scientifiques de l’Ecole Normale Supérieure, volume 57, pages 61–111. Société mathématique
de France, 1940.
[11] Lawrence Page, Sergey Brin, Rajeev Motwani, and Terry Winograd. The pagerank citation ranking:
bringing order to the web. 1999.
[12] Ronny Lempel and Shlomo Moran. The stochastic approach for link-structure analysis (salsa) and the tkc
effect. Computer Networks, 33(1):387–401, 2000.
[13] Sahand Negahban, Sewoong Oh, and Devavrat Shah. Iterative ranking from pair-wise comparisons. In
Advances in Neural Information Processing Systems, pages 2474–2482, 2012.
[14] Jacob Steinhardt and Percy Liang. Learning fast-mixing models for structured prediction. In ICML’15,
2015.
[15] Krishna B Athreya and Örjan Stenflo. Perfect sampling for doeblin chains. Sankhyā: The Indian Journal
of Statistics, pages 763–777, 2003.
[16] Fan Chung. The heat kernel as the pagerank of a graph. Proceedings of the National Academy of Sciences,
104(50):19735–19740, 2007.
[17] Christina E Lee, Asuman Ozdaglar, and Devavrat Shah. Computing the stationary distribution locally. In
Advances in Neural Information Processing Systems, pages 1376–1384, 2013.
[18] Lubos Takac and Michal Zabovsky. Data analysis in public social networks. In International. Scientific
Conf. & Workshop Present Day Trends of Innovations, 2012.
[19] Alan Mislove, Massimiliano Marcon, Krishna P. Gummadi, Peter Druschel, and Bobby Bhattacharjee.
Measurement and Analysis of Online Social Networks. In Proceedings of the 5th ACM/Usenix Internet
Measurement Conference (IMC’07), San Diego, CA, October 2007.
[20] Stanford network analysis platform (snap). http://http://snap.stanford.edu/. Accessed:
2014-02-11.
[21] Paolo Boldi, Marco Rosa, Massimo Santini, and Sebastiano Vigna. Layered label propagation: A multiresolution coordinate-free ordering for compressing social networks. In ACM WWW’11, 2011.
[22] Laboratory for web algorithmics. http://law.di.unimi.it/datasets.php. Accessed: 201402-11.

10

