Non-Markovian Monte Carlo on Directed Graphs

TL;DR

This paper introduces Non-Markovian Monte Carlo (NMMC), a novel framework for sampling and inference on directed graphs, overcoming limitations of traditional MCMC methods that require undirected or reciprocal edges.

Contribution

The paper develops the first NMMC framework for directed graphs, enabling arbitrary distribution sampling and eigenvector centrality estimation through history-dependent random walks.

Findings

Successfully achieves any target distribution on directed graphs.

Estimates eigenvector centrality in a distributed manner.

Numerical experiments confirm theoretical results.

Abstract

Markov Chain Monte Carlo (MCMC) has been the de facto technique for sampling and inference of large graphs such as online social networks. At the heart of MCMC lies the ability to construct an ergodic Markov chain that attains any given stationary distribution $\boldsymbol{\pi}$, often in the form of random walks or crawling agents on the graph. Most of the works around MCMC, however, presume that the graph is undirected or has reciprocal edges, and become inapplicable when the graph is directed and non-reciprocal. Here we develop a similar framework for directed graphs, which we call Non-Markovian Monte Carlo (NMMC), by establishing a mapping to convert $\boldsymbol{\pi}$ into the quasi-stationary distribution of a carefully constructed transient Markov chain on an extended state space. As applications, we demonstrate how to achieve any given distribution $\boldsymbol{\pi}$ on a directed graph and estimate the eigenvector centrality using a set of non-Markovian, history-dependent random walks on the same graph in a distributed manner. We also provide numerical results on various real-world directed graphs to confirm our theoretical findings, and present several practical enhancements to make our NMMC method ready for practical use in most directed graphs. To the best of our knowledge, the proposed NMMC framework for directed graphs is the first of its kind, unlocking all the limitations set by the standard MCMC methods for undirected graphs.

Tables

Table 1. Table 1. Statistics of the datasets

	$\#$ nodes	$\#$ edges	$\#$ nodes (LSCC)	$\#$ edges (LSCC)
Gnutella	8,846	31,839	3,234	13,453
Slashdot	82,168	948,464	71,307	912,381
Wiki-Talk	2,394,385	5,021,410	111,881	1,477,893
Amazon	400,727	3,200,440	380,167	3,069,889

Equations

$$\sum _{i\in \mathcal{N}}x(i)A_{ij}=\lambda x(j)\text{ for all}j\in \mathcal{N}⟺\mathbf{xA}=\lambda \mathbf{x},$$

$$\nu (j)=\mathbb{P}\{X_t=j|X_0\sim \mathbf{\nu },X_t\ne 0\}$$

$$\underset{t\to \infty }{\lim }\mathbb{P}\{X_t=j|X_0\sim {\mathbf{\mu }}_0,X_t\ne 0\}={\nu }^{\prime }(j)$$

$$\nu (j)=\frac{\mathbb{P}\{X_1=j|X_0\sim \mathbf{\nu }\}}{\mathbb{P}\{X_1\ne 0|X_0\sim \mathbf{\nu }\}}=\frac{{\sum }_{i\in \mathcal{N}}\nu (i)P_{ij}}{1-{\sum }_{i\in \mathcal{N}}\nu (i)P_{i0}},$$

$$\sum _{i\in \mathcal{N}}\nu (i)P_{ij}=\lambda \nu (j),$$

$$\lambda =\mathbb{P}\{X_1\ne 0|X_0\sim \mathbf{\nu }\}=1-\sum _{i\in \mathcal{N}}\nu (i)P_{i0}.$$

$$\sum _{i\in \mathcal{N}}\pi (i){\tilde{P}}_{ij}=\sum _{i\in \mathcal{N}}\pi (i)Q_{ij}{\gamma }_{ij}=\kappa \pi (j)$$

$${\gamma }_{ij}\propto \frac{\pi (j)}{\pi (i)}\cdot \frac{{\alpha }_{ij}}{Q_{ij}},$$

$$\sum _{i\in \mathcal{N}}\pi (i)P_{ij}=\sum _{i\in S_j}\pi (i)Q_{ij}{\gamma }_{ij}=\sum _{i\in S_j}\frac{1}{c}\pi (j){\alpha }_{ij}=\frac{1}{c}\pi (j),$$

$${\gamma }_{ij}\propto \frac{\pi (j)}{\pi (i)}\cdot {\left[Q_{ij}|S_j|\right]}^{-1}=:b_{ij},$$

$${\gamma }_{ij}\propto b_{ij}=\frac{\pi (j)}{\pi (i)}\cdot \frac{d_i^{+}}{d_j^{-}},$$

$$\sum _{i\in \mathcal{N}}\nu (i)A_{ij}=(\lambda c)\cdot \nu (j)⟺\mathbf{\nu }\mathbf{A}=(\lambda c)\mathbf{\nu },$$

$$\underset{t\to \infty }{\lim }\frac{1}{t}\sum _{k=1}^{t}f(X_k)={\mathbb{E}}_{\mathbf{\pi }}\{f\}\text{mbox}a.s.,$$

$$\underset{t\to \infty }{\lim }\mathbb{E}\left\{f(X_t)|X_0\sim {\mathbf{\mu }}_0,X_t\ne 0\right\}={\mathbb{E}}_{\mathbf{\pi }}\{f\},$$

$$\nu (j)=\sum _{i\in \mathcal{N}}\nu (i)\left[P_{ij}+P_{i0}\nu (j)\right],j\in \mathcal{N}.$$

$$P_{ij}^{\mathbf{\nu }}:=P_{ij}+P_{i0}\nu (j)={\tilde{P}}_{ij}+P_{i0}\nu (j),i,j\in \mathcal{N}.$$

$${\hat{\mathbf{\mu }}}_t=\frac{1}{t+1}\sum _{k=0}^t{\delta }_{Z_k},$$

$${\hat{\mu }}_t(i)=\frac{1}{t+1}\sum _{k=0}^t{\delta }_{Z_k}(i)=\frac{1}{t+1}\sum _{k=0}^t{1}_{\{Z_k=i\}},$$

$$\mathbb{P}\{Z_{t+1}=j|{\mathcal{F}}_t\}=\sum _{i\in \mathcal{N}}{P}_{ij}^{{\hat{\mathbf{\mu }}}_t}\cdot 1_{\{Z_t=i\}}$$

$${\hat{\mathbf{\mu }}}_t=\frac{{\sum }_{k=0}^t{w}_k{\delta }_{Z_k}}{{\sum }_{k=0}^t{w}_k}=\sum _{k=0}^t{\eta }_t(k){\delta }_{Z_k},$$

$${\eta }_t(k):=\frac{w_k}{{\sum }_{k^{\prime }=0}^t{w}_{k^{\prime }}},k=0,1,2,\dots ,t.$$

$$\underset{t\to \infty }{\lim }{\hat{\mathbf{\mu }}}_t=\mathbf{\nu }\text{mbox}a.s.,\text{mbox}and\underset{t\to \infty }{\lim }\mathbb{E}\{{\hat{\mathbf{\mu }}}_t\}=\mathbf{\nu },$$

$$W_t{\hat{\mathbf{\mu }}}_t=W_{t-1}{\hat{\mathbf{\mu }}}_{t-1}+w_t{\delta }_{Z_t},$$

$${\hat{\mathbf{\mu }}}_t$$

$$=\left[1-{\eta }_t(t)\right]{\hat{\mathbf{\mu }}}_{t-1}+{\eta }_t(t){\delta }_{Z_t},$$

$${\mathbb{E}}_{{\mathbf{\pi }}^{\prime }}\left\{f\frac{\mathbf{\pi }}{{\mathbf{\pi }}^{\prime }}\right\}=\sum _{i\in \mathcal{N}}f(i)\frac{\pi (i)}{{\pi }^{\prime }(i)}{\pi }^{\prime }(i)={\mathbb{E}}_{\mathbf{\pi }}\{f\}.$$

$$\Vert {\hat{\mathbf{\mu }}}_t-\mathbf{\nu }{\Vert }_{TV}\triangleq \underset{E\subseteq \mathcal{N}}{\max }|{\hat{\mu }}_t(E)-\nu (E)|=\frac{1}{2}\sum _{i\in \mathcal{N}}|{\hat{\mu }}_t(i)-\nu (i)|.$$

$$\sqrt{t}\left(S_t/t-{\mathbb{E}}_{\mathbf{\pi }}\{f\}\right)\stackrel{d}{⟶}N(0,{\sigma }^2),\text{as}t\to \infty ,$$

$${\sigma }^2=\underset{t\to \infty }{\lim }\frac{1}{t}\mathrm{Var}\left(S_t\right)={\mathrm{Var}}_{\mathbf{\pi }}\{f\}+2\sum _{k=1}^{\infty }{\mathrm{Cov}}_{\mathbf{\pi }}\left(f(X_0),f(X_k)\right),$$

$${\gamma }_{ij}=\min \left\{1,\frac{\pi (j)Q_{ji}}{\pi (i)Q_{ij}}\right\},$$

$$\pi (i){\gamma }_{ij}{Q}_{ij}=\pi (j){\gamma }_{ji}{Q}_{ji},$$

Full text

Directed graphs, sampling, eigenvector centrality, transient Markov chains, non-reversible Markov chains, quasi-stationary distribution, history-dependent random walks

††copyright: acmcopyright††journal: POMACS††journalyear: 2019††publicationmonth: 1††ccs: Mathematics of computing Probabilistic algorithms††ccs: Mathematics of computing Stochastic processes

1. Introduction

Markov Chain Monte Carlo (MCMC) is a general technique for sampling from an arbitrary probability distribution by constructing a Markov chain that achieves the target distribution as its stationary distribution (Liu, 2004). By definition, every MCMC method requires the so-called balance condition, i.e., $\pi(j)\!=\!\sum_{i}\pi(i)P_{ij}$ for all $j$, where $P_{ij}$ is the transition probability from state $i$ to $j$ and $\bm{\pi}\!=\![\pi(1),\pi(2),\ldots,\pi(n)]$ is the stationary distribution. Most popular MCMC methods such as Metropolis-Hastings algorithm and Gibbs sampler, often referred to as Glauber dynamics, however, require a simple yet more restrictive condition, called detailed balance condition, i.e., $\pi(i)P_{ij}\!=\!\pi(j)P_{ji}$ for all $i,j$ (Hastings, 1970; Peskun, 1973; Liu, 2004). The detailed balance condition is generally preserved under a proposal and acceptance framework for state transitions. Assuming the current state $i$, a transition to state $j$ is proposed with probability $Q_{ij}$, which is the transition probability of an irreducible (proposal) Markov chain, and then accepted with acceptance probability $\gamma_{ij}$. This two-step process ensures the convergence to the desired stationary distribution $\bm{\pi}$ via achieving the detailed balance condition, i.e., $\pi(i)Q_{ij}\gamma_{ij}\!=\!\pi(j)Q_{ji}\gamma_{ji}$ for all $i,j$ 111While the proposal and acceptance framework is mostly used for the Metropolis-Hastings algorithm, (random scan) Gibbs sampler also fits into the framework where the conditional distribution of the Gibbs sampler is considered as the proposal distribution along with the acceptance probability of one. In this sense, the Gibbs sampler is a special case of the Metropolis-Hastings algorithm. Furthermore, when it comes to binary spin systems, it is essentially the Barker’s method (Barker, 1965) that was compared with the Metropolis algorithm under the proposal and acceptance framework in (Hastings, 1970; Peskun, 1973).. The detailed balance condition imposes a serious restriction, named ‘reciprocity’ condition, that $Q_{ij}\!>\!0$ if and only if $Q_{ji}\!>\!0$ (Diaconis and Saloff-Coste, 1998; Liu, 2004).

The detailed balance condition is a local condition that eases the construction of a Markov chain in MCMC methods, but is also criticized as a source of slow convergence to equilibrium and high sampling variability because every transition should be balanced out with its reverse transition. Thus, some recent advances in MCMC methods have been made to remove the requirement of the detailed balance condition in order to accelerate the convergence speed or improve sampling efficiency. A popular class of strategies is to transform a given reversible Markov chain to its related nonreversible chain by using the idea of lifting (Chen et al., 1999; Diaconis et al., 2000; Fernandes and Weigel, 2011; Turitsyn et al., 2011; Schram and Barkema, 2015; Sakai and Hukushima, 2016; Vucelja, 2016), nonbacktracking (Neal, 2004; Alon et al., 2007; Diaconis and Miclo, 2013), or vorticity (Sun et al., 2010; Chen and Hwang, 2013; Ichiki and Ohzeki, 2013; Bierkens, 2016; Rey-Bellet and Spiliopoulos, 2016; Kaiser et al., 2017). The former two ideas are based on expanding the state space to give rise to persistence of motion or directionality in state transitions. The other idea is that for the transition matrix of a given reversible chain, which is symmetric with respect to $\bm{\pi}$, a proper antisymmetric matrix with $\bm{\pi}$ (or more matrices) is added to introduce a non-trivial probability flux flowing over a cycle in one direction. They all, however, require the reciprocity condition, which is carried over from the original reversible chain before the transformation.

On the other hand, MCMC also fits into the problem of sampling from a target probability distribution over a large, unknown graph, or a huge state space whose exploration can only be done by traversing along the edges between nodes (states). The sampling problem with such limited access to the space often arises when sampling and estimating structural properties of various types of large networks, including peer-to-peer networks and online social networks, whose complete information is generally unavailable in advance. MCMC has then been a key enabler to obtain a sequence of samples by constructing an ergodic Markov chain or ‘random walk’ on a given graph, whose equilibrium distribution equals the target distribution, with only a limited and local knowledge of the graph structure (Stutzbach et al., 2009; Ribeiro and Towsley, 2010; Gjoka et al., 2011; Lee et al., 2012; Hardiman and Katzir, 2013; Avrachenkov et al., 2016; Zhou et al., 2015, 2016; Xu et al., 2014, 2017; Lee et al., 2017).

Most of the MCMC methods for constructing random walks on graphs, however, become inapplicable to directed graphs, though many real-world networks can still only be represented by directed graphs (Newman, 2003, 2010; Easley and Kleinberg, 2010; Timár et al., 2017). The famous examples are Web graphs, online social networks having follower-followee relationship, neuronal and metabolic networks, to name a few. The directed graphs naturally violate the reciprocity condition, i.e., there is at least one directed edge $(i,j)$ for which $Q_{ij}\!>\!0$ but $Q_{ji}\!=\!0$ due to the absence of its reverse directed edge $(j,i)$. The lack of reciprocity condition simply prohibits us from using the usual MCMC methods and even obtaining a closed-form expression of the stationary distribution of a given Markov chain on a directed graph. Thus, the stationary distribution has been merely estimated based on the visit counts obtained by running the Markov chain on a directed graph without knowing the explicit form of the distribution (Bar-Yossef et al., 2000; Henzinger et al., 2000; Rusmevichientong et al., 2001). This is also the very reason why popular social-network-analysis metrics defined on a directed graph such as PageRank and eigenvector centrality are generally computed and estimated by the power method, assuming that a global knowledge of the whole graph structure is available a priori (Newman, 2010).

In this paper, we present a general algorithmic framework called Non-Markovian Monte Carlo (NMMC) that constructs a set of non-Markovian random walks whose historical empirical distribution converges to any desired probability distribution $\bm{\pi}$ over a directed, non-reciprocal graph, using only local information on the graph structure. The main idea is to establish a non-trivial mapping to convert the desired distribution into the so-called quasi-stationary distribution (QSD) of a duly constructed transient Markov chain on an extended state space. We then leverage history-dependent random walks or reinforced walks to estimate the QSD “on the fly”, which does not require the entire topological information. This novel construction leads to a large degree of freedom and does not require the reciprocity condition that has been implicitly or explicitly assumed in almost all MCMC methods, and thus makes our approach the only viable solution to directed graphs. As a byproduct, our NMMC framework also allows us to estimate the eigenvector centrality of a given directed graph, one of the key metrics used in ranking or identifying influential nodes in epidemics or information spreading, by using the same set of history-dependent random-walk crawling agents.

The rest of this paper is structured as follows. In Section 2, we collect several terms and definitions on nonreversible Markov chains on directed graphs as well as on the eigenvector centrality (EVC) of a directed graph, to be used throughout the paper. In Section 3, we provide necessary background on the QSD, its key properties and characterizations, and then formally establish a mapping to convert a given distribution $\bm{\pi}$ and the EVC of a graph into the QSD of a carefully constructed transient Markov chain. In Section 4, we first explain recent developments in the literature to estimate the QSD and then build our NMMC method with theoretical guarantees to achieve a given $\bm{\pi}$ and to estimate the EVC using random crawlers on directed graphs. We also discuss several ways to improve the speed of convergence and the practicability including dynamically adjusted algorithmic parameters based on run-time estimates on the fly. Section 5 presents our extensive simulation results using a set of real-life directed graphs to numerically support our theoretical findings. In Section 6, we thoroughly compare our NMMC method to the relevant literature, ranging from the traditional MCMC methods to the PageRank algorithm. We also present how to relax (potentially impractical) assumptions with our NMMC method, which require the in-degree information and the strong connectedness of the underlying directed graph, and demonstrate the feasibility of our NMMC method as a truly working solution to most directed graphs in practice. We provide related prior studies around sampling of directed graphs in Section 7 and conclude with a brief summary capturing our novel contributions in Section 8.

2. Preliminaries

2.1. Nonreversible Markov Chains on Directed Graphs

Let $\mathcal{G}\!=\!(\mathcal{N},\mathcal{E})$ be a finite, aperiodic, strongly connected graph with a set of nodes $\mathcal{N}\!=\!\{1,2,\ldots,n\}$ and a set of directed edges $\mathcal{E}$. The graph $\mathcal{G}$ is defined by an $n\times n$ adjacency matrix $\mathbf{A}$, with elements $A_{ij}\!=\!1$ if there is an edge from node $i$ to $j$ (directed), i.e., $(i,j)\in E$, and $A_{ij}\!=\!0$ otherwise. Here $\mathbf{A}$ is asymmetric, namely, there is at least one directed edge $(i,j)$ for which $A_{ij}\!=\!1$ but $A_{ji}\!=\!0$ (non-reciprocal). Let $d^{+}_{i}\!=\!\sum_{j\in\mathcal{N}}A_{ij}$ be the out-degree of node $i$, i.e., the number of out-bound edges emanating from $i$. Similarly, let $d_{i}^{-}\!=\!\sum_{j\in\mathcal{N}}A_{ji}$ be the in-degree of node $i$ (the number of edges pointing to node $i$). We shall demonstrate how the ‘strong connectedness’ of the underlying graph $\mathcal{G}$ can be relaxed in Section 6.4.

Consider an ergodic Markov chain on $\mathcal{N}$ with its transition matrix $\mathbf{Q}\!=\![Q_{ij}]$ such that $Q_{ij}\!>\!0$ for every $(i,j)\!\in\!\mathcal{E}$, or $A_{ij}\!=\!1$, which is irreducible and aperiodic on $\mathcal{N}$. The irreducibility of $\mathbf{Q}$ is granted from the connectivity of the graph $\mathcal{G}$, or the irreducibility of the adjacency matrix $\mathbf{A}$, in that for any $i,j\!\in\!\mathcal{N}$, there exists a path following directed edges from $i$ to $j$, i.e., there exists $k\!>\!0$ with $[\mathbf{A}^{k}]_{ij}\!>\!0$. Similarly for the aperiodicity of $\mathbf{Q}$. Then, since $\mathbf{A}$ is asymmetric, the ergodic Markov chain is nonreversible without satisfying the reciprocity condition. Its unique stationary distribution $\bm{\mu}\!=\![\mu(1),\mu(2),\cdots,\mu(n)]$, which satisfies $\mu(j)\!=\!\sum_{i\in\mathcal{N}}\mu(i)Q_{ij}$ for all $j\!\in\!\mathcal{N}$, is also simply unknown for a general graph $\mathcal{G}$. For example, if we define a simple random walk on $\mathcal{G}$, it is an ergodic, nonreversible Markov chain on $\mathcal{N}$ with its transition matrix $\mathbf{Q}^{\textsf{srw}}\!=\![Q^{\textsf{srw}}_{ij}]$ given by $Q^{\textsf{srw}}_{ij}\!=\!A_{ij}/d_{i}^{+}$, where $d_{i}^{+}\!>\!0$ for all $i$. There is, however, no closed-form expression of the stationary distribution $\bm{\mu}$ for the simple random walk on a general directed graph.

2.2. Eigenvector Centrality

In addition to sampling from a target probability distribution $\bm{\pi}$ on a directed graph, another metric of interest to estimate on the directed graph is the so-called eigenvector centrality (EVC). Specifically, the EVC of node $i\!\in\!\mathcal{N}$, say $x(i)$, represents the notion of importance (centrality) of node $i$ and is proportional to the sum of the EVC values of nodes pointing to $i$, i.e.,

[TABLE]

with some positive constant $\lambda\!>\!0$, where $\mathbf{x}$ is a row vector $\mathbf{x}\!=\![x(1),x(2),\cdots,x(n)]$. Since the adjacency matrix $\mathbf{A}$ is irreducible, the Perron-Frobenius theorem (Meyer, 2000) asserts that there exists a unique solution $\mathbf{x}$ to (1), which is the left leading eigenvector of $\mathbf{A}$ with all positive entries $x(i)\!>\!0$, $i\!=\!1,2,\ldots,n$, corresponding to the largest eigenvalue $\lambda\!=\!\lambda(\mathbf{A})\!>\!0$, i.e., the spectral radius of $\mathbf{A}$. Throughout this paper, we refer to this left leading eigenvector $\mathbf{x}$ as the EVC of $\mathbf{A}$. Since $\beta\mathbf{x}$ ($\beta\neq 0$) can also be a leading eigenvector for $\lambda(\mathbf{A})$, we also assume that $\mathbf{x}$ is normalized such that $\sum_{i\in\mathcal{N}}x(i)\!=\!1$.

The EVC finds many important applications. It has been playing a crucial role in ranking nodes in social network analysis and identifying influential nodes (users) in epidemics or information spreading (Canright and Engø-Monsen, 2006; Carreras et al., 2007; Newman, 2010; Easley and Kleinberg, 2010; Lu et al., 2016), and the reciprocal of its corresponding eigenvalue $1/\lambda(\mathbf{A})$ is known to be the epidemic threshold of a general graph for SIR/SIS epidemic processes (Ganesh et al., 2005; Draief et al., 2008; Chakrabarti et al., 2008; Van Mieghem et al., 2009; Chen et al., 2016). Again, we note that there is no closed-form expression of the EVC $\mathbf{x}$. While the EVC has also been used as a main kernel for the maximal entropy random walk in the statistical physics (Burda et al., 2009; Sinatra et al., 2011), its practical construction based only on local information still remains an open issue.

Even worse, unlike the stationary distribution $\bm{\mu}$ of an ergodic Markov chain, e.g., the simple random walk, satisfying $\bm{\mu}\!=\!\bm{\mu}\mathbf{Q}$ with its (stochastic) transition matrix $\mathbf{Q}$, the EVC $\mathbf{x}$ in (1) cannot be made as a stationary distribution of an ergodic Markov chain (or a random walk) unless the adjacency matrix $\mathbf{A}$ can be made stochastic by some constant multiple, i.e., $d_{i}^{+}\!=\!d$ for all $i$. Even if this is the case, the solution $\mathbf{x}$ to (1) defies any closed-form expression since the in-degrees $d_{i}^{-}$ can be all different, reducing to the same situation with an ergodic Markov chain on a directed graph whose stationary distribution $\bm{\mu}$ is unknown. Thus, the only viable approach to compute the EVC $\mathbf{x}$ is again to rely on the power method for which the entire adjacency matrix $\mathbf{A}$ needs to be known a priori. To the best of our knowledge, it remains open or unknown as to how (or whether or not it is possible) to build crawling algorithms using only local information towards estimating the EVC $\mathbf{x}$.

3. NMMC: Mapping to QSD

The main goal of our proposed NMMC framework is to achieve any arbitrary probability distribution $\bm{\pi}$ on a general directed graph $\mathcal{G}$ and also to estimate the EVC of its adjacency matrix $\mathbf{A}$ using a set of random walks in a distributed manner without requiring the entire matrix $\mathbf{A}$ a priori. The main enabler of our NMMC framework to tackle this challenge is to convert any given ergodic, nonreversible Markov chain (e.g., simple random walk) on $\mathcal{N}$ with $\mathbf{Q}\!=\![Q_{ij}]_{i,j\in\mathcal{N}}$ and unknown $\bm{\mu}$ into a related suitable transient Markov chain on an augmented state space $\mathcal{N}^{+}\!:=\!\mathcal{N}\cup\{0\}$, where the node set $\mathcal{N}$ is a set of transient (non-absorbing) states and ‘0’ is an (artificially added) absorbing state. This mapping is made based solely on local information in a way that the resulting quasi-stationary distribution (QSD) becomes the target distribution $\bm{\pi}$. Later in Section 4 we present how to efficiently estimate the required QSD using a set of non-Markovian, history-dependent random walks, which is yet another integral part of our NMMC framework. In what follows, we first collect several definitions and basic equations we need for later use.

3.1. A Primer on QSD

Suppose that we have a discrete time transient Markov chain $\{X_{t}\}_{t\geq 0}$ on a finite state space $\mathcal{N}^{+}\!:=\!\mathcal{N}\cup\{0\}$ with its transition matrix $\mathbf{P}\!=\![P_{ij}]_{i,j\in N^{+}}$, where $\mathcal{N}\!=\!\{1,2,\cdots,n\}$ is the set of transient (non-absorbing) states and ‘0’ is the absorbing state with $P_{00}\!=\!1$, i.e., once the chain hits [math], it remains there forever. We have $P_{i0}\!>\!0$ for at least one $i\!\in\!\mathcal{N}$ to ensure that the chain is transient. Clearly, $\mathbf{P}$ is a stochastic matrix on $\mathcal{N}^{+}$. We assume that this Markov chain is irreducible and aperiodic on the set of non-absorbing states in the sense that $\mathbb{P}\{X_{t}\!=\!j~{}\mbox{for some}~{}t\!\geq\!1~{}|X_{0}\!=\!i\}\!>\!0$ and $\mbox{g.c.d.}\{t\!>\!0:P^{t}_{ii}\!>\!0\}\!=\!1$ for all $i,j\!\in\!\mathcal{N}$, where $P^{t}_{ii}\!=\!\mathbb{P}\{X_{t}\!=\!i~{}|X_{0}\!=\!i\}$ is a $t$-step transition probability.222Equivalently, this is often termed as $\mathcal{N}$ being an irreducible class (van Doorn and Pollett, 2009, 2013; Benaïm and Cloez, 2015). Then, we have $\lim_{t\to\infty}X_{t}\!=\!0$ almost surely (i.e., absorption is certain eventually). If we restrict the stochastic matrix $\mathbf{P}$ only on $\mathcal{N}$, the resulting $n\times n$ matrix $\mathbf{\tilde{P}}\!=\![\tilde{P}_{ij}]_{i,j\in\mathcal{N}}$ is naturally a sub-stochastic matrix.

Consider the dynamics of the transient Markov chain $\{X_{t}\}$ at any given time $t$ by conditioning on the event that the chain is still surviving, i.e., the chain has not hit the absorbing site ‘0’ by time $t$. These conditional probabilities $\mathbb{P}\{X_{t}\!\in\!\cdot~{}|X_{t}\!\in\!\mathcal{N}\}\!=\!\mathbb{P}\{X_{t}\!\in\!\cdot~{}|X_{t}\!\neq\!0\}$ are well-defined time-evolving full probability measures fully supported on the set of non-absorbing states $\mathcal{N}$. Now, we can consider the stationary probability measure with respect to these conditional probability measures; in other words, the invariant distribution $\bm{\nu}=[\nu(1),\nu(2),\cdots,\nu(n)]$ under the Markovian dynamics conditioned to survive. This stationary distribution $\bm{\nu}$ is called the quasi-stationary distribution (QSD) (Méléard and Villemonais, 2012; van Doorn and Pollett, 2009, 2013; Groisman and Jonckheere, 2013; Benaïm and Cloez, 2015; Grigorescu and Kang, 2013; Blanchet et al., 2016), satisfying

[TABLE]

for all non-absorbing states $j\!\in\!\mathcal{N}$ and for all time $t\!\geq\!0$, where $X_{0}\sim\bm{\nu}$ means that the initial state $X_{0}$ is drawn from the distribution $\bm{\nu}$.

Exactly in the same way as we approximate the stationary distribution of an ergodic Markov chain from its long-term limiting behavior, we can take a similar line of approach to acquire the desired QSD. Specifically, by taking the limit on the conditional distribution at time $t$ in (2) starting from an initial distribution $\bm{\mu}_{0}$, we obtain

[TABLE]

for all $j\!\in\!\mathcal{N}$. If there exists an initial distribution $\bm{\mu}_{0}$ such that $\bm{\nu}^{\prime}:=[\nu^{\prime}(1),\nu^{\prime}(2),\cdots,\nu^{\prime}(n)]$ exists, we call this limiting distribution a quasi limiting distribution (QLD). It then follows that under our finite state space setting with irreducible and aperiodic assumptions, the QSD $\bm{\nu}$ and the QLD $\bm{\nu}^{\prime}$ coincide, i.e., $\bm{\nu}\!=\!\bm{\nu}^{\prime}$, and they are uniquely determined for any initial distribution $\bm{\mu}_{0}$ (van Doorn and Pollett, 2013).

We can rewrite (2) with $t\!=\!1$ as, for any $j\!\in\!\mathcal{N}$,

[TABLE]

which then reads as

[TABLE]

where

[TABLE]

In other words, since $P_{ij}\!=\!\tilde{P}_{ij}$ for $i,j\!\in\!\mathcal{N}$, we can rewrite (4) as $\bm{\nu}\mathbf{\tilde{P}}\!=\!\lambda\bm{\nu}$, where $\mathbf{\tilde{P}}$ is the sub-stochastic matrix, which is a restriction on $\mathcal{N}$ from the original transition matrix $\mathbf{P}$ on $\mathcal{N}^{+}$. Then, by the Perron-Frobenius theorem for irreducible, non-negative matrices (not necessarily stochastic) 333The Perron-Frobenius theorem asserts that the leading (left) eigenvector of any irreducible, non-negative matrix corresponding to its largest eigenvalue is the only one with every component being positive (Meyer, 2000)., the QSD $\bm{\nu}$ must be the uniquely determined left principal eigenvector of $\mathbf{\tilde{P}}$ with $\nu(j)\!>\!0$ for all $j\!\in\!\mathcal{N}$, corresponding to the maximal eigenvalue $\lambda\!>\!0$. Note that $1-\lambda=\sum_{i\in\mathcal{N}}\nu(i)P_{i0}>0$ is the one-step absorption probability starting from the QSD $\bm{\nu}$.

3.2. Achieving $\bm{\pi}$  via QSD

We next show how to construct a transient Markov chain $\mathbf{P}\!=\![P_{ij}]$ on the augmented state space $\mathcal{N}^{+}\!=\!\mathcal{N}\!\cup\!\{0\}$ with a virtual absorbing site ‘0’ added to the original state space $\mathcal{N}$, such that the resulting QSD $\bm{\nu}$ becomes any desired target probability distribution $\bm{\pi}\!=\![\pi(1),\pi(2),\cdots,\pi(n)]$ with $\pi(i)\!>\!0$ for $i\!\in\!\mathcal{N}$, i.e., $\bm{\nu}\!=\!\bm{\pi}$. This mapping, as we shall show shortly, is always possible as long as the underlying nonreversible Markov chain $\mathbf{Q}\!=\![Q_{ij}]$ on the original state space $\mathcal{N}$ (or the directed graph $\mathcal{G}$) is ergodic. Note that letting $\bm{\mu}\!=\![\mu(1),\mu(2),\cdots,\mu(n)]$ be the stationary distribution of the ergodic chain $\mathbf{Q}$, satisfying $\bm{\mu}\!=\!\bm{\mu}\mathbf{Q}$, we do not require any knowledge on $\bm{\mu}$, i.e., $\bm{\mu}$ is unknown in our construction, as illustrated in Figure 1(a).

Adopting the proposal and acceptance framework as in usual MCMC methods, we first use the transition probability $Q_{ij}$ as a proposal distribution from the current state (node) $i$ to $j$ ($i,j\!\in\!\mathcal{N}$). This proposed move is then ‘accepted’ with probability $\gamma_{ij}$, which will be determined shortly. Unlike usual MCMC methods, however, with probability $1\!-\!\gamma_{ij}$, this move is ‘rejected’, in which case the chain gets absorbed to ‘0’ and stays there forever. We thus have $P_{ij}=Q_{ij}\gamma_{ij}$ for $i,j\!\in\!\mathcal{N}$, and $P_{i0}=\sum_{j\in\mathcal{N}}Q_{ij}(1-\gamma_{ij})$ for each $i\!\in\!\mathcal{N}$, with $P_{00}\!=\!1$. Recall that $\tilde{P}_{ij}\!=\!P_{ij}$ for $i,j\!\in\!\mathcal{N}$. See for Figure 1(b) for an illustration.

Despite the similar structure, our proposed method is fundamentally different from the MCMC methods, including the famous Metropolis-Hastings (MH) algorithm. The crucial difference is that our proposed chain $\mathbf{Q}$ is nonreversible on a directed graph without the reciprocity condition, and goes to the absorbing site ‘0’ if rejected. In contrast, in the MH algorithm the proposed chain is on a undirected graph (typically reversible), and if rejected, it simply stays at the current node and tries again in the next time step.

The next step is to choose the acceptance probability $\gamma_{ij}\in(0,1]$ such that

[TABLE]

is satisfied for all $j\!\in\!\mathcal{N}$, i.e., $\bm{\pi}\mathbf{\tilde{P}}\!=\!\kappa\bm{\pi}$ for some $\kappa\!>\!0$. Then, by the Perron-Frobenius theorem, the distribution $\bm{\pi}$ is the unique left leading eigenvector of $\mathbf{\tilde{P}}$ and thus coincides with the unique QSD $\bm{\nu}$ of the transient chain $\mathbf{P}$, corresponding to the largest eigenvalue $\lambda\!=\!\kappa$. The key aspects here are (i) how to choose $\gamma_{ij}$ using only locally available information and (ii) how to efficiently estimate or approximate the resulting QSD, or to obtain samples from the QSD in a distributed manner. In what follows, we present the specific form of $\gamma_{ij}$ that maps any desired $\bm{\pi}$ on a given directed graph to the QSD $\bm{\nu}$ of the suitably constructed transient chain $\mathbf{P}$. The case of (ii) using crawling walks on a directed graph will be discussed later in Section 4.

For a given (proposed) ergodic, nonreversible chain $\mathbf{Q}$, we define $S_{j}\!:=\!\{k\!\in\!\mathcal{N}~{}|~{}Q_{kj}\!>\!0\}$ for each $j\!\in\!\mathcal{N}$. This is the set of states (nodes) from which the proposed chain $\mathbf{Q}$ can enter into state $j$ in one step with positive probability. For each $S_{j}$, consider a probability distribution $\{\alpha_{kj}\}_{k\in S_{j}}$ such that $\alpha_{kj}\!>\!0$ for all $k\!\in\!S_{j}$ and $\sum_{k\in S_{j}}\alpha_{kj}\!=\!1$. Note that the set $S_{j}$ is the support of the probability distribution $\{\alpha_{\cdot j}\}$. We then have the following.

Theorem 3.1.

For any choice of $\{\alpha_{kj}\}_{k\in S_{j}}$ on each $S_{j}$, $j\!\in\!\mathcal{N}$, if we set

[TABLE]

such that $\gamma_{ij}\leq 1$ for any pair $(i,j)$ with $Q_{ij}>0$, then the target distribution $\bm{\pi}$ becomes the QSD $\bm{\nu}$ of the transient chain $\mathbf{P}$, where $P_{ij}\!=\!Q_{ij}\gamma_{ij}$.

Proof.

Let $\bm{\pi}$ be the given target distribution on $\mathcal{N}$. Fix $j\!\in\!\mathcal{N}$. Choose any arbitrary probability distribution $\{\alpha_{kj}\}_{k\in S_{j}}$ on the support $S_{j}$. By definition of $S_{j}$, we observe that for any $i\not\in S_{j}$, i.e., $Q_{ij}=0$, the choice of $\gamma_{ij}$ is irrelevant since the proposed move from $i$ to $j$ will never take place. We can rewrite $\gamma_{ij}$ in (7) as $\gamma_{ij}=\frac{1}{c}\cdot\frac{\pi(j)\alpha_{ij}}{\pi(i)Q_{ij}}$ with some suitable scaling constant $c>1$ such that $\gamma_{ij}\leq 1$. Now, observe that

[TABLE]

where the last equality is from $\sum_{k\in S_{j}}\alpha_{kj}\!=\!1$. Since this holds for all $j\!\in\!\mathcal{N}$ and $\tilde{P}_{ij}\!=\!P_{ij}$ for $i,j\!\in\!\mathcal{N}$, we have $\bm{\pi}\mathbf{\tilde{P}}\!=\!(1/c)\bm{\pi}$ with $1/c\!>\!0$. In view of the Perron-Frobenius theorem, $\bm{\pi}$ is the only positive eigenvector of $\mathbf{\tilde{P}}$ and is the leading left eigenvector, corresponding to the maximal eigenvalue $1/c$.

On the other hand, as explained in Section 3.1, from the Perron-Frobenius theorem we have the existence of a unique QSD $\bm{\nu}$ of the transient chain $\mathbf{P}$ that verifies $\bm{\nu}\mathbf{\tilde{P}}\!=\!\lambda\bm{\nu}$, with $\lambda$ in (5). Therefore, the unique QSD $\bm{\nu}$ must coincide with our desired probability distribution $\bm{\pi}$, with $\lambda\!=\!1/c\in(0,1)$. ∎

As a special case, we can set $\alpha_{kj}\!=\!1/|S_{j}|\!=\![\sum_{k\in\mathcal{N}}\mathbf{1}_{\{Q_{kj}>0\}}]^{-1}$ for all $k\!\in\!S_{j}$, and then have the following expression for $\gamma_{ij}$.

Corollary 3.2.

Consider

[TABLE]

such that $\gamma_{ij}\leq 1$ for each $(i,j)$ with $Q_{ij}\!>\!0$. Then, the QSD $\bm{\nu}$ of $\mathbf{P}$ becomes identical to the target distribution $\bm{\pi}$.

While the probability distribution $\{\alpha_{\cdot j}\}$ on the support $S_{j}$ can be made arbitrary as shown in Theorem 3.1, we consider $\gamma_{ij}$ in (9) throughout the rest of the paper for its simplicity. Note that the acceptance probability $\gamma_{ij}$ in (9) requires $\pi(j)/\pi(i)$, which is essentially the desired target distribution up to multiplicative constants, the proposal probability $Q_{ij}$ (clearly known), and $|S_{j}|$ (the number of possible transitions into node $j$). These are all local information around $i,j$ when the chain attempts to make a move from $i$ to $j$.

As a special case, we consider $\mathbf{Q}\!=\!\mathbf{Q}^{\textsf{srw}}$, the transition matrix of the simple random walk on a given directed graph $\mathcal{G}$. In this case, $Q_{ij}\!=\!A_{ij}/d_{i}^{+}$ and $|S_{j}|\!=\!d_{j}^{-}$, the in-degree of node $j$. This is the only required local information, not exact knowledge (identities) of the nodes whose outgoing edges point to each node. Such in-degree can be obtained without much difficulty in online social networks having follower-followee relationship such as Google++, Twitter, and Instagram, since the in-degree is usually available as part of a user profile, which is simply the number of followers (Vesdapunt and Garcia-Molina, 2016). It is worth noting that the social networks are often much more stringent on retrieving the IDs of inlinks and outlinks of each user than looking up user profiles (Vesdapunt and Garcia-Molina, 2016). We shall even demonstrate in Section 6.3 that an online estimation of the in-degree of each node would suffice for our proposed NMMC method, even if such in-degree information may not be available.

Thanks to Theorem 3.1 and Corollary 3.2, without knowing the stationary distribution $\bm{\mu}$ of the simple random walk, which is unknown in general, we can achieve any desired distribution $\bm{\pi}$ by estimating the QSD of the transient chain $\mathbf{P}\!=\![P_{ij}]$ on $\mathcal{N}^{+}$ with $P_{ij}\!=\!\gamma_{ij}Q_{ij}\!=\!\gamma_{ij}A_{ij}/d_{i}^{+}$, where

[TABLE]

with some scaling constant $c$ such that $b_{ij}/c\leq 1$ for all $(i,j)\!\in\!\mathcal{E}$. Then, if we want to achieve the uniform distribution $\bm{\pi}$ on $\mathcal{G}$, i.e., $\pi(i)\!=\!1/n$ for all $i\!\in\!\mathcal{N}$, then $c\!\geq\!\max_{(i,j)\in\mathcal{E}}d_{i}^{+}/d_{j}^{-}$. Since $d_{j}^{-}\!\geq\!1$ for any $j\in\mathcal{N}$ (since the graph is irreducible), one easy choice would be $c\!=\!\max_{i\in\mathcal{N}}d_{i}^{+}$, i.e., the maximum out-degree of the graph. When such a global constant is known in advance, our Theorem 3.1, or Corollary 3.2, applies as is. In case that such constant is not available, later in Section 4.2, we will demonstrate that a real-time (dynamic) estimation of the global constant can also be safely used.

3.3. Mapping EVC to QSD

In addition to achieving any desired distribution $\bm{\pi}$ as a QSD, we demonstrate that we can go one step further to estimate the EVC $\mathbf{x}$ of the adjacency matrix $\mathbf{A}$ of a given directed graph $\mathcal{G}$ in a similar fashion. As explained in Section 2.2, by definition, the EVC $\mathbf{x}$ is the (normalized) left leading eigenvector of the adjacency matrix $\mathbf{A}$ with the corresponding maximal eigenvalue $\lambda(A)$. Then, by simply adjusting $\gamma_{ij}$ using local information, we can equate the resulting QSD with the EVC $\mathbf{x}$ of $\mathbf{A}$, as shown next.

Corollary 3.3.

If we set $\gamma_{ij}\!=\!A_{ij}/(cQ_{ij})$ for any $(i,j)$ pair with $Q_{ij}\!>\!0$, where $c$ is some suitable scaling constant such that $\gamma_{ij}\!\leq\!1$, then the QSD $\bm{\nu}$ of the transient $\mathbf{P}$ with elements $P_{ij}\!=\!Q_{ij}\gamma_{ij}\!=\!A_{ij}/c$ for $i,j\!\in\!\mathcal{N}$, becomes the EVC $\mathbf{x}$ of the adjacency matrix $\mathbf{A}$.

Proof.

Using $\gamma_{ij}\!=\!A_{ij}/(cQ_{ij})$, and from (4) with $P_{ij}$ replaced by $\tilde{P}_{ij}$, we observe that

[TABLE]

with $\lambda$ given in (5), i.e., the resulting QSD $\bm{\nu}$ is the left eigenvector of the adjacency matrix $\mathbf{A}$ with $\nu(i)\!>\!0$, corresponding to the positive eigenvalue $\lambda c$. The assertion then follows from the Perron-Frobenius theorem. ∎

Corollary 3.3 says that the well known EVC $\mathbf{x}$ of the adjacency matrix $A$ of a given directed graph $\mathcal{G}$ can now be mapped to the unique QSD $\bm{\nu}$ of a transient Markov chain $\mathbf{P}$ driven by a general nonreversible chain $\mathbf{Q}$ using only local information for $\gamma_{ij}$, and that the corresponding largest eigenvalue $\lambda(\mathbf{A})\!=\!\lambda c$ becomes the inverse of the epidemic threshold of the given graph (See Section 2.2.) For example, when the simple random walk is used, i.e., $Q_{ij}\!=\!Q^{\textsf{srw}}_{ij}\!=\!A_{ij}/d_{i}^{+}$, we see that $\gamma_{ij}\!=\!d_{i}^{+}/c$, where $c\!=\!\max_{i\in\mathcal{N}}d_{i}^{+}$ can be set as the maximum out-degree of the graph.

So far, we have shown how to map any desired distribution $\pi$ and the EVC $\mathbf{x}$ on a directed graph onto the QSD of some transient Markov chain $\mathbf{P}$, which is obtained from any ergodic, nonreversible chain $\mathbf{Q}$, with the simple random walk as an example. Now, the question is how to efficiently simulate or estimate the given QSD, or to obtain samples from the given QSD, via the same underlying nonreversible chain $\mathbf{Q}$ used for graph exploration. In the next section, we present several methods to simulate or estimate the QSD and complete our NMMC method using a set of non-Markovian random walks (or crawling agents) to accomplish the goal.

4. NMMC: QSD Estimation

4.1. Non-Markovian Random Walks to Approximate QSD

We are often interested in calculating the expectation of a function $f:\mathcal{N}\!\to\!\mathbb{R}$ defined on the state space $\mathcal{N}$ with respect to a desired distribution $\bm{\pi}$, i.e., $\mathbb{E}_{\bm{\pi}}\{f\}\!=\!\sum_{i\in\mathcal{N}}f(i)\pi(i)$, by constructing an ergodic Markov chain $\{X_{t}\}$ that has the desired $\bm{\pi}$ as its equilibrium distribution. This is one of the main applications of the MCMC methods. By the ergodic theorem, we have

[TABLE]

for any initial distribution for $X_{0}$. Thus, one can directly utilize this defining relationship by simulating the ergodic chain $\{X_{t}\}$ for a long time and taking the time average in order to estimate or approximate $\mathbb{E}_{\bm{\pi}}\{f\}$.

In case of the QSD $\bm{\nu}$ of a transient Markov chain $\{X_{t}\}$ on $\mathcal{N}^{+}\!=\!\mathcal{N}\cup\{0\}$, we would also consider the long-time averages to estimate $\mathbb{E}_{\bm{\pi}}\{f\}$ driven by the desired QSD $\bm{\nu}\!=\!\bm{\pi}$ (or the QSD $\bm{\pi}$ itself) in view of

[TABLE]

which is from (3) and the identity between QSD and QLD in the limit for any initial distribution $\bm{\mu}_{0}$. Unfortunately, however, this does not help much. In order to rely on the regularity inherent in the limiting relationship in (13), we would have to generate many trajectories for $\{X_{t}\}$ for a long time, only to find out that most of them get absorbed or trapped to the absorbing site ‘0’ since $\mathbb{P}\{X_{t}\!\neq\!0\}$ goes to zero exponentially fast, and we are left with very few trajectories that survive for a long time and possibly contribute to useful samples for the estimation. In short, this brute-force approach simply would not work for numerical estimation of the QSD in general.

In the mathematics literature, one practical approach in approximating the QSD is based on a historical urn-based process or a reinforced random walk on the state space. Although such a process loses the Markovian property due to the dependency on its own history or trajectory for the next move, using the past trajectory of the process is amenable to practical implementation. First, observe that from (4) and (5), the QSD $\bm{\nu}$ of the transient chain $\mathbf{P}$ satisfies

[TABLE]

This can be interpreted as an ergodic Markov chain on the state space $\mathcal{N}$ with its transition matrix $\mathbf{P}^{\bm{\nu}}\!=\![P^{\bm{\nu}}_{ij}]$ given by

[TABLE]

That is to say, assuming that $\bm{\nu}$ is known, the process evolves according to $P_{ij}\!=\!\tilde{P}_{ij}$ until it gets absorbed, upon which it instantaneously ‘redistributes’ to a new node in $\mathcal{N}$ according to $\bm{\nu}$. Let $\Phi(\bm{\nu})$ be the invariant (stationary) distribution of this ergodic Markov chain with its transition matrix $\mathbf{P}^{\bm{\nu}}\!=\![P^{\bm{\nu}}_{ij}]$, indexed with $\bm{\nu}$. Then, we can see that the QSD $\bm{\nu}$ in (14) becomes precisely the fixed point of the mapping from $\bm{\nu}$ to $\Phi(\bm{\nu})$. In other words, $\bm{\nu}$ is the QSD for $\mathbf{P}$ if and only if $\Phi(\bm{\nu})\!=\!\bm{\nu}$. The construction of such an ergodic chain is, however, practically infeasible in estimating the target QSD $\bm{\nu}\!=\!\bm{\pi}$, since lack of the entire topological information prohibits direct access into the state space $\mathcal{N}$ according to $\bm{\pi}$, or $\bm{\pi}$ is simply unknown, e.g., for EVC of $\mathbf{A}$.

The above observation nonetheless motivated one of the earliest efforts in approximating the QSD of a given transient chain by Aldous et al. (Aldous et al., 1988) to replace $\bm{\nu}$ (unknown) in (14) with the (known) empirical distribution of the history of the original process.444In the literature, there are some study on Fleming-Viot (FV) processes and the evolution of the associated spatial empirical measures converging to the QSD (Grigorescu and Kang, 2004; Asselah et al., 2011; Grigorescu and Kang, 2012; Méléard and Villemonais, 2012; Grigorescu and Kang, 2013; Groisman and Jonckheere, 2013). It is essentially to replace the unknown $\bm{\nu}$ in (14) by the ‘spatial’ occupancy measure of a large number of particles undergoing identical dynamics with mutual interaction upon particle absorption. To be more precise, in our graph setting, it can be described as follows. There are $K$ particles or agents, each of which moves according to $\mathbf{P}$ of the transient Markov chain on $\mathcal{N}^{+}$. When any of the particles gets absorbed by ‘0’, it is instantaneously relocated to one of the remaining particles, chosen uniformly at random, on $\mathcal{N}$. Then, the empirical measure of the positions of the particles converges to the QSD as both the number of particles $K$ and time $t$ go to infinity. The number of particles required under the FV method, however, should be much larger than the size of the graph (i.e., $K\!\gg\!n$), making the FV method impractical to our scenario. Thus, we do not consider the FV method in this paper. To be precise, we define its historical empirical distribution $\hat{\bm{\mu}}_{t}$ by time $t$ of the underlying stochastic process $\{Z_{t}\}_{t\geq 0}$ as

[TABLE]

with $\hat{\bm{\mu}}_{t}\!=\![\hat{\mu}_{t}(1),\cdots,\hat{\mu}_{t}(n)]$, where $\hat{\mu}_{t}(i)$ is simply given by

[TABLE]

which records the relative visit frequency of $\{Z_{k}\}$ to the state $i\!\in\!\mathcal{N}$ over the first $t\!+\!1$ time steps. Note that $\hat{\mu}_{t}(i)$ is random while $\sum_{i\in\mathcal{N}}\hat{\mu}_{t}(i)\!=\!1$ for each $t$, i.e., $\hat{\bm{\mu}}_{t}$ is in fact a probability distribution (vector) on $\mathcal{N}$ for each $t$.

Consider the natural filtration $\{\mathcal{F}_{t}\}_{t\geq 0}$ associated with the historical process $\{Z_{t}\}$, i.e., $\mathcal{F}_{t}:=\sigma\{Z_{k}\,|\,k\leq t\}$. We then have

[TABLE]

for any given $j\!\in\!\mathcal{N}$, where $P_{ij}^{\hat{\bm{\mu}}_{t}}$ is from (15) with $\bm{\nu}$ replaced by $\hat{\bm{\mu}}_{t}$. That is to say, given the full history $\mathcal{F}_{t}$ of the process $\{Z_{t}\}$ up to time $t$, the process will choose its next state $Z_{t+1}$ by the original transition matrix $\mathbf{P}$ if not getting absorbed. When it hits the absorbing state ‘0’, it will instantaneously jump to another state in $\mathcal{N}$ with probability determined by the historical empirical distribution $\hat{\bm{\mu}}_{t}\!=\!\frac{1}{t+1}\sum_{k=0}^{t}\delta_{Z_{k}}$, observed by the current time instant $t$ and thus readily available to us.

This method has then been recently generalized in (Benaïm and Cloez, 2015; Benaïm et al., 2018) by utilizing the stochastic approximation theory and its induced differential equations developed for reinforced random walks (Benaïm, 1997; Pemantle, 2007), allowing different forms of the historical empirical distribution $\hat{\bm{\mu}}_{t}$ that substitutes $\bm{\nu}$ in (15). Specifically, the distribution $\hat{\bm{\mu}}_{t}$ is generalized as

[TABLE]

with

[TABLE]

Here, $\hat{\mu}_{t}(i)$ indicates the visit frequency to node $i$ by $\{Z_{k}\}$ up to time $t$, each of which is weighted by $\eta_{t}(k)$. This makes the redistribution dynamics of the process $\{Z_{t}\}$ different, while the other remains the same as above. When the process hits the absorbing state at time $t\!+\!1$, it relocates directly to one of the previously-visited nodes according to $\hat{\bm{\mu}}_{t}$ in (17). This can also be interpreted as choosing $Z_{t+1}\!:=\!Z_{k}$ with probability $\eta_{t}(k)$ in (18), $k\!=\!0,1,\ldots,t$. Setting $w_{k}\!=\!1$, or $\eta_{t}(k)\!=\!1/(t\!+\!1)$, for all $k$ in (17) reduces to (16).

The resulting process $\{Z_{t}\}$ is clearly non-Markovian and history-dependent. It is shown in (Benaïm and Cloez, 2015; Benaïm et al., 2018) that for a large class of non-decreasing weight functions $w_{k}$, or sequences $\{\eta_{t}(t)\}_{t\geq 0}$,

[TABLE]

for any arbitrary initial state $Z_{0}$ with $\hat{\bm{\mu}}_{0}\!=\!\delta_{Z_{0}}$. Note that the allowed weight function $w_{k}$ includes, for example, $w_{k}\!=\!1$ (constant), $w_{k}\!=\!k^{a}$ (polynomial), and $2^{k^{b}}$ (sub-exponential), where $a\!>\!1$ and $b\!\in\!(0,1)$. These weight functions make the sequence $\{\eta_{t}(t)\}_{t\geq 0}$ non-increasing in $t$ and satisfy the required conditions in (Benaïm and Cloez, 2015; Benaïm et al., 2018). We adopt this history-based process with reinforcement to estimate the desired QSD $\bm{\nu}\!=\!\bm{\pi}$ as an integral part of our NMMC method.

The basic operation of our NMMC method to achieve any desired distribution $\bm{\pi}$ on a directed graph $\mathcal{G}$ or to estimate the EVC $\mathbf{x}$ of the adjacency matrix $\mathbf{A}$ can be summarized as follows. Under our NMMC method, a set of non-Markovian random-walk agents are launched to crawl the graph, and each of them independently moves over $\mathcal{G}$ by choosing the next node $Z_{t+1}$ at the current node $Z_{t}\!=\!i$ according to the transition probability $P_{ij}\!=\!Q_{ij}\gamma_{ij}$, before it gets absorbed and relocated. In other words, a move from node $i$ to $j$ is proposed with probability $Q_{ij}$ following a directed edge $(i,j)$ on $\mathcal{G}$. This move is then accepted with probability $\gamma_{ij}$ as given in Theorem 3.1 and Corollary 3.2 for achieving $\bm{\pi}$ (or Corollary 3.3 for estimating the EVC $\mathbf{x}$). Otherwise, it is rejected in which case the agent is instantaneously relocated (or redistributed) to one of the previously-visited nodes according to the historical empirical distribution $\hat{\bm{\mu}}_{t}$ in (17). Figure 1 illustrates the overall operation of our NMMC method.

We note that the history-dependent/reinforced random walks are quite relevant and amenable to practical implementation on most popular online social networks, since the crawler (e.g., a simple random walk) continues to collect the ID (e.g., URL) of each node visited as it explores the graph and can always relocate/redistribute to one of the previously-visited nodes, while keeping track of which nodes have been visited up until now. The historical empirical distribution governs not only the dynamics of each crawling agent, but it also enables us to estimate the target QSD $\bm{\nu}\!=\!\bm{\pi}$. Therefore, our NMMC method provides the same benefit and fruitful denouement that we cherished from the standard MCMC methods that are rooted in the ergodic theorem in (12).

4.2. Performance Enhancements

So far, we have demonstrated how any ergodic (inherently nonreversible) Markov chain $\mathbf{Q}$ with unknown $\bm{\mu}$ can be used to achieve a given desired distribution $\bm{\pi}$ on a directed graph $\mathcal{G}$ and to estimate the EVC $\mathbf{x}$ of the adjacency matrix $\mathbf{A}$. It is essentially achieved by combining our mapping $\bm{\pi}$ or $\mathbf{x}$ to the QSD $\bm{\nu}$ of a related transient Markov chain and the recent development in the mathematics literature as to estimating the QSD via a set of non-Markov chain, history-dependent random walks.

As shall be seen later in Section 5, the direct application of the non-Markovian random walks with historical empirical distribution $\hat{\bm{\mu}}_{t}$ may not be satisfactory due to a possible imbalance between their diffusion and redistribution, both of which are the necessary integral components. In this section, we present several ways to make the crawling walks achieve a right balance between their diffusion and redistribution, which translates into speeding up the convergence of $\hat{\bm{\mu}}_{t}$ to the desired QSD $\bm{\nu}\!=\!\bm{\pi}$ (or $\bm{\nu}\!=\!\mathbf{x}$).

Initial empirical measure $\hat{\bm{\mu}}_{0}$. We next revisit the historical empirical distribution $\hat{\bm{\mu}}_{t}$ in (17). Letting $W_{t}\!:=\!\sum^{t}_{k=0}w_{k}$ for a given positive weight sequence, $\{w_{k}\}_{k\geq 0}$, we observe that (17) can be written as, for $t\geq 1$,

[TABLE]

which then leads to

[TABLE]

from the identities $W_{t}\!=\!W_{t-1}\!+\!w_{t}$ and $\eta_{t}(t)\!=\!w_{t}/W_{t}$. This implies that the historical empirical distribution $\hat{\bm{\mu}}_{t}$ can be computed recursively based on (21). Although the convergence of $\hat{\bm{\mu}}_{t}$ to the QSD $\bm{\nu}$ was shown for any arbitrary initial state $Z_{0}$ with $\hat{\bm{\mu}}_{0}\!=\!\delta_{Z_{0}}$ in (Benaïm and Cloez, 2015; Benaïm et al., 2018), its proof based on the stochastic approximation theory and the induced differential equations in fact does not require the initial distribution $\hat{\bm{\mu}}_{0}$ in (21) to be in the form of $\delta_{Z_{0}}$, where $\hat{\bm{\mu}}_{0}$ has a value of 1 at the element $Z_{0}$ and zero at all other elements, but can be made arbitrary as long as $\hat{\bm{\mu}}_{0}$ is a probability distribution (vector) on $\mathcal{N}$.

This freedom on the choice of $\hat{\bm{\mu}}_{0}$ gives us a great degree of flexibility. Clearly, the choice of $\hat{\bm{\mu}}_{0}$ affects the whole evolution of the process $\{Z_{t}\}$ for all time $t$ via redistribution mechanism, and thus can be made to our advantage. For example, a judicious choice of $\hat{\bm{\mu}}_{0}$ can translate into a situation where a crawling agent $\{Z_{t}\}$ can start with an ‘already-established’ history to some degree or more generally, with any arbitrary (normalized) visit counts assigned to a subset of nodes on the graph to begin with. Thus, our proposed NMMC method can be generalized with arbitrary $\hat{\bm{\mu}}_{0}$ and its overall procedure is now summarized in Algorithm 1. Here the proposed chain $\mathbf{Q}$, the initial empirical measure $\hat{\bm{\mu}}_{0}$, and the (unnormalized) weight function $w_{k}$ should be specified as inputs into Algorithm 1. Note that the initial position of a crawling agent $Z_{0}$ is randomly chosen according to $\hat{\bm{\mu}}_{0}$.

Dynamic and adaptive choice of parameters. Our proposed NMMC method is mathematically guaranteed but its speed of convergence depends on the form of the acceptance probability $\gamma_{ij}$, along with $Q_{ij}$ and the target QSD, which interweaves diffusion and redistribution – two seemingly opposite (yet integral) components of the method. As we have seen from Section 3, the specific form of $\gamma_{ij}$ is different depending on the choice of $\mathbf{Q}\!=\![Q_{ij}]$ and what to estimate via the target QSD, i.e., the target distribution $\bm{\pi}$ or the EVC $\mathbf{x}$ of the adjacency matrix $\mathbf{A}$. It is particularly desirable to maintain some reasonable values of $\gamma_{ij}$ over almost every $(i,j)\in\mathcal{E}$, or more precisely, transition pair $(i,j)$ with $Q_{ij}\!>\!0$, so as to strike the right balance between diffusion and redistribution to achieve faster convergence.

We take a closer look at the acceptance probability $\gamma_{ij}$ as a function of $Q_{ij}$ and the target QSD. Recall that from Theorem 3.1 and Corollary 3.2, it is given by $\gamma_{ij}\!=\!b_{ij}/c$ in (9) for achieving (or sampling from) a given distribution $\bm{\pi}$ on the directed graph $\mathcal{G}$, and the corresponding normalizing constant $c$ can be obtained based on $c\!=\!\max_{i,j}b_{ij}$ (the smallest possible value for $c$). Similarly, for estimating the EVC $\mathbf{x}$ of the adjacency matrix $\mathbf{A}$, it can be chosen as $c\!=\!\max_{i,j}1/Q_{ij}$ for $\gamma_{ij}\!=\!A_{ij}/(cQ_{ij})$ from Corollary 3.2. The global constant $c$ over the graph, however, may not be available in advance or could take a long time to estimate, although its real-time (dynamic) estimation can be incorporated as shall be explained shortly. Furthermore, the resulting values of $\gamma_{ij}$ may vary over $(i,j)\!\in\!\mathcal{E}$ and they can also be small for some $(i,j)$ pairs. Such small acceptance probabilities in turn make the agent get absorbed to ‘0’ and redistributed to the already visited nodes. The redistribution mechanism is of course an important and necessary component of our algorithm to shape up the historical empirical distribution $\hat{\bm{\mu}}_{t}$ toward the target QSD. But, if it is too strong (very small $\gamma_{ij}$), the agent may virtually get stuck on a small set of visited nodes before diffusing out to a potentially vast unvisited portion of the graph (state space) for exploration.

To remedy this possibly slow diffusion problem, we observe that the target QSD $\bm{\nu}\!=\!\bm{\pi}$ when achieving any desired distribution $\bm{\pi}$ is completely controllable to our benefit. To be precise, if it is of primary concern to estimate $\mathbb{E}_{\bm{\pi}}\{f\}\!=\!\sum_{i\in\mathcal{N}}f(i)\pi(i)$ for a given sampling function $f:\mathcal{N}\!\to\!\mathbb{R}$, which is often the case for ‘undirected’ graph sampling (Stutzbach et al., 2009; Ribeiro and Towsley, 2010; Gjoka et al., 2011; Lee et al., 2012; Hardiman and Katzir, 2013; Avrachenkov et al., 2016), the target distribution can be controlled together with a so-called importance sampling method in MCMC (Liu, 2004). With an alternative sampling distribution $\bm{\pi}^{\prime}$, it is to estimate

[TABLE]

Thus, we can choose the target distribution $\bm{\pi}$ along with $\mathbf{Q}$ so that it translates into more balanced $\gamma_{ij}$. As shall be demonstrated later in Section 5, with $\mathbf{Q}\!=\!\mathbf{Q}^{\textsf{srw}}$, setting the target $\bm{\pi}$ to be proportional to node in-degrees, i.e., $\pi(i)\propto d^{-}_{i},i\!\in\!\mathcal{N}$, turns out to be very effective in achieving fast convergence speed. This improvement is somewhat expected since the simple random walk would visit high in-degree nodes more frequently, which helps the historical empirical distribution $\hat{\bm{\mu}}$ get closer to this target $\bm{\pi}$ quickly.

To further alleviate the slow diffusion problem and also enable online estimation of the normalizing constant $c$, we allow $c$ to be time-varying, say $c_{t}$, so that it can be learned through graph exploration by the crawling agent. For ease of exposition, we focus on achieving the desired $\bm{\pi}$ over the graph, but everything can easily be carried over to estimating the EVC $\mathbf{x}$ of $\mathbf{A}$. We set the acceptance probability $\gamma_{ij}(t)\!:=\!b_{ij}/c_{t}$, which is updated in a way that $c_{t}\!:=\!b_{ij}$ at every time $t$ if $c_{t-1}\!<\!b_{ij}$ for a pair of the current node $i$ and its neighbor $j$ as a potential next node, and $c_{t}\!:=\!c_{t-1}$ if otherwise. In other words, the agent keeps track of the maximum $c_{t}$ over all the values of $b_{ij}$ for node pairs that have been discovered up to time $t$. Note that $c_{0}\!=\!1$ initially.

We can take one step further to foster the initial diffusion speed to a greater extent, by intentionally slowing down the growth of $c_{t}$ for higher values of $\gamma_{ij}(t)$ (being closer to one). Specifically, with probability $p\!>\!0$, we perform the aforementioned update of $c_{t}$, while with probability $1\!-\!p$, the time-varying normalizing constant $c_{t}\!=\!c_{t-1}$ is kept the same even if $c_{t-1}\!<\!b_{ij}$, in which case the proposed move to $j$ is always accepted. The entire procedure of ‘dynamic’ NMMC method is summarized in Algorithm 2, which replaces the lines 3–10 of Algorithm 1. Note that the resulting acceptance probability becomes $\gamma_{ij}(t)\!=\!\min\{1,b_{ij}/c_{t}\}$ as seen from Algorithm 2.

While the NMMC method with time-varying $c_{t}$ is seemingly heuristic, it is in fact well supported by Theorem 3.1 and our observation on the initial empirical measure $\hat{\bm{\mu}}_{0}$, i.e., the theoretical guarantee remains intact. To see this, we define a stopping time as $t^{*}\!:=\!\min\{t\!\geq 0\!~{}|~{}c_{t}\!=\!c\}$. Since the graph is finite, this $t^{*}$ is almost surely finite for $p\!>\!0$. We would run the NMMC method with $c_{t}$ until $t^{*}$, and then reset the time clock to zero at $t^{*}$ and treat the historical empirical distribution $\hat{\bm{\mu}}_{t^{*}}$ to be newly set as $\hat{\bm{\mu}}_{0}$. It then follows from our observation pertaining to $\hat{\bm{\mu}}_{0}$, along with Theorem 3.1, that the historical empirical measure $\hat{\bm{\mu}}_{t}$ converges to the given target QSD $\bm{\nu}\!=\!\bm{\pi}$. In addition to the theoretical guarantees, we next demonstrate that our NMMC method with a suitable choice of the $c_{t}$-updating probability $p$ indeed greatly improves the speed of convergence to the target QSD $\bm{\nu}\!=\!\bm{\pi}$ by effectively combating the slow-diffusion problem triggered by overly frequent redistribution to already-visited sites.

5. Simulation Results

In this section, we present simulation results to support our theoretical findings. To this end, we consider four real-world network datasets available in the SNAP repository (Leskovec and Krevl, 2014) as test-case directed graphs, whose statistics are summarized in Table 1. The numbers of nodes and edges in the largest strongly connected component (LSCC) of each graph are given in the fourth and fifth columns respectively. Note that self-loops, if present, e.g., in the original dataset of Slashdot graph, are removed for our simulations.

The primary focus here is, through numerical simulations, to confirm that the historical empirical distribution $\hat{\bm{\mu}}_{t}$ by our NMMC method in Algorithms 1 and 2 (the original one and its dynamic extension) converges to the target QSD $\bm{\nu}\!=\!\bm{\pi}$ (or $\bm{\nu}\!=\!\mathbf{x}$) and to evaluate the speed of convergence of $\hat{\bm{\mu}}_{t}$ to $\bm{\nu}$. As a performance metric, we use the total variation distance (TVD)555Note that our purpose here is completely different from the typical use of TVD for characterizing the speed of convergence of an ergodic, reversible Markov chain on a finite state-space, say $\mathbf{Q}$, to its stationary distribution, say $\bm{\mu}$. To be precise, for any initial $\bm{\mu}_{0}$, $\|\bm{\mu}_{0}\mathbf{Q}^{t}-\bm{\mu}\|_{TV}\leq\kappa\rho^{t}$ for some constant $\kappa$, where $\rho$ is the second largest eigenvalue modulus of $\mathbf{Q}$ with $\rho\!<\!1$ (Brémaud, 1999). This exponential convergence speed is for the ‘marginal’ distribution $\bm{\mu}_{0}\mathbf{Q}^{t}$ at time $t$ to $\bm{\mu}$ as $t$ increases. In contrast, we are interested in the convergence of the historical empirical distribution $\hat{\bm{\mu}}_{t}$, which encompasses the entire history of the underlying process with proper time weights. between the historical empirical distribution $\hat{\bm{\mu}}_{t}$ and the target QSD $\bm{\nu}$, which is given by

[TABLE]

The inputs of Algorithms 1 and 2 are chosen as follows. The initial position $Z_{0}$ of each agent is randomly determined, with the initial empirical measure $\hat{\bm{\mu}}_{0}\!=\!\delta_{Z_{0}}$. For the target QSD $\bm{\nu}$ to be mapped and estimated, we consider (i) uniform distribution, i.e., $\pi(j)\!=\!1/n$, (ii) in-degree distribution, i.e., $\pi(j)\!\propto\!d^{-}_{j}$, and (iii) the EVC $\mathbf{x}$ of the adjacency matrix $\mathbf{A}$. We also consider $\mathbf{Q}\!=\!\mathbf{Q}^{\textsf{srw}}$, where $Q_{ij}\!=\!A_{ij}/d^{+}_{i}$ for $i,j\!\in\!\mathcal{N}$. We use the LSCC of each graph for simulations of all these cases. In order to understand how different choices of unnormalized weight functions $w_{k}$ would affect the speed of convergence, we vary weight functions as $w_{k}\!=\!1,k^{1},k^{3},k^{5},k^{10},2^{\sqrt{k}}$, with increasing tendency to relocate to recently-visited nodes upon redistribution. For each simulation we use a set of independent crawling agents, which are independently moving over the graph. While each agent maintains its own historical empirical distribution $\hat{\bm{\mu}}_{t}$ for redistribution in case of individual absorption, all the historical empirical distributions are combined and normalized as a unified distribution for the convergence test.

We first present simulation results for Gnutella graph. We use 100 independent crawling agents for all the results. Figure 2(a) shows the TVD results when achieving a uniform distribution ${\mathbf{u}}$ under the NMMC method in Algorithm 1, where the global constant $c$ is known a priori. Recall that for achieving $\bm{\pi}\!=\!{\mathbf{u}}$, the acceptance probability $\gamma_{ij}$ in (9) becomes $\gamma_{ij}\!=\!b_{ij}/c$ with $b_{ij}\!=\!d^{+}_{i}/d^{-}_{j}$ and $c\!=\!\max_{i,j}b_{ij}$. This global constant $c$ can be as large as $c\!=\!\max_{i}d_{i}^{+}$ and this was the case in this simulation. Figure 2(a) shows that $\hat{\bm{\mu}}_{t}$ is getting closer to the target QSD $\bm{\nu}\!=\!\bm{\pi}$ for all considered weight functions $w_{k}$. In addition, Figure 2(a) exhibits that heavier emphasis on latter visited sites leads to faster convergence of $\hat{\bm{\mu}}_{t}$ to $\bm{\pi}$. We would expect that it is attributed to higher weights of being relocated toward more recently-visited nodes that gives rise to faster diffusion on $\mathcal{G}$ and quicker exploration of the graph. Nonetheless, the overall performance may be considered unsatisfactory as the resulting TVD values are still high. This means that the redistributions still take place too often, resulting in an imbalance between diffusion and redistribution.

In Figure 2(b)–(d), we present the TVD results for estimating the same uniform ${\mathbf{u}}$ under the dynamic NMMC method in Algorithm 2 while varying the value of the $c_{t}$-updating probability $p$ from 1 to 0.01. It turns out that, with the choice of $p\!=\!0.01$, the TVD is reduced by more than 0.1 at every time step when compared to that of $p\!=\!1$. To achieve the same precision in terms of the TVD, the number of time steps $t$ required for $p\!=\!1$ would be much larger than that of $p\!=\!0.01$ by at least an order of magnitude. Observing a clear decreasing trend with $p$, one may think that $p\!=\!0$ would be the best. However, when $p\!=\!0$, $c_{t}$ is never updated over $t$ as seen from our Algorithm 2, i.e., $c_{t}\!=\!c_{0}\!=\!1$ for all $t$. Still, it is possible that $\gamma_{ij}(t)\!=\!\min\{1,b_{ij}/c_{t}\}\!=\!\min\{1,b_{ij}\}$ is less than one, meaning that a proposed move is not always accepted. Nonetheless, since $c_{t}$ never approaches the true value of $c$ when $p\!=\!0$, there is no theoretical guarantee for the convergence either. Accordingly, the resulting TVD does not decrease as shown in Figure 3(a), where $w_{k}\!=\!1$ is considered. Furthermore, in view of faster diffusion with smaller $p$ leading to a better speed of convergence, we can take another extreme case by simply running the simple random walk all the time, which always diffuses without any redistribution mechanism. It turns out that the resulting TVD again does not decrease as shown in Figure 3(b). The observations so far collectively indicate an importance of having a right balance between diffusion and redistribution to shape up the historical empirical distribution $\hat{\bm{\mu}}_{t}$ toward the target QSD and to achieve faster convergence.

We further take a closer look into the cases with $w_{k}\!=\!k^{\alpha}$ and $w_{k}\!=\!2^{\sqrt{k}}$. While the case with $w_{k}\!=\!2^{\sqrt{k}}$ is the best with no dynamic update of $c$ as seen from Figure 2(a), it becomes worse than the case with $w_{k}\!=\!k^{\alpha}$ when using the dynamic NMMC method with $c_{t}$-updating probability $p$ as shown in Figure 2(b)–(c). The case with $w_{k}\!=\!2^{\sqrt{k}}$ when combined with the dynamic updates of $c_{t}$ places over-emphasis on diffusion, which causes more harm than help. This again signifies the importance of striking the right balance between diffusion and redistribution of the reinforced walk $Z_{t}$ toward faster convergence, although the presence of such a crossover would depend on the choice of the directed graph $\mathcal{G}$, target QSD $\bm{\nu}$ and the proposal chain $\mathbf{Q}$. Furthermore, we have observed that the case with $w_{k}\!=\!k^{\alpha}$ is always better with the use of $c_{t}$-updating probability $p$ than the case with $w_{k}\!=\!2^{\sqrt{k}}$, when achieving the in-degree distribution, i.e., $\pi(j)\!\propto\!d^{-}_{j}$, and the EVC $\mathbf{x}$ of the adjacency matrix $\mathbf{A}$. We omit the results for brevity and exclude the case with $w_{k}\!=\!2^{\sqrt{k}}$ for the subsequent numerical evaluations. Note that the acceptance probability $\gamma_{ij}$ in (9) becomes $\gamma_{ij}\!=\!d^{+}_{i}/(cd^{-}_{i})$ with $c\!=\!\max_{i}d^{+}_{i}/d^{-}_{i}$ for the former, and $\gamma_{ij}\!=\!d^{+}_{i}/c$ with $c\!=\!\max_{i}d^{+}_{i}$ for the latter.

We instead report the results under the (unnormalized) weight function $w_{k}\!=\!k^{\alpha}$ with different values of $\alpha$ and $p\!=\!0.01$ for all three cases of the target QSD in Figure 2(d)–(h). For achieving $\bm{\pi}\!=\!{\mathbf{u}}$, the case with $\alpha\!=\!1$ exhibits the fastest decaying rate in its TVD curve initially, but later becomes slower than the other cases. Similar tendency with $\alpha\!=\!1$ is observed for achieving $\pi(j)\!\propto\!d^{-}_{j}$ and estimating the EVC $\mathbf{x}$, as shown in Figure 2(e)-(h). On the other hand, we see that the case with $\alpha\!=\!10$ eventually becomes the best with the fastest decaying rate, as we can see from the inset figure of Figure 2(d). We also observe that the eventual winner is the case with $\alpha\!=\!5$ for achieving $\pi(j)\!\propto\!d^{-}_{j}$, while the case with $\alpha\!=\!3$ eventually appears to be the best for estimating the EVC $\mathbf{x}$. This set of simulation results indicate that target QSD $\bm{\nu}$ clearly affects the extent of diffusion to achieve the right balance between diffusion and redistribution.

In addition, we observe that TVD values decay to zero polynomially fast, as evidenced in Figure 2(f), Figure 2(h), and the inset figure of Figure 2(d), all drawn on a log-log scale. We have also observed such a polynomial decaying behavior for the other choices of parameters. This behavior can be explained as follows.

Suppose that we have a trajectory of an ergodic Markov chain, say $\{X_{t}\}_{t\geq 0}$, with its stationary distribution $\bm{\pi}$. For a function $f:\mathcal{N}\!\to\!\mathbb{R}$, letting $S_{t}\!:=\!\sum_{k=1}^{t}f(X_{k})$, $S_{t}/t$ converges to $\mathbb{E}_{\bm{\pi}}\{f\}$ almost surely as seen from (12) as $t$ grows. In addition, the standard Central Limit Theorem (CLT) for ergodic Markov chains says that for any bounded function $f$ and any initial distribution on $X_{0}$,

[TABLE]

where $N(0,\sigma^{2})$ is a Gaussian random variable with zero mean and variance $\sigma^{2}$. The (asymptotic) variance $\sigma^{2}$ is given by (Liu, 2004; Lee et al., 2012; Avrachenkov et al., 2016)

[TABLE]

where the first term in RHS is the marginal variance of $f$ with respect to $\bm{\pi}$ and the second term in RHS is the lag-$k$ autocovariance of the stationary sequence $\{f(X_{t})\}$. From (24), we see that the fluctuations of $S_{t}/t$ around the mean value $\mathbb{E}_{\bm{\pi}}\{f\}$ are roughly within the interval of length of order $\sigma/\sqrt{t}$. Even when $X_{t}$’s are $i.i.d.$ over $t$, each having distribution $\bm{\pi}$ (i.e., perfect independent samples all the time), the asymptotic order of convergence from $S_{t}/t$ to $\mathbb{E}_{\pi}\{f\}$ still remains the same as $O(1/\sqrt{t})$, while $\sigma^{2}$ reduces to $\mathrm{Var}_{\bm{\pi}}\{f\}$.

What we have in this paper is the convergence of $\hat{\bm{\mu}}_{t}\!\in\!\Omega$ to the target $\bm{\nu}\!\in\!\Omega$ almost surely and the typical size of the fluctuations of the difference $\hat{\bm{\mu}}_{t}-\bm{\nu}$, where $\Omega\!=\!\{\bm{\nu}\in[0,1]^{n}|\sum_{i=1}^{n}\nu(i)\!=\!1\}$ is the set of all probability distributions on $\mathcal{N}$. If the underlying process on $\mathcal{N}$ were to be a Markov chain, by resorting to similar CLT, we would obtain similar $O(1/\sqrt{t})$ convergence speed where the typical size of the fluctuations is now characterized by that of a $n$-dimensional multivariate Gaussian random vector with some $n\times n$ covariance matrix. However, the underlying process $Z_{t}$ in our framework is history-dependent and thus non-Markovian, so we would not expect anything faster than $O(1/\sqrt{t})$ for its error as seen from Markovian or $i.i.d.$ cases. Nonetheless, the figures drawn on a log-log scale in our setting demonstrate that the error decays still polynomially fast $O(t^{-\beta})$. These observations are further supported by the mathematical results in (Benaïm and Cloez, 2015), which apply to the static NMMC method. Specifically, Corollary 1.3 in (Benaïm and Cloez, 2015) implies that if one considers a class of polynomial weight functions $w_{k}\!=\!k^{\alpha}$, $\alpha\!\geq\!0$, then the difference $\hat{\bm{\mu}}_{t}-\bm{\nu}$ in the $l^{1}$ norm is bounded above by $Ct^{-\theta}$ almost surely for some (random) constant $C\!>\!0$, where $\theta\!\leq\!\frac{1}{2}$.666In other words, for every possible sample path $\omega$ of the underlying stochastic process $\{Z_{t}\}$ by the crawler, there exists some random constant $C(\omega)$ such that the difference $\hat{\bm{\mu}}_{t}-\bm{\nu}$ in the $l^{1}$ norm is bounded above by $C(\omega)t^{-\theta}$ almost surely. Note also that $\theta$ depends on the extremal value of the spectral functional, which requires the knowledge of the entire spectrum of the underlying graph structure. We refer to (Benaïm and Cloez, 2015) for more details. In addition, Theorem 1.4 in (Benaïm and Cloez, 2015) shows that the CLT from $\hat{\bm{\mu}}_{t}$ to $\bm{\nu}$ holds with respect to a rescaled time $\eta_{t}(t)^{-1}\!\sim\!t$, as it would hold for Markovian or $i.i.d.$ cases.

We now turn our attention to the simulation results under Slashdot and Wiki-Talk graphs. As the graphs are much bigger than the Gnutella graph, we properly increase the number of agents from 100 to 2000. We observe the similar trends as in the Gnutella graph and only present the TVD results for the cases with $w_{k}\!=\!k^{\alpha}$ ($\alpha\!=\!0,1,3,5,10$) when using the dynamic NMMC method with $c_{t}$-updating probability $p\!=\!0.01$ in Figures 4 and 5, which consistently outperforms the other choices of $p\!=\!1,0.1$ for all three cases of the target QSD. We still observe a polynomial decrease in the convergence speed of $\hat{\bm{\mu}}_{t}$, which is asymptotically linear on a log-log scale. We see that the case with $\alpha\!=\!3$ or $\alpha\!=\!5$ exhibits the fastest decaying rate in its TVD curve for achieving $\bm{\pi}\!=\!{\mathbf{u}}$ under both Slashdot and Wiki-Talk graphs and for estimating the EVC $\mathbf{x}$ under Slashdot graph. On the other hand, the case with $\alpha\!=\!1$ becomes the best with the fastest decaying rate for the other cases. This again shows that to strike a balance between diffusion and redistribution, the extent of diffusion, which can be controlled by the choice of $\alpha$, should be properly chosen based on the underlying graph structure and the target QSD. We also note that the convergence becomes much faster when the target QSD is set to be the EVC $\mathbf{x}$ of $\mathbf{A}$ or the in-degree distribution, i.e., $\pi(j)\!\propto\!d^{-}_{j}$, as compared to the case of the uniform distribution $\bm{\pi}\!=\!{\mathbf{u}}$.

We finally report the TVD results of the dynamic NMMC method with $c_{t}$-updating probability $p\!=\!0.01$ under Amazon graph for the cases with $w_{k}\!=\!k^{\alpha}$ ($\alpha\!=\!0,1,3,5,10$) when achieving $\pi(i)\propto d^{-}_{i}$. Since the size of the LSCC of the Amazon graph is over $\times 100$ bigger than that of the Gnutella graph, we use $10^{4}$ agents. We again observe the polynomial decay in the convergence speed of $\hat{\bm{\mu}}_{t}$. We also see that the case with $w_{k}\!=\!1$ or $w_{k}\!=\!k$ leads to the fastest decaying rate. This suggests that too strong a bias in weights toward more recently-visited nodes can have an adversary effect, leading to an imbalance between diffusion and redistribution.

6. Discussion and Practical Considerations

In this section, we thoroughly discuss the comparison between our NMMC method and the relevant literature, ranging from the MCMC methods to the PageRank algorithm. We also address practical issues with our NMMC method and demonstrate that our NMMC method is readily applicable to most direct graphs in practice.

6.1. MCMC vs. NMMC

We first discuss the key differences between the popular MCMC method and our NMMC method for sampling from a given probability distribution $\bm{\pi}$ on a graph $\mathcal{G}$. In essence, the MCMC method is about how to build a Markov chain $\mathbf{P}$ on a given state space (e.g., graph) such that it attains any desired stationary distribution $\bm{\pi}$, i.e., $\bm{\pi}=\bm{\pi}\mathbf{P}$ on the same space (Liu, 2004; Richey, 2010). The Metropolis-Hastings (MH) algorithm, undeniably the most famous one in the MCMC literature, can be described as follows. For a given target distribution $\bm{\pi}$ to achieve on a graph, consider an underlying Markov chain $X_{t}$ on this graph with its transition matrix $\mathbf{Q}=[Q_{ij}]$. At the current state $X_{t}=i$, the chain proposes to move to $X_{t+1}=j$ with probability $Q_{ij}$, which is then accepted with probability

[TABLE]

and rejected with probability $1-\gamma_{ij}$ in which case the chain stays at the same state, i.e., $X_{t+1}=i$. Clearly, it then immediately follows that

[TABLE]

for all $i,j\in\mathcal{N}$, i.e., the resulting stochastic matrix $\mathbf{P}=[P_{ij}]$ with $P_{ij}=\gamma_{ij}Q_{ij}$ for $i\neq j$ and $P_{ii}=1-\sum_{k\neq i}P_{ik}$ becomes reversible with respect to $\bm{\pi}$, thus satisfying the original equation $\bm{\pi}=\bm{\pi}\mathbf{P}$. This method has proven to be extremely versatile and powerful, using only local information to find $\gamma_{ij}$ with any transition matrix $Q_{ij}$ on the given state space. From an algebraic point of view, this amounts to constructing a reversible stochastic matrix $\mathbf{P}$ such that its leading left-eigenvector matches a given $\bm{\pi}$. As mentioned in Section 1, again, the most important yet implicit assumption in (25) is that $Q_{ji}>0$ whenever $Q_{ij}>0$, i.e., nodes $i$ and $j$ have reciprocal relationship. If not, i.e., there exists a directed edge from $i$ to $j$ but not vice versa, we would have $Q_{ji}=0$ even if $Q_{ij}>0$, for which the MH algorithm no longer applies. What we have established in this paper is to show how we can remove the requirements of (i) being stochastic and (ii) reciprocity in the underlying process.

One may ask what if our NMMC method is applied to an ‘undirected’ graph and/or how it would be compared to the MH algorithm under the undirected graph for graph sampling. We again emphasize that our NMMC method is the first of its kind, applicable to directed graphs or in general when the proposed move is non-reciprocal. Our NMMC method is not meant to be compared with the MH algorithm, as the MH cannot even be applied to directed graphs and our algorithm is not an improvement over the MH algorithm.

Nonetheless, for the sake of completeness, we provide simulation results of the MH algorithm and our NMMC method when applied to an undirected graph. To this end, we use the undirected version of the Gnutella graph (its strongly connected component). For the MH algorithm, we consider the following two popular random-walk versions, since it still depends on the choice of $\mathbf{Q}$. One is to use $Q_{ij}\!=\!1/d_{\max}$ if $(i,j)\!\in\!\mathcal{E}$ and $Q_{ij}\!=\!0$ if $(i,j)\!\not\in\!\mathcal{E}$ and $i\!\neq\!j$, with $Q_{ii}\!=\!1-d_{i}/d_{\max}$, where $d_{i}$ is the degree of node $i$ and $d_{\max}$ is the maximum degree of the undirected graph (Bar-Yossef et al., 2000; Boyd et al., 2004). We call this version MH-MAX. The other is to use the simple random walk for $\mathbf{Q}$, where $Q_{ij}\!=\!1/d_{i}$ if $(i,j)\in\mathcal{E}$ and $Q_{ij}\!=\!0$, otherwise (Stutzbach et al., 2009; Gjoka et al., 2011; Lee et al., 2012). We call this version MH-SRW. Note that the transition matrix of the simple random walk is still used as the proposed chain $\mathbf{Q}$ under our NMMC method, but the resulting acceptance probability in (10) is now given by $\gamma_{ij}\propto\frac{\pi(j)}{\pi(i)}\frac{d_{i}}{d_{j}}$, since the underlying graph is undirected.

We here consider $\bm{\pi}\!=\!{\mathbf{u}}$ and use 100 random walks. Letting $\{X_{t}\}$ be the trajectory by the MH-MAX (or MH-SRW), we consider $\hat{x}(i)\!:=\!\frac{1}{t}\sum_{k=1}^{t}\mathbf{1}_{\{X_{k}=i\}}$ to estimate $\pi(i)\!=\!1/n$. For performance metrics, we use the TVD and also the normalized root mean square error (NRMSE), which is defined as $\sqrt{\mathbb{E}\{(\hat{x}(i)\!-\!1/n)^{2}\}}/(1/n)$ for each $i$. We report the average of the NRMSE over all $i\!\in\!\mathcal{N}$. Figure 7 shows the TVD and NRMSE results (on a log-log scale) of MH-MAX, MH-SRW, and our NMMC method with $w_{k}\!=\!k^{\alpha}$ ($\alpha\!=\!1,3,5,10$) and $p\!=\!0.01$. We see that all the TVD and NRMSE values decay to zero polynomially fast as we discussed in Section 5.

6.2. PageRank vs. NMMC

Our NMMC method may be reminiscent of the famous PageRank algorithm and its associated ‘random surfer’ model due to the seemingly common nature of following successive outgoing links at random and occasionally performing a random jump to a node in the graph. They are, however, fundamentally different, as explained below.

The random surfer model was introduced in the seminal work by Brin and Page (Brin and Page, 1998; Page et al., 1999) in order to interpret the PageRank vector as the stationary distribution, say $\bm{\pi^{*}}$, of an ‘imaginary’ random walk (or random surfer) on a Web graph. Namely, a random surfer keeps clicking successive outgoing links uniformly at random and occasionally gets bored of following links then jumps or teleports to a random page (by entering a new destination in the browser’s URL line). In the long run, the proportion of the time the random surfer spends on a given page becomes the PageRank score of the page. We again emphasize that the random surfer model was just a random-walk “interpretation” of the PageRank algorithm, not an actual algorithmic operation implementable by actual crawlers on the graph. The rationale behind such a random-walk interpretation is to enable one to resort to the results of random walks on a graph, or the theory of Markov chains, in order to understand the properties of the PageRank algorithm such as convergence properties as in the seminal work (Brin and Page, 1998; Page et al., 1999) and follow-up studies (Haveliwala and Kamvar, 2003; Langville and Meyer, 2006).

However, we are not concerned with speeding up the PageRank algorithm nor devising an efficient numerical method for computing $\bm{\pi^{*}}$. We are rather focused on a somewhat reverse problem, which is to obtain samples from a desired probability distribution $\bm{\pi}$ over a directed graph without requiring the global knowledge of the graph structure. It is also worth noting that unlike the operation of a random surfer/walk in the random surfer model that is to teleport to any random node (even including not-yet-visited ones), our NMMC method is only limited to making a random jump (as a process of redistribution) to one of already-visited nodes, which can be readily implemented in practice without acquiring the global knowledge. We do not assume a full access to the state space, which would be necessary to realize a teleportation in the random surfer model if one wishes to jump to any arbitrary page (associated with its own URL) in a Web graph, and would be equivalent to assuming that one can draw $i.i.d.$ samples directly from the space.

6.3. Online In-degree Estimation

We validate the feasibility of our NMMC method with an online estimate of the in-degree of each node, which is possibly latent information. As explained in 3.2, the in-degree information is often readily available as part of a user profile in online social networks such as Google++, Twitter, and Instagram. However, in other types of directed graphs like Web graphs, the in-degree information may not be available but needs to be discovered. To make our NMMC method applicable to the latter, we use a simple way of estimating the in-degree of node $i$. Each crawling agent, in the stage of proposing a movement from $i$ to $j$, can discover the presence of an incoming edge from $i$ into $j$. Whenever a new incoming edge is discovered, the agent increments the value of the in-degree of $j$, say $\hat{d}^{-}_{j}$, with $\hat{d}^{-}_{j}\!:=\!1$ initially. All the agents share the estimate of the in-degree of each node. We report the resulting performance of our NMMC method in Figure 8 and observe that our NMMC method performs well even with the simple online in-degree estimation. This manifests the robustness and versatility of our NMMC method in practice.

6.4. Going Beyond Strongly Connected Directed Graphs

We have assumed that the underlying graph $\mathcal{G}$ is strongly connected so far. While this assumption has been made mainly for ease of exposition, it may not hold in practice. We thus discuss how we can relax this assumption so that our NMMC method becomes ready for practical use. We first note that the key assumption for the NMMC method is essentially the irreducibility of the ‘proposed’ Markov chain $\mathbf{Q}$ on transient states $\mathcal{N}$, as explained in Section 3. This may not be equivalent to assuming that the underlying directed graph be strongly connected for sampling from an arbitrary probability distribution $\bm{\pi}$. It is worth noting that for the estimation of the EVC $\mathbf{x}$ of the directed graph, the underlying graph needs to be strongly connected by definition, or to properly define the eigenvector centrality as a valid measure (Meyer, 2000). In this case, the assumption of strongly connected graphs is not a limitation at all. Thus, we below focus only on achieving $\bm{\pi}$.

If the directed graph $\mathcal{G}$ is not strongly connected, the support of $\bm{\pi}$ would need to be on the set of all the nodes, say $\mathcal{N}_{S}$, that are “reachable” by traversing directional edges when starting from an initial node, or a set of initial nodes, say $S$, that are reachable among themselves. Note that the “unreachable” nodes $\mathcal{N}\setminus\mathcal{N}_{S}$ are simply invisible and can safely be ignored. The technical changes are then how to deal with (i) ‘dangling’ nodes, which are the nodes with no forward links, and (ii) a set of nodes that are reachable but without outgoing edges to the other nodes in $\mathcal{N}_{S}$. Note that these two problems were resolved in the PageRank algorithm by adopting the teleportation to any node in the graph, yet to numerically compute the (unknown) PageRank scores. While such global teleportation is infeasible for sampling in practice, its personalized PageRank version is practically implementable and can be adopted into our NMMC method.777The personalized PageRank scores are still unknown and need to be computed by the PageRank algorithm (or the power method), as long as the PageRank scores themselves are concerned.

Fix the set of initial nodes $S$. The state space $\mathcal{N}$ then reduces to $\mathcal{N}_{S}$. For the proposed chain $\mathbf{Q}$ on $\mathcal{N}_{S}$, we consider the following transition matrix: For $i,j\in\mathcal{N}_{S}$,

[TABLE]

for some constant $p\!\in\!(0,1)$. At each node $i$, with probability $p$, the agent follows one of its outgoing edges (if exists) uniformly at random, and jumps to one of the nodes in $S$ uniformly with probability $1-p$. If node $i$ is a dangling node with no outgoing edges ($d_{i}^{+}\!=\!0$), then the agent always jumps to a random node in $S$. It is not difficult to see the irreducibility of $\mathbf{Q}$ on $\mathcal{N}_{S}$ for our NMMC method. Then, the acceptance probability $\gamma_{ij}$ in (9) can properly be obtained. In Figure 9, we present the TVD results for estimating the uniform ${\mathbf{u}}$ on $\mathcal{N}_{S}$ using the dynamic NMMC method with $\mathbf{Q}$ in (26), where $\mathcal{N}_{S}$ is the set of nodes that are reachable from any node in the LSCC of Gnutella graph and $|\mathcal{N}_{S}|\!=\!8,566$. We here use 300 independent crawling agents starting from 300 fixed, distinct nodes $S$ and $p\!=\!0.95$. Figure 9 demonstrates the feasibility of our NMMC method even when the underlying graph is not strongly connected.

6.5. Query Cost vs. Performance for Graph Sampling

So far, all the simulation results have been shown with respect to the number of time steps, since they have the primary purpose of numerically validating the convergence of our NMMC method and evaluating its convergence speed. Recall that the ergodic theorem in (12) and the CLT in (24) hold for ergodic Markov chains with respect to the number of time steps. Nonetheless, when it comes to graph sampling via crawling, ‘query cost’ becomes practically important. In general, crawling real-world networks is done through public yet restrictive API interfaces provided by the networks or their Web interfaces (by following URLs and hyperlink structures). Such external access is often rate-limited due to security concerns, privacy protection, infrastructure burden and commercial purposes (Vesdapunt and Garcia-Molina, 2016; Xu et al., 2017). For example, Twitter allows only 15 API requests to retrieve IDs of a user’s followers or IDs of the friends that the specific user is following every 15 minutes (Twitter, [n. d.]).

From a practical point of view, it would be informative to see the performance of our NMMC method with respect to the query cost in terms of the API calls (or HTTP requests) required. We below present the performance of NMMC method with respect to the query cost under Gnutella and Slashdot graphs. To this end, we adopt the definition of the query cost commonly used in the literature (Gjoka et al., 2011; Ribeiro et al., 2012; Zhou et al., 2015, 2016), which is the number of unique queries required. In other words, API calls are spent only when visiting a node (or querying its page/profile) for the first time, and revisiting the same node since then (duplicate queries) can be done immediately by retrieving the nodal information from local cache without consuming the query rate limit. This indicates that the redistribution mechanism of our NMMC method is achieved at no additional cost as long as query cost is concerned. Thus, the ‘due’ redistributions are only advantageous to shape up the historical empirical distribution $\hat{\bm{\mu}}_{t}$ toward the target QSD. The local cache can also be easily shared among multiple crawling agents.

In Figures 10 and 11, we present the TVD results of the dynamic NMMC method with $k^{\alpha}$ ($\alpha\!=\!0,1,3,5,10$) and $p\!=\!0.01$ under Gnutella and Slashdot graphs, respectively, as the query cost varies. Here, the query cost is defined as the number of unique queries normalized by the graph size (LSCC size) and given in percentage. The same numbers of agents are used as before for Figures 2 and 4, respectively. We see that the choice of $w_{k}\!=\!k$ exhibits the overall fastest decaying rate in its TVD curve for all three cases of the target QSD under both Gnutella and Slashdot graphs. This also indicates that higher redistribution weights toward more recently-visited nodes is no longer beneficial in improving the convergence speed when the query cost is concerned.

In a similar vein, we showed the TVD results of our NMMC method until $10^{4}$ time steps in Figure 2, which is greater than the LSCC size of the Gnutella graph. Considering 100 agents used, one might expect that the entire graph is already visited by that time. However, it has not been the case, since there are a lot of redistributions that are counted toward the number of time steps, but not toward the number of unique queries. The multiple agents can also visit the already-visited nodes many times, which is again not counted toward the number of unique queries. This is confirmed in Figure 12, which shows the number of unique queries, which is normalized by the LSCC size, for all three cases of the target QSD on the Gnutella graph using the dynamic NMMC method with $k^{\alpha}$ ($\alpha\!=\!0,1,3,5,10$) and $p\!=\!0.01$.

6.6. Comparison with the Latest Approach for Graph Sampling

Lastly, we provide a comparison with the latest approach in (Ribeiro et al., 2012) for sampling directed graphs. Among only a few studies available in the literature on sampling from directed graphs, which shall be explained later in Section 7, the latest practical solution is the one in (Ribeiro et al., 2012), which is called Directed Unbiased Random Walk (DURW) algorithm. While our NMMC method is capable of achieving any distribution $\bm{\pi}$ on a directed graph and estimating its EVC $\mathbf{x}$, which is much more than just achieving ‘uniform’ samples from the directed graph, the comparison of our NMMC method to the DURW algorithm would be interesting in the context of uniform graph sampling and is provided below. We here focus on estimating $\bm{\pi}\!=\!{\mathbf{u}}$.

The DURW algorithm has two main parts. The first part is to transform the underlying directed graph to an undirected graph. At each newly visited node, it keeps outgoing edges only to ‘not-yet-visited’ nodes and converts them as undirected edges. This edge-pruning idea was originally proposed in (Bar-Yossef et al., 2000). Upon the transformed, undirected graph, the second part is to use (i) a weighted random walk with random jumps to everywhere in the graph and (ii) the importance sampling. Specifically, the weighted random walk at node $i$ chooses one of its neighbors uniformly at random, with probability $\frac{d_{i}}{d_{i}+w}$, and jumps (or teleports) to any random one in the graph, with probability $\frac{w}{d_{i}+w}$, where $w$ is the weight parameter for choosing a random jump and $d_{i}$ is the degree of $i$. While the former operation consumes a unit cost if a new node is sampled, the latter consumes the cost $c\!>\!1$. For a directed social network, assuming that its user-ID space is accessible, the cost $c$ of a random jump is the number of API calls required on average to find a valid ID for sampling. It clearly depends on how sparse the user-ID space may be, which is unknown and uncontrollable. Note that the random jumps may not be feasible at all in practice, e.g., for sampling Web graphs. In addition, letting $\bm{\pi^{\prime}}$ be the stationary distribution of the weighted random walk with random jumps, one can see that $\pi^{\prime}(i)$ is proportional to $d_{i}\!+\!w$. The importance sampling is then applied to correct the bias from the non-uniform $\bm{\pi^{\prime}}$. Thus, letting $\{X_{t}\}$ be the trajectory by the DURW algorithm, it is to use the ratio estimator $\hat{x}(i)\!:=\!\frac{\sum_{k=1}^{t}\mathbf{1}_{\{X_{k}=i\}}/\pi^{\prime}(X_{k})}{\sum_{k=1}^{t}1/\pi^{\prime}(X_{k})}$ to estimate $\bm{\pi}\!=\!{\mathbf{u}}$. Note that, instead of $\bm{\pi^{\prime}}$, its unnormalized version can be used in practice. We refer to (Ribeiro et al., 2012) for more details.

Figure 13 shows the TVD and NRMSE results of our NMMC method with $w_{k}\!=\!k^{\alpha}$ ($\alpha\!=\!0,1,3,5,10$) and $p\!=\!0.01$, and the DURW algorithm with various choices of $w$ and $c$ under Gnutella graph, which are from (Ribeiro et al., 2012). For example, the choice of $w\!=\!0.1$ behaves more like the pure simple random walk, while the larger values of $w$ put more weights on random jumps. In addition, $c\!=\!10$ means 10 API calls spent on average to sample a valid ID (or to teleport to a random node). Similarly for $c\!=\!77$. In addition, we here consider a single random crawler. It is because in the DURW algorithm, each crawler has to build its own state space (or the transformed, undirected graph) as it explores the graph, which is clearly time-varying and depends on its sample path, hence the state space can never be shared with other crawlers. The query cost is again defined as the number of unique queries normalized by the graph size (LSCC size) and given in percentage. We report the average of the NRMSE over all $i\!\in\!\mathcal{N}$.

We see that our NMMC method generally performs better than the DURW algorithm. A possible explanation would be that the way to transform the underlying directed graph to an undirected graph makes the transformed one ‘tree-like’. Thus, traversing the transformed graph could get slower. While random jumps can speed up the graph traversal, the high cost of each random jump eventually outweighs the benefit of moving to everywhere in the graph. We also observe that carefully choosing the value of $w$ with given (yet unknown) cost $c$, the DURW algorithm would be comparable to our NMMC method. In other words, the DURW algorithm requires a good choice of $w$ a priori against unknown $c$, before running the random crawler to explore a directed graph by ‘undirectifying’ directed edges on the fly. In contrast, our NMMC method performs well for a wide range of values of $\alpha$ for $w_{k}\!=\!k^{\alpha}$.

We finally emphasize that while we take (directed) graph sampling as a primary application of our NMMC method in this paper, the NMMC method can potentially find much broader scope of applications beyond the graph sampling as we have witnessed the widespread applications of the MH algorithm in many disciplines. The potential broader impact stems from the versatile applicability of the NMMC method for generating any arbitrary probability distribution over the state spaces violating the reciprocity condition, where the MH algorithm and other standard MCMC methods are no longer applicable.

7. Related Work

There have been a few studies on achieving uniform samples from directed graphs, ranging from traditional Web graphs (Bar-Yossef et al., 2000; Henzinger et al., 2000; Rusmevichientong et al., 2001) to recent online social networks (Wang et al., 2010; Ribeiro et al., 2012; Gabielkov et al., 2014). The key underlying methods have been mainly (explicitly or implicitly) the MH algorithm and importance sampling for Markov chains as in (22). The ‘non-uniform’ random walks, along with the importance sampling, are typically the simple random walks and (personalized) PageRank walks (Henzinger et al., 2000; Rusmevichientong et al., 2001). Unfortunately, however, their stationary distributions are simply unknown, and the only viable approach in (Henzinger et al., 2000; Rusmevichientong et al., 2001) has been to empirically “estimate” the stationary distribution (e.g., based on the visit counts obtained by running a separate random walk), which is still known to suffer from the bias toward the nodes with high in-degree (Bar-Yossef and Gurevich, 2008). Another simple way to get around the difficulty is to just convert every directed edge to an undirected one (excluding the duplicate ones) (Rusmevichientong et al., 2001; Wang et al., 2010; Gabielkov et al., 2014), which technically has noting to do with sampling from directed graphs.

A more practical method was proposed in (Bar-Yossef et al., 2000) to construct an undirected graph on-the-fly from a target directed graph for sampling, only considering the inbound links for each node discovered before the node is visited for the first time, and then to apply, in essence, the MH algorithm with $Q_{ij}\!=\!1/\max_{k}d^{+}_{k}$ over the transformed undirected graph to obtain uniform samples. In short, the edge conversion is done with edge pruning. This idea has been adopted later in (Ribeiro et al., 2012) for a weighted random walk with ‘random jumps’ (to anywhere in the undirected graph), equipped with the importance sampling due to its non-uniform stationary distribution (on the undirected graph). In contrast to all these, without knowing the stationary distributions or requiring any transformation of the directed graph to another undirected version, our NMMC method is able to sample from any target distribution over the directed graph with theoretical guarantees. We recently came to know the work by Avrachenkov et al. (Avrachenkov et al., 2018) in which they briefly mention that their proposed reinforcement learning-based sampling technique can potentially be applied to strongly connected directed graphs, but without any theoretical and practical demonstrations. We here emphasize that our NMMC method is even capable of sampling from a given distribution $\bm{\pi}$ over the set of reachable nodes $\mathcal{N}_{S}$ of the directed graph, even if the graph is not strongly connected, not to mention the capability of estimating the EVC $\mathbf{x}$ of the (strongly connected) directed graph.

8. Conclusion

In this paper, we have shown how to go beyond the limitation set by the current MCMC methods when state spaces violate the reciprocity condition and the global knowledge is unavailable. Unlike the most studies in the literature that aim to improve a given MCMC algorithm or to propose a better one within the same framework, our approach is fundamentally different and unique, in that it works directly over any nonreversible proposed chain such as simple random walks on directed graphs without sacrificing their distributed nature in the implementation. Our NMMC framework builds upon our careful mapping from the proposed chain to the target distribution on a set of transient states, and entails the machinery of the quasi-stationary distribution of a suitably constructed transient chain and the induced random walks with reinforcement to relocate to positions from their past history, which are amenable to distributed implementation using only locally available information. We have also devised a number of ways to improve the overall speed of convergence and to enhance the practicability, by utilizing online estimates of in-degree and the required global constant and by relaxing the assumption of the strong connectedness of the underlying directed graph, to make our NMMC method readily applicable to most directed graphs in reality.

9. Acknowledgments

The authors thank the anonymous reviewers and especially our shepherd, Konstantin Avrachenkov, for many constructive comments and suggestions that greatly improved the quality of this paper. Do Young Eun is supported in part by National Science Foundation under Grant Nos. CNS-1423151 and CNS-1824518.