Self-Repellent Random Walks on General Graphs - Achieving Minimal Sampling Variance via Nonlinear Markov Chains (Extended Abstract)

1 Introduction

Random walk-based techniques are a staple in statistics and learning theory. Markov chains such as the Metropolis Hastings random walk, designed to achieve any given target probability distribution as its stationary measure, are widely used as Markov chain Monte Carlo (MCMC) samplers and in distributed optimization via stochastic gradient descent [Sun et al., 2018; Hu et al., 2022]. The local nature of the information required to compute state transition probabilities means that the algorithms scale well and are robustly implementable over state spaces such as large graphs/networks with general topologies. However, classic Markov chains can often be victims of limitations set by the underlying topology of the state space (communication matrix or adjacency matrix of the underlying network structure) leading to correlated samples which can negatively affect the estimator performance. It has also been well established that the time-reversibility requirement for the classical MCMC samplers is one of the causes for their slow convergence (see Section 1 [Andrieu and Livingstone, 2021]). One way in which this problem has been approached in the literature is via construction of non-reversible versions of the base Markov chain [Diaconis et al., 2000; Turitsyn et al., 2011; Chen and Hwang, 2013; Ma et al., 2016; Thin et al., 2020], which is often done by inducing some form of non-backtracking behaviour, that is, avoiding states most recently visited by the random walker [Alon et al., 2007]. This involves the random walker interacting with some of its own past history, and has been shown to possess better efficiency than the original base Markov chain in the sense of the MCMC estimator achieving a smaller asymptotic variance [Neal, 2004; Lee et al., 2012]. Since these non-backtracking based methods only utilize the most recent history of the random walker and are still provably more efficient, it is natural to consider the design of protocols where the random walker interacts with its entire past history to speed up its diffusion and increase its sampling efficiency, especially for sampling over discrete state spaces. This is the approach taken in our paper.

Let G(N, E) be an undirected, connected graph where N ≜ {1, · · · ,N} denotes the set of nodes and E denotes the set of edges, where we say (i, j) ∈ E if there is an edge between nodes i, j ∈ N. We use A = [aij]i,j∈N to represent the adjacency matrix of the graph, where aij > 0 if

(i, j) ∈ E, and zero otherwise; N(i) ≜ {j ∈ N | (i, j) ∈ E} refers to the set of neighbors of node i; deg(i) ≜ ∑j∈N aij will refer to the degree of each node i ∈ N. Denote by Σ the N-dimensional probability simplex over N, with Int(Σ) denoting its interior, and let P ≜ [Pij]i,j∈N be the transition probability matrix of an ergodic, time-reversible Markov chain over N, with its stationary distribution µ ≜ [µi]i∈N. Without loss of generality, we assume Pij > 0 if and only if aij > 0. In this setup, we design Self-Repellent Random Walks (SRRWs) on general graphs 1 indexed by a tunable parameter α ≥ 0, all of which can sample from µ ∈ Σ, and then study their sampling 'efficiency' as a function of α (with α = 0 being equivalent to the baseline Markov chain with transition kernel P).

The SRRW transition kernel: Consider the Markov chain kernel (transition matrix) K[x] ≜ [Kij [x]]i,j∈N , whose transition probabilities are mappings Kij : Σ → [0, 1], given by

for any probability vector x ≜ [xi]i∈N ∈ Σ. Here, {rµi }i∈N is a family of positive functions rµi : [0, 1] → R+ parameterized by µi, with rµi (xi) decreasing in xi ∈ [0, 1] and rµi (µi) = C, for all i ∈ N. 2 Transition probability kernels defined in this fashion, taking probability distributions as argument, are called 'nonlinear' Markov kernels [Andrieu et al., 2007; Andrieu et al., 2011], as opposed to classical Markov chains with kernels P that are often interpreted as linear operators - the transition probabilities at each step being independent of x (i.e., the case where rµi (·) is a constant function).

Stochastic processes utilizing nonlinear Markov kernels are called nonlinear Markov chains, and can be simulated/generated using self-interacting Markov chains (SIMCs) (see [Del Moral and Miclo, 2004; Del Moral and Miclo, 2006; Moral and Doucet, 2010]. Let {Xn}n≥0 be a random walker over N, and let xn be its occupational measure or historical empirical distribution up to time n ≥ 0, written as

where δXk is the delta measure whose Xk'th entry is one and the rest are zero, thus recording the position of the random walker at time k ≥ 0. The process {Xn}n≥0 becomes a SIMC if at each time step n ≥ 0, the random walker makes transitions according to some nonlinear kernel K[xn], not necessarily as defined in (1). We say that the process {Xn}n≥0 is a SRRW if it is a SIMC with K[x] as in (1). We use the term self-repellent since at each time step, the transition probability to a node j ∈ N is proportional to rµj ([xn]j) where [xn]j = 1n ∑nk=1 1{Xk=j}, and is thus a decreasing function of the visit count to j ∈ N. In other words, the walker is less likely to move to a node that has been visited more often so far (thus self-repellent), as shown in Figure 1.

1 We consider graphs because they represent a generalization of (discrete) finite state spaces by imposing a communication (adjacency) matrix. The existence of an edge between two nodes (states) represents a non-zero probability of state transitions between the two nodes.

2 As we shall see later, xi will be directly proportional to the visit count to any node ∈ N, since x ∈ Σ will be the empirical distribution of the self-repellent randomwalk.

When Pij ∝ aij in (1) for each i ∈ N, the SRRW is a self-repellent version of the well-known simple random walk (SRW) procedure, with the target distribution being proportional to the degree of the nodes, that is, µi ∝ deg(i) for all i ∈ N. Like most general MCMC procedures, the SRRW kernel can also be defined for any given sampling distribution µ ∈ Int(Σ), for instance, by setting P to be the transition matrix of a Metropolis Hastings Random Walk (MHRW) with stationary distribution µ. For example if µi = 1/N, that is µ = 1N 1 - the uniform distribution over the set of nodes N, then we can choose Pij = min { 1 deg(i) , 1 deg(j) } for all (i, j) ∈ E, with Pii = 1 - ∑j̸=i Pij. The matrix P defined in this manner is the MHRW kernel with the uniform distribution as its stationary measure, and is among the most commonly used kernels for unbiased graph sampling [Lee et al., 2012; Li et al., 2015] and distributed optimization [Sun et al., 2018]. The elegance in the Metropolis Hastings algorithm and the key to its widespread adaptation lies in the fact that at each time step, the entries of µ need only to be known for the neighbouring nodes of the random walkers (that is, only local information required), and only up to a constant multiple. This property ensures a robust, scalable implementation of the MHRW, since global constants are often unknown for large networks a priori.

Our SRRW construction begins with rµi (·) taking a polynomial form for all i ∈ N, given by

where the parameter α ≥ 0 can be perceived as the strength of the self-repellence mechanism designed into the SRRW transition kernel. Similar to the MHRW transition kernel, only the local information regarding entries of µ needs to be known at any given time step, and up to only a constant multiple. For convenience, we formalize this property as scaleinvariance (S.I.): an SRRW kernel possesses S.I. if for all i, j ∈ N

• (i) Computing Kij[x] only requires knowing µk for k ∈ N(i), and only up to a constant multiple for any i ∈ N;

Indeed, we show in [Doshi et al., 2023, Appendix D] that out of all possible forms for the functions rµi (·), only the polynomial form as in (3) possesses the S.I. property. Henceforth, we restrict ourselves to the polynomial form of rµi (·).

2 Our Contributions

2.1 Almost Sure Convergence and CLT

We show that given any MCMC kernel P which samples from a target distribution µ, the corresponding SRRW is asymptotically more efficient as a random walk-based sampler. We do this by first showing that

We then provide second-order convergence results in the form of a central limit theorem (CLT); that is, we show that there exists an asymptotic co-variance matrix V(α) ∈ RN×N parameterized by α ≥ 0, such that

We obtain these results by first viewing the SRRW as a stochastic approximation (SA) algorithm with statedependent noise [Harold J. Kushner, 1997; Fort, 2015], allowing us to form a connection between the stochastic process and a deterministic system of ordinary differential equations (ODEs). We establish global convergence results for this ODE system, which lay the foundations for proving the almost sure convergence of the stochastic SRRW process.

2.2 Co-Variance Ordering in Parameter α

For any α ≥ 0 we derive the exact form of V(α) in terms of α and the spectrum (eigenvalues and eigenvectors) of P. This allows us to show that kernels parameterized by larger α are asymptotically more efficient samplers; that is, they achieve smaller sampling variance. This is done by showing that the asymptotic covariance matrices follow a Loewner ordering: 3

In other words, as long as the numerical/computational stability of the random walk implementation can be ensured, larger values of α are always more favourable in terms of achieving a smaller (asymptotic) sampling variance. In [Doshi et al., 2023, Corollary 4.7], we also derive an upper bound on the ratio of its sampling variance to that of the baseline Markov chain, and show that this upper bound goes down to zero as α → ∞ with speed O(1/α). This is surprising because asymptotically for large enough α, the SRRW, which is a stochastic process whose trajectories are constrained by 'walking' on the underlying communication matrix of the network, achieves smaller sampling variance than an i.i.d. sampler 4 whose variance is always a constant positive value.

2.3 Application to Distributed Stochastic Optimization

We study a family of distributed stochastic optimization algorithms, known as token algorithms, where gradients are sampled by a token traversing a network of agents in random-walk fashion [Sun et al., 2018; Hu et al., 2022; Even, 2023]. Typically, these random-walks are chosen to be Markov chains that asymptotically sample from the desired target distribution µ, and play a critical role in the convergence of the optimization iterates. Our paper [Hu et al., 2024], as illustrated in Figure 2, takes a novel approach by replacing the standard linear Markovian token by SRRW - a nonlinear Markov chain. The SRRW-driven token algorithm is described as follows.

3 Matrices A, B follow the Loewner ordering A <L B if A ̸= B and B - A is positive semi-definite.

4 This corresponds to a sampler that can visit any node i with probability µi independent of its previous position at any given time. Clearly, in the graph setting, this requires the sampler to 'jump' to any other node by ignoring the underlying network structure altogether - something which random walkers on general graphs are not permitted to do.

where {βn} and {γn} are step size sequences decreasing to zero, and K[x] is the SRRW kernel (1) with rµi as in (3). Iterations (7a) and (7b) (i.e. without considering (7c)) correspond to a generalized SRRW process, with now xn being a weighted empirical measure. Setting γn = 1/(n + 1) coincides with the unweighted case discussed earlier. The H(θn,Xn) terms in the update rule (7c) driving the optimization iterates {θn} embeds gradient information such that solving EX∼µ[H(θ∗ , X)] = 0 is equivalent to obtaining a local minimizer for the distributed optimization problem. This is because (7c) corresponds to the more general stochastic approximation (SA) iteration, and numerous stochastic optimization algorithms which can be expressed as token algorithms, such as stochastic gradient descent (SGD) [Hu et al., 2022] and stochastic heavy ball (SHB) [Gadat et al., 2018; Li et al., 2022], can be expressed as a special case of SA. Thus, since the update rule (7) is essentially a SA algorithm with SRRW driven noise sequence, we call it SA-SRRW. The detailed algorithmic setup and model assumptions on SASRRW can be found in [Hu et al., 2024, Sections 1 and 2].

The update rule (7) is in fact an example of a two-timescale stochastic approximation with state dependent (controlled) Markov noise. For any α ≥ 0, we prove almost sure convergence of the optimization iterates θn, as well as a CLT. That is, we show that

where Vθ(α) is the asymptotic co-variance matrix for the iterates {θn}.

Our key result involves the case when βn = o(γn), for which we say that θn is on a slower timescale compared to xn. In this case, in addition to proving a co-variance ordering as in Section 2.2, i.e.,

we also prove that Vθ(α) decreases to zero at a rate of O(1/α2 ) [Hu et al., 2024, Proposition 3.4]. This is especially surprising, since we achieve an amplification in performance gain when using the SRRW for distributed stochastic optimization, as compared to the Monte Carlo sampling application where the performance gain in α is in O(1/α).

3 Conclusion

Our work [Doshi et al., 2023] introduces the SRRW as a dropin replacement for any Markov chain with stationary distribution µ, and shows that it achieves better asymptotic performance for Monte Carlo sampling tasks as the degree of repellency α increases. Our follow-up work [Hu et al., 2024] extends the scope from sampling to distributed stochastic optimization (and stochastic approximations in general), showing that certain combinations of step-sizes also achieve accelerated performance in α. Numerically tests show that the asymptotic performance guarantee is also realized for finite time steps, as can be seen in Figure 3. Our work is an instance where the asymptotic analysis approach allows the design of better algorithms despite the usage of unconventional noise sequences such as nonlinear Markov chains like the SRRW, for which traditional finite-time analytical approaches fall short, thus advocating their wider adoption.

Acknowledgments

This work [Doshi et al., 2023] was done primarily while Vishwaraj Doshi was with the Operations Research Graduate Program, North Carolina State University. This work was supported in part by National Science Foundation under Grant Nos. CNS-2007423 and IIS-1910749.