# Maximal dispersion of adaptive random walks

Gabriele Di Bona  School of Mathematical Sciences, Queen Mary University of
London, London E1 4NS, England  Leonardo Di Gaetano  Department of Network and
Data Science, Central European University, 1100 Vienna, Austria  Vito Latora
School of Mathematical Sciences, Queen Mary University of London, London E1
4NS, England  Dipartimento di Fisica ed Astronomia, Università di Catania and
INFN, I-95123 Catania, Italy  Complexity Science Hub Vienna (CSHV), Vienna,
Austria  Francesco Coghi<EMAIL_ADDRESS>Nordita, KTH Royal Institute
of Technology and Stockholm University, Hannes Alfvéns väg 12, SE-106 91
Stockholm, Sweden

###### Abstract

Random walks are the most versatile tool to explore a complex network. These
include maximum entropy random walks (MERWs), which are maximally dispersing
and therefore play a key role as they optimize information spreading. However,
building a MERW comes at the cost of knowing beforehand the global structure
of the network to be explored. Here, we propose an adaptive random walk (ARW),
which instead maximizes dispersion by updating its transition rule on the
local information collected while exploring the network. We show how to derive
the ARW via a large-deviation representation of the MERW as a finite-time rare
event of an unbiased random walk and study its dynamics on synthetic and real-
world networks.

In the last decade we have assisted to a wave of new works studying extreme
(rare) events associated to dynamical processes evolving on complex networks
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17]. The necessity to
understand unlikely events comes from the fact that, although rare, their
appearance determines the future of the system under study, which may be
potentially catastrophic, e.g., earthquakes [1]. In this context, researchers
have focused on random walks and their load and flow properties on nodes and
edges [2, 3, 15, 17] as models of traffic in transportation [11, 16] and
communication networks [4], or on epidemic models and extinction events [6],
or again on general order-disorder [8] and percolation transitions [9, 10] to
corroborate the robustness of networks. In these settings, rare events are
often driven by internal noise, e.g., inherent randomness in a random walk,
rather than environmental changes, and their understanding could provide us
with control mechanisms to keep away from harmful scenarios [5, 6, 8, 10].

Among all dynamical processes evolving on networks, discrete-time biased and
unbiased random walks (URWs) certainly stand out as simple and insightful
models of diffusion processes on discrete topologies [18]. URWs are local
processes that hop, at each time-step, on a neighbouring node picked from a
uniform probability distribution, locally maximizing the entropy rate, or
diffusivity, over the network. On the other hand, it is often of interest to
globally maximize the spreading of a random walk over the network, such that
it assumes uniform probability distribution over all paths. A process that
achieves this has been proposed about a decade ago [19, 20], named maximum
entropy random walk (MERW), and so far used in a myriad of applications:
testing network robustness [21, 22] and navigability [23, 24, 25, 26],
detecting connectivities and communities [27, 28, 29, 30], benchmarking
optimal mass transportation strategies [31], predicting disease associations
[32], assessing molecular trapping [33] and neutral quasispecies evolution in
biology [34] to name just a few. Such a MERW, although characterized by local
movements, i.e., it jumps to neighbouring nodes, it is non-local in the choice
of the next visited node. Indeed, it needs to know the topology of the entire
network to globally maximize the entropy rate [19, 20, 35].

The dynamics and thermodynamics of the MERW have already been studied [20, 35]
and although, in certain applications, it appears that it is not always
advantageous to consider MERWs as they may localize in portions of the network
[20], these processes are known to be maximally dispersing and hitherto
constitute the best test bench for the mixing properties of other processes
[19, 35]. Although on one hand MERW are extremely useful processes, on the
other hand, as mentioned, they require global knowledge of the embedding
network. For this reason, several attempts have been made to construct optimal
random walks with only local information [35] obtaining almost maximal-entropy
random walks with step probabilities proportional to a power of the degree of
the target node.

In this manuscript, we propose an adaptive random walk (ARW) that locally
updates its stepping rule based on the structure of the explored network.
Without requiring any prior knowledge of the whole topology, the ARW
outperforms the MERW, as it is maximally dispersive on every portion of the
network visited and not only on the whole graph. Via a bridge between large
deviation theory and network science, a MERW can be seen as a rare event of an
URW. We hence exploit this to construct the ARW as a single-trajectory rare-
event sampling algorithm [36] that only makes use of local information
available to adapt itself—changing its transition probabilities step by
step—in order to best spread on the network. As mentioned, our framework
heavily builds on large deviations, not only as a theory to understand and
predict the appearance of rare events, but also as a designing tool of the
ARW. In the following we show that, while the random walks move and discover
the underlying network, the ARW has an entropy production rate closer to the
maximal one on the visited portion of the graph if compared with an URW and a
MERW. Moreover, when the network is fully explored, the ARW and the MERW have
similar mixing properties and, in the thermodynamic limit, they become the
same process.

We start by considering an URW $X=(X_{1},X_{2},\cdots,X_{n})$ on a finite
connected and undirected graph $G=(V,E)$ characterized by a set of vertices
$V$ and a set of edges $E$. The topology of the graph is encoded in the
adjacency matrix $A=\left\\{a_{ij}\right\\}$ having components $a_{ij}=1$ if
the nodes $i$ and $j$ are connected, and $a_{ij}=0$ otherwise. The URW
dynamics is determined by the stochastic transition matrix $\Pi$, having
components

$\pi_{ij}=\frac{a_{ij}}{k_{i}}\ ,$ (1)

where $k_{i}=\sum_{j\in V}a_{ij}$ is the degree of node $i$, describing the
probability of the URW to move from $X_{\ell}=i$ at time $\ell$ to
$X_{\ell+1}=j$ at time $\ell+1$. We focus on the dynamical observable

$C_{n}=\frac{1}{n}\sum_{\ell=1}^{n}\ln k_{X_{\ell}}\ ,$ (2)

which is purely additive and, as such, only depends on the visited nodes at
each time step. Because of the ergodicity of the URW, guaranteed by the
properties of the graph $G$, $C_{n}$ converges with probability $1$ to the
ergodic average

$\left\langle\ln k\right\rangle_{\rho}=\sum_{i\in V}\rho_{i}\ln k_{i}\coloneqq
c^{*}\ ,$ (3)

where $\rho=\\{\rho_{i}\\}$ is the stationary distribution of the URW $X$, in
the long-time limit $n\rightarrow\infty$.

The observable $C_{n}$ in Eq. (21) plays a pivotal role in this manuscript.
Indeed, up to boundary terms that do not influence our discussion and thanks
to Eq. (1), we can rewrite
$C_{n}=-\ln(\mathbb{P}(X_{1},X_{2},\cdots,X_{n}))/n$, i.e., as a function of
the probability of the trajectory of an URW. Taking the long-time limit of the
average over all paths of $C_{n}$ we get the so-called Kolmogorov–Sinai
entropy production rate $h$ [37, 19, 35], i.e.,

$h=\lim_{n\rightarrow\infty}\left\langle C_{n}\right\rangle$ (4)

which, for an URW, takes the form in Eq. (3), and can be interpreted as the
mean information per step generated by the random walk. Furthermore, for
general ergodic random walks (not necessarily unbiased), the form in Eq. (4)
is known [37], with abuse of notation, to be given by

$h=-\sum_{i,j}\rho_{i}\pi_{ij}\ln\pi_{ij}\ .$ (5)

It is then clear that the observable $C_{n}$ is a random variable of the
random walk process that represents the fluctuating version of $h$, viz. the
fluctuating trajectory entropy [38] (sometimes called self-process [39]).

The finite-time fluctuating nature of $C_{n}$ around $c^{*}$ is what interest
us here. A complete understanding of the fluctuations is given by the
probability density $P_{n}(c)\coloneqq P(C_{n}=c)$, which is known to have the
large deviation form

$P_{n}(c)=e^{-nI(c)+o(n)},$ (6)

where $o(n)$ denotes corrections smaller than linear in $n$. The focus thus
moves onto studying the non-negative large deviation rate function $I$ given
by the leading exponential behaviour in Eq. (6), characterized by having a
unique zero in $c^{*}$. The rate function can be calculated by means of the
so-called Gärtner–Ellis theorem [40, 39, 41], which states that $I$ is given
by the Legendre–Fenchel transform of the scaled cumulant generating function
(SCGF)

$\Psi(s)=\lim_{n\rightarrow\infty}\frac{1}{n}\ln\mathbb{E}\left[e^{nsC_{n}}\right]\
,$ (7)

as this last is differentiable for a finite graph [41]. In particular, as the
URW is an ergodic Markov process, the SCGF can be obtained as

$\Psi(s)=\ln\zeta_{s},$ (8)

where $\zeta_{s}$ is the dominant eigenvalue of the so-called tilted matrix
$\tilde{\Pi}_{s}=\left\\{\left(\tilde{\pi}_{s}\right)_{ij}\right\\}$, with
components

$\left(\tilde{\pi}_{s}\right)_{ij}=\pi_{ij}e^{s\ln
k_{i}}=\pi_{ij}k_{i}^{s}=a_{ij}k_{i}^{s-1}\ .$ (9)

Hence, the likelihood of fluctuations can be studied focusing on the SCGF
$\Psi$ rather than the rate function $I$. In Fig. 1(a) we show the form of
$\Psi$ for the observable $C_{n}$ in Eq. (21) obtained by solving the spectral
eigenvalue problem in Eq. (8) numerically. As expected, the function looks
differentiable and since $\Psi(0)=0$ the normalisation of $P_{n}(c)$ is
guaranteed.

Figure 1: (a) SCGF $\Psi(s)$ of the observable $C_{n}$ in Eq. (21). (b)
Entropy production rate $h(s)$ of the driven process in Eq. (10); this shows a
maximum for $s=1$.

However, calculating the likelihood of fluctuations is only a first step
towards the prediction and control of rare events. It is indeed important to
also understand how these extreme events are created in time. In this context,
we construct the driven process [42, 43, 44] associated with a given
fluctuation $C_{n}=c$ that here is a locally-biased version of the URW [43,
12] whose transition probability matrix is given by

$\left(\pi_{s}\right)_{ij}=\frac{\left(\tilde{\pi}_{s}\right)_{ij}r_{s}(j)}{r_{s}(i)\zeta_{s}}=\frac{a_{ij}k_{i}^{s-1}r_{s}(j)}{r_{s}(i)\zeta_{s}}\
,$ (10)

where $r_{s}$ is the right eigenvector associated with $\zeta_{s}$. The driven
process is still Markovian and ergodic and can be interpreted as the effective
dynamics of the subset of paths of the URW leading to a fluctuation $C_{n}=c$
[43, 12, 13]; to match such a fluctuation, the Laplace parameter $s$ must
satisfy the relation $s=I^{\prime}(c)$, or

$c=\Psi^{\prime}(s),$ (11)

see [45]. In equilibrium statistical mechanics, this is similar to fixing the
equilibrium energy of the canonical ensemble by fixing its inverse temperature
(played here by the parameter $s$) to the derivative of the entropy [42].
Eventually, the entropy rate of the driven process can be obtained taking Eq.
(10) and plugging it into Eq. (5) and can be expressed in terms of the SCGF
[12] as

$h(s)=\Psi(s)+(1-s)\Psi^{\prime}(s)\ .$ (12)

This expression can be analytically proved to have a global maximum for $s=1$,
i.e.,

$h(1)=\Psi(1)=\ln\zeta_{1}\ ,$ (13)

shown in Fig. 1(b), where $\zeta_{1}$ is the dominant eigenvalue of the
adjacency matrix $A$. Replacing $s=1$ in the driven process, Eq. (10) gives
the MERW

$\left(\pi_{1}\right)_{ij}=\frac{a_{ij}r_{1}(j)}{r_{1}(i)\zeta_{1}}\ ,$ (14)

allowing us to interpret a MERW on a network as a biased random walk creating
a particular rare event fluctuation—given by replacing $s=1$ in Eq. (11)—of an
URW [12].

This result shows, on the one hand, that we can sample a particular rare event
of an URW by simulating a MERW and, on the other hand, that a MERW can in
principle be obtained from an URW by opportunely conditioning on a certain
rare event of the observable $C_{n}$ in Eq. (21). The latter observation is
key to introduce our adaptive random walk (ARW), which through successive
local adaptations of an URW converges to a MERW on the whole network. In fact,
as we will show in the following, such ARW reaches maximum entropy while
exploring the network, i.e., much before the entire graph has been visited.

To construct the ARW we develop here an algorithm based on a rare-event
sampling scheme firstly proposed in the context of reinforcement learning
methods for risk assessment in Markov chain models [46, 47], later extended to
continuous-time processes [48, 49], and recently adapted to the study of
random walks on networks [36]. In our algorithm, the random walk updates its
transition probability matrix at each time step, in order to typically visit a
specific rare event of the URW, obtained fixing $s=1$ in Eq. (11). We refer
the reader to the Supplemental Material (SM) for a general formulation of the
sampling algorithm valid for all additive observables and $s\in\mathbb{R}$.

Formally, the ARW is a discrete-time process $Y=(Y_{1},Y_{2},\cdots,Y_{n})$,
where $Y_{n}\in V$ is the position of the ARW on the graph at time $n$. The
core of the ARW is based on an adaptive power method to solve the following
dominant eigenvalue equation

$\tilde{\Pi}_{1}r_{1}=\zeta_{1}r_{1}\ ,$ (15)

which is known to be cardinal to construct the MERW in Eq. (14). In fact, in a
large deviation context, one seeks to solve Eq. (15) to determine the SCGF
$\Psi(1)$ via Eq. (8). As previously discussed, this SCGF is linked to the
probability of the rare event $c=\Psi^{\prime}(1)$, typically generated by a
MERW.

More in detail, our adaptive power method simulates single Markov chain
transitions with importance sampling [50]. In particular, supposing that the
ARW is located on the node $i$ at time $n$, i.e., $Y_{n}=i$, the next step is
proposed according to an estimate of the MERW in Eq. (14) that reads

$\left(\pi_{1}^{(n)}\right)_{ij}=\frac{1}{Z}\frac{a_{ij}r_{1}^{(n)}(j)}{r_{1}^{(n)}(i)r_{1}^{(n)}(i_{0})}=\frac{a_{ij}r_{1}^{(n)}(j)}{\sum_{j^{\prime}\in
V}a_{ij^{\prime}}r_{1}^{(n)}(j^{\prime})}\ ,$ (16)

where $Z$ is the normalisation factor, $i_{0}$ is an a-priori fixed node, and
$r_{1}^{(n)}$ is the $n$-th time estimate of the eigenvector centrality [36].
This last is given by the stochastic-approximation [51] formula

$\begin{split}&r_{1}^{(n+1)}(i)=r_{1}^{(n)}(i)+\\\ &\hskip
28.45274pta(n)\mathbf{1}_{Y_{n}=i}\left(\frac{\sum_{j\in
V}a_{ij}r_{1}^{(n)}(j)}{r_{1}^{(n)}(i_{0})}-r_{1}^{(n)}(i)\right)\end{split}$
(17)

based on the learning rate $a(n)$ and on an asynchronous update via the
indicator function $\mathbf{1}_{Y_{n}=i}$. In other words, the $i$-th
component of the eigenvector centrality is updated only when the process $Y$,
at the $n$-th time step, has visited node $i$. Also note that, in the
following, all ARW simulations are obtained by fixing a learning rate
$a(n)=1/((n+1)^{0.1})$ in Eq. (17). For technical details on the adaptive
scheme we address the reader to the SM.

Figure 2: (a-c) The entropy production rate $h$, calculated as in Eq. (19) of
single trajectories of ARW, MERW and URW (solid lines) is compared to the
corresponding optimal entropy production rate $\bar{h}$ on the discovered
graph (dashed lines). (d-f) Median (solid lines) and first and third quartiles
(shaded area) of the normalized differences $(\bar{h}-h)/\bar{h}$ over an
ensemble of $1000$ trajectories of ARW, MERW, and URW. Results are shown for
random walks running on the giant connected component of an Erdös–Rényi random
graph with $1000$ nodes and average degree $3$ in (a,d), a Barabasi–Albert
with $1000$ nodes and $m=4$ in (b,e), and an air transportation network with
$3618$ nodes and $14142$ links in (c,f).

The ARW is randomly initialized on a node of the network with a normalized
random right eigenvector $r_{1}^{(0)}$ and evolves according to the fully
local rules in Eq. (16) and Eq. (17). At each time step it tends to optimize
the spreading—aiming to the maximum entropy production—on the portion of the
network visited. In the long-time limit, it will eventually converge to the
MERW of Eq. (14) as

$r_{1}^{(n)}\rightarrow r_{1},\hskip
17.07182ptr_{1}^{(n)}(i_{0})\rightarrow\zeta_{1},\hskip
8.5359pt\text{and}\hskip 8.5359ptZ\rightarrow 1\ .$ (18)

We insist on the fact that differently from the MERW, the ARW does not need to
know the full topology of the network since the beginning, as it learns it on
the run. Thanks to this, its entropy production rate stays always close to the
maximum rate on the visited portion of the graph. This is drastically
different from a MERW which does not maximize the entropy while it explores
the graph, but reaches optimal spreading only when the whole network has been
visited.

We show this in Fig. 2(a)–(c) where we compare the spreading properties of
single trajectories of ARW, MERW, and URW on two network models, namely
Erdös–Rényi, and Barabasi–Albert, and on a real-world man-made network. The
last network describes an air transportation system, with each node being a
city and two nodes are connected if at least one airplane flew between the two
cities in the time window [1st Nov 2000 - 31st Oct 2001] [52]. Each of the
three processes is initialized on a randomly selected node of the network, and
evolves according to its transition probability matrix. As the walker moves to
(previously unvisited) nodes, we calculate its entropy production rate $h$
(solid line) and compare it with the optimal $\bar{h}$ (dashed line) given by
the logarithm of the dominant eigenvalue of the adjacency matrix associated
with the portion of graph made by all (and only) previously visited links. The
entropy production rate $h$ of each process is calculated via a modified
version of Eq. (5), which is

$h(M)=-\sum_{(i,j)\in E(n)}\rho_{i}(M)\pi_{ij}(M)\ln\pi_{ij}(M)\ .$ (19)

This takes into account the number of visited links $M=|E(n)|$ up to time $n$,
where $E(n)$ is the set of visited links,
$\pi_{ij}(M)=\pi_{ij}/(\sum_{(i,j^{\prime})\in E(n)}\pi_{ij^{\prime}})$ if the
link $(i,j)$ is in $E(n)$, and $0$ otherwise, and the stationary distribution
$\rho(M)$ is calculated as the left eigenvector of
$\Pi(M)=\left\\{\pi_{ij}(M)\right\\}$. As the solid line is closer to the
corresponding dashed one for the ARW with respect to the MERW and the URW, we
notice that the ARW reaches earlier—while exploring the network—best spreading
performances. When the whole graph has been visited, the ARW performance is
comparable to that of the MERW with the entropy production rate in Eq. (13),
whereas the URW converges to Eq. (5).

Moreover, in Fig. 2(d)–(f) we plot the median (solid line) of the relative
difference between the entropy production rate $h$ and the optimal $\bar{h}$,
together with first and third quartiles (shaded area). The median and the
quartiles are calculated over an ensemble of $1000$ trajectories and are used
in place of the mean and standard deviation because of the unknown
distribution of $h$ around $\bar{h}$ (certainly not symmetric as
$h\leq\bar{h}$). This gives further evidence of the fact that the ARW performs
better than the MERW (and the URW) at maximizing the spreading while
discovering the structure of the graph.

Fig. 1 of the SM shows the entropy production rates of ARW, MERW, and URW on
the graph induced by the visited nodes, i.e, also considering links that have
not yet been visited. In this case, the ARW is still outperforming MERW and
URW on the Erdös–Rényi graph, but is only marginally better on the scale-free
and the real-world network. Therefore, this other study confirms the good
choice of comparison for the entropy production rates in Fig. 2: this needs to
be calculated on the portion of network only composed by the visited links.

The optimal dispersion of the fully-local ARW on the visited links of the
network comes at a price: the ARW takes longer than URW and MERW to cover the
whole graph. This is a consequence of two properties of the dynamics of the
ARW: (i) in the beginning, the ARW spends a certain time localised in the
region where it has been initialised—we refer to this time window as a ‘warm-
up’ time—(ii) once the warm-up is over, the ARW visits the network with a
typical exploration time. During the warm-up time, our process finely tunes
the eigenvector centrality and the transition probability matrix to set off an
efficient exploration of the network. In Fig. 2 of the SM we show how many
time steps are needed, on average, to discover new links, analysing an
ensemble of 1000 trajectories of URW, MERW, and ARW for each graph used in
Fig. 2. We notice that, in all the cases, the ARW remains localised in the
first few visited links of the network during the warm-up phase. This is also
the main reason why the network coverage time—the number of time steps to
visit all links—is, on average, $10^{2}-10^{4}$ time steps longer for an ARW
than for a MERW or an URW. We compare the ARW coverage time with URW and MERW
in Fig. 3

Figure 3: Coverage time $T$ as a function of $L$, i.e., the number of links in
the giant connected component of Erdös-Rényi graphs with average degree $3$
and increasing size for ARW, MERW, and URW.

where we plot $T$ as a function of the number of links $L$ in the giant
connected component of Erdös–Rényi graphs with average degree $3$ and
increasing number of nodes. We find that that the coverage time $T$ is related
to the total number of links $L$ by the equation

$T\propto L^{\alpha}\ ,$ (20)

where $\alpha$ is the scaling exponent. Remarkably, after the initial warm
up—which is evident by the overall upward shift of the ARW power-law fit—the
scaling exponents of ARW and MERW are similar ($\alpha_{A}\approx 2.9$,
$\alpha_{M}\approx 2.7$), but larger than the one of the URW which has no
global constraint to satisfy ($\alpha_{U}\approx 1.3$).

As a last remark, we point out that in our studies the initial conditions for
the ARW—in particular, the initial right eigenvector $r_{1}^{(0)}$—play a key
role in determining the initial warm-up time and the overall accuracy of the
ARW in optimizing the spreading while exploring the network. In Fig. 3 of the
SM we collect several plots that compare the spreading and coverage time
performances of another adaptive random walk, dubbed $\overline{\text{ARW}}$,
with MERW, URW, and the ARW discussed so far. The $\overline{\text{ARW}}$
evolves in time according to the same algorithm of the ARW, but differently
from the latter it is initialized with a non-normalized right eigenvector.
Noticeably, the warm-up time—therefore, also the coverage time—is greatly
reduced with respect to the ARW at the expense of a spreading that is
comparable with a MERW and hence not optimal while exploring the network. We
hope that our work will stimulate further investigations to explore such
trade-off between optimal spreading accuracy and exploration times.

In this manuscript, we have proposed a random walk process, named ARW, that
has optimal spreading properties while exploring a network. As such, it
outperforms the well-known MERW which, requiring the knowledge of the whole
network beforehand, is non-local in the choice of the next visited node. The
ARW is instead based on a random walk that by means of only local rules adapts
itself to be a maximum entropy random walk while visiting the network, without
requiring any prior knowledge of it. The adaptive step is the fundamental
ingredient of the process. Via a large-deviation tilting on the fluctuating
trajectory entropy observable, the adaptive step allows the ARW to typically
observe a maximum entropy production rate, which qualifies as a rare event of
the URW. We have also shown two key features of the ARW dynamics: the initial
warm-up time—independent on the size of the graph—that allows the ARW to set
off an optimal spreading on the network, and the scaling exponent of the
network coverage time, which is similar to that of the MERW. Besides the
theoretical novelty driven by a large deviation study of random-walk rare
events, we believe that the ARW might find useful applications to network
information spreading [53, 54] in all such cases where no (or only partial)
prior knowledge on the network is available, or when the network is changing
in time [55, 56, 57]. It could also be used to study dispersion properties of
other dynamical processes on real networks, e.g., aiming at optimal
exploration in congested networks [58, 59]. Furthermore, the algorithm at the
core of the ARW could also be used to sample other rare-event fluctuations
associated with any additive observables of random walks [36].

##  I Codes

All the codes used in the manuscript are available at
https://github.com/gabriele-di-bona/ARW.

##  II Supplemental Material

##  III Adaptive random walk algorithm

In this Section we give details on the adaptive random walk process
$Y=(Y_{\ell})_{\ell=1}^{n}$ which is based on an efficient algorithm to sample
large-deviation rare events. Differently from the main text, we will not
restrict to the particular case of a MERW, but we will describe the numerical
scheme in its full generality. However, all the results in the main text (MT)
can be obtained fixing $f(X_{\ell})=k_{X_{\ell}}$ in Eq. (21) and replacing
$s=1$ in (22) and following.

We consider a Markov chain $X=(X_{\ell})_{\ell=1}^{n}$ evolving in a finite
discrete state space $\Gamma$ of $N$ states (this is the graph $G=(V,E)$ in
the MT) according to the (irreducible and aperiodic) transition matrix $\Pi$
with elements $\pi_{ij}$ that characterize the probability of going from a
state $X_{\ell}=i$ at time $\ell$ to a state $X_{\ell+1}=j$ at time $\ell+1$.
Associated with the Markov chain $X$, we consider a purely time-additive
observable of the general form

$C_{n}=\frac{1}{n}\sum_{\ell=1}^{n}f(X_{\ell})\ ,$ (21)

where $f$ is any function of the state. As mentioned in the MT, one can
estimate the likelihood of fluctuations of $C_{n}$ by calculating the scaled
cumulant generating function (SCGF) in Eq. (7)MT (and Eq. (8)MT) and by
Legendre–Fenchel transforming it to obtain the large deviation rate function,
i.e., the leading exponential behaviour of the probability density in Eq.
(6)MT. However, this procedure is known to be difficult since it involves the
calculation of the dominant (Perron–Frobenius) eigenvalue $\zeta_{s}$ of the
non-negative tilted matrix
$\tilde{\Pi}_{s}=\left\\{\left(\tilde{\pi}_{s}\right)_{ij}\right\\}$ having
elements

$\left(\tilde{\pi}_{s}\right)_{ij}=\pi_{ij}e^{sf(i)}\ ,$ (22)

see [39, 41, 43, 36] for further details. This problem can be mathematically
formalized as finding the solution of the following spectral equation

$\tilde{\Pi}_{s}r_{s}=\zeta_{s}r_{s}\ ,$ (23)

where along with $\zeta_{s}$ we have the right eigenvector $r_{s}$. Both these
are fundamental to define the driven process, which is characterized by a
transition matrix whose components are

$\left(\pi_{s}\right)_{ij}=\frac{\left(\tilde{\pi}_{s}\right)_{ij}r_{s}(j)}{r_{s}(i)\zeta_{s}}=\frac{\pi_{ij}e^{sf(i)}r_{s}(j)}{r_{s}(i)\zeta_{s}}\
,$ (24)

as in the first equality of Eq. (10)MT.

In the following, we show a numerical method that approximates the dominant
right eigenvector $r_{s}$ and, consequently, the dominant eigenvalue
$\zeta_{s}$, and so the SCGF $\Psi$ too (Eq. (8)MT), by making use of an
adaptive stochastic power-method scheme. The core of the algorithm is rooted
on works of reinforcement learning for risk-sensitive control of Markov chains
[60, 46, 47, 61] and was also recently adapted to estimate large deviation
functions of continuous-time diffusion processes [48, 49]. Here, we repropose
a discussion of the method that appeared in [36].

A bare power method scheme recursively multiplies an arbitrary vector
$v_{0}\coloneqq r_{s}^{(0)}$ (the superscript $(0)$ refers to the zero-th step
in the numerical scheme) by the transition matrix $\tilde{\Pi}_{s}$. In this
case, we can write

$\zeta_{s}^{(n+1)}r_{s}^{(n+1)}=\tilde{\Pi}_{s}^{n+1}r_{s}^{(0)}=\tilde{\Pi}_{s}r_{s}^{(n)}$
(25)

where $r_{s}^{(n+1)}$ is the approximation at the $(n+1)$-th step of the
dominant right eigenvector normalized by the $(n+1)$-th estimate of the
dominant eigenvalue $\zeta_{s}^{(n+1)}$. This can be calculated as

$\zeta_{s}^{(n)}=r_{s}^{(n)}(i_{0})\ ,$ (26)

where, a priori, the component $i_{0}$ is arbitrarily chosen, but here we fix
it to be the maximum component of $r_{s}^{(n)}$ in magnitude such that the
right eigenvector is normalized according to its infinity norm111The infinity
norm of a vector is the maximum of the vector’s absolute-value components.. As
$n\rightarrow\infty$ we then have

$r_{s}^{(n)}\rightarrow r_{s}\hskip 28.45274pt\text{and}\hskip
28.45274ptr_{s}^{(n)}(i_{0})\rightarrow\zeta_{s}\ .$ (27)

Numerically, it is known that the power method scheme in Eq. (25) is
particularly inefficient for large (and very connected) state spaces, because
at each step one matrix multiplication needs to be calculated. To overcome
this problem, as we will see in the following, one may consider to implement a
numerical scheme that simulates single Markov chain transitions rather than
keeping track of the full structure of the state space in the matrix
multiplication. These transitions could be naively proposed according to the
tilted matrix in Eq. (22). Notwithstanding this, due to the presence of the
exponential factor, such an approach could lead to divergences. For this
reason, one considers importance sampling [50] and simulates the $(n+1)$-th
transition step according to the $n$-th estimate of the driven process in Eq.
(24), that is

$\left.(\pi_{s}^{(n)})_{ij}=\frac{\pi_{ij}e^{sf(i)}r^{(n)}_{s}(j)}{r^{(n)}_{s}(i)r^{(n)}_{s}(i_{0})}\middle/\left(\sum_{j^{\prime}=1}^{N}\frac{\pi_{ij^{\prime}}e^{sf(i)}r^{(n)}_{s}(j^{\prime})}{r^{(n)}_{s}(i)r^{(n)}_{s}(i_{0})}\right)\right.=\frac{\pi_{ij}r^{(n)}_{s}(j)}{\sum_{j^{\prime}=1}^{N}\pi_{ij^{\prime}}r^{(n)}_{s}(j^{\prime})}\
.$ (28)

With a simulated transition from $i$ to $j$ the estimate of the right
eigenvector in Eq. (25) therefore reads

$r^{(n+1)}_{s}(i)=\frac{e^{sf(i)}}{r^{(n)}_{s}(i_{0})}r^{(n)}_{s}(j)\frac{\pi_{ij}}{(\pi_{s})^{(n)}_{ij}}=\frac{e^{sf(i)}\sum_{j^{\prime}=1}^{N}\pi_{ij^{\prime}}r_{s}^{(n)}(j^{\prime})}{r_{s}^{(n)}(i_{0})}\
,$ (29)

where we normalize the right eigenvector with $r^{(n)}_{s}(i_{0})$ ($i_{0}$
is, once again, the maximum component of $r_{s}^{(n)}$ in magnitude) and, as
usual in importance sampling, we unbias the transition multiplying by the
likelihood ratio222Discrete Radon–Nikodym derivative.
$\pi_{ij}/(\pi_{s}^{(n)})_{ij}$. Notice that it is only in the limit
$n\rightarrow\infty$ that the estimate in Eq. (28) tends to the driven process
of (24). Indeed, the normalisation appearing in (28) has a pivotal role in
correctly estimating the driven process transition probabilities, and it is
only in the limit $n\rightarrow\infty$ that such normalisation converges to
$1$.

Eventually, convergence of Eq. (29), as in Eq. (27), is realized making use of
a stochastic approximation scheme. Stochastic approximation methods [51, 62]
are discrete-time (stochastic) recursive schemes similar, in form, to gradient
descent methods. They make use of an annealing schedule, also called learning
rate, to quantify, at each time step, ‘how much’ of the new information coming
from an algorithm is used to update the value of an observable of interest in
the next step. In our case, this observable of interest is the dominant right-
eigenvector, and the stochastic approximation scheme applied to the algorithm
in Eq. (29) reads

$r^{(n+1)}_{s}(i)=r^{(n)}_{s}(i)+a(n)\mathbf{1}_{Y_{n}=i}\left(\frac{e^{sf(i)}\sum_{j=1}^{N}\pi_{ij}r_{s}^{(n)}(j)}{r_{s}^{(n)}(i_{0})}-r^{(n)}_{s}(i)\right)\
,$ (30)

where $\mathbf{1}_{Y_{n}=i}$ makes the update asynchronous, viz. it selects
only the state $i$ at time step $n$, and $a$ is the annealing schedule, which
should satisfy

$\lim_{n\rightarrow\infty}\sum_{\ell=1}^{n}a(\ell)=\infty\hskip
56.9055pt\lim_{n\rightarrow\infty}\sum_{\ell=1}^{n}a^{2}(\ell)<\infty.$ (31)

A proof of the convergence of this numerical method is given in [51], where it
is shown that the numerical scheme tracks the solution of a particular
nonautonomous differential equation.

The remarkable advantage of this adaptive stochastic numerical scheme, based
on the combination of spectral methods with important sampling, is that it
allows us to estimate the SCGF $\Psi$ by simulating a single trajectory of a
Markov chain that learns on the run the particular rare event of interest $c$
of Eq. (21), linked to $s$ by $\Psi^{\prime}(s)=c$. For instance, we simulate
the ARW described in the MT with Eq. (8) and Eq. (10) by fixing $s=1$, i.e.,
Eq. (16)MT and Eq. (17)MT, to adaptively obtain a random walk with maximum
entropy production on every visited portion of the network. Moreover, in a
large deviation context, one can estimate the full SCGF by discretizing with a
fine mesh the tilting parameter space, i.e., $[\dots,s-\Delta s,s,s+\Delta
s,\dots]$, and running the algorithm for every point in the grid. Eventually,
the estimated SCGF is calculated by interpolating all the values obtained.

Generally speaking, a direct implementation of the numerical scheme presented
above, although extremely efficient from a computational point of view, could
still lead to noisy results in the calculation of the dominant eigenvalue and
right eigenvector in absence of good initial guesses for $r_{s}$ and
$\zeta_{s}$. In the simulations used for Fig. (2)MT and Fig. (3)MT we have
initialized the vector $r_{s}$ with random components sampled from the $[0,1]$
uniform distribution and then normalized the vector according to the $L^{1}$
norm. We will detail this point in Section VI. We also address the reader to
Chapter 4 of [36] for a more general discussion.

###  III.0.1 Pseudocode

For full transparency of the work done, we propose in the following a
pseudocode for the ARW, obtained by implementing the algorithm presented in
the previous Section III for the particular case of a transition matrix $\Pi$
associated with an URW and tilting parameter $s=1$. We address the reader to
the GitHub folder (see CODES Section above) where we have collected all the
codes used to draft the manuscript.

1

Input : $\Pi$, $N$ (actual size of the graph), $k$ (degree list), $\alpha$
(learning rate exponent), Max (maximum number of iterations)

Output : right dominant eigenvector: $r_{1}$, dominant eigenvalue: $\zeta_{1}$

2

// Initialisation

3

4$r_{1}=$ $L^{1}$-norm rnd list ;

5 $\rho=0$ list initialisation of the (non-normalized) empirical occupation
measure ;

6 $\text{start}=$ rnd($N$) ;

7 $\rho(\text{start})\mathrel{{+}{=}}1$ ;

8 $i_{0}$ = rnd($N$) ;

9

10for  _$n\leftarrow 0$ to Max_ do

11

// setting the learning rate

12 $a(n)=\frac{1}{(1+n)^{\alpha}}$ ;

13 if  _$n!=0$_ then

14 start $=$ next ;

// Empirical occupation measure

15 $\rho(\text{start})\mathrel{{+}{=}}1$ ;

16 for  _$j\leftarrow 1$ to $N$_ do

// step 1: driven process update in Eq. (28).

17
$\pi^{(n)}_{\text{start}j}=\frac{\pi^{(n-1)}_{\text{start}j}r_{1}(j)}{\sum_{j^{\prime}=1}^{N}\pi^{(n-1)}_{\text{start}j^{\prime}}r_{1}(j^{\prime})}$
;

18

19 end for

20

21 end if

22

// step 2: pick a random index according to the probability vector
$\pi^{(n)}_{\text{start}}$ (biased random walk).

23 next $=$ rnd($N$,$\pi^{(n)}_{\text{start}}$) ;

24

// step 3: asynchronous update of the right eigenvector in Eq. (30).

25
$r_{1}(\text{start})=r_{1}(\text{start})+a(n)\left(\frac{k(\text{start})\sum_{j\in\partial\text{start}}r_{1}(j)}{k_{\text{start}}r_{1}(i_{0})}-r_{1}(\text{start})\right)$
;

26

27 $i_{0}=\max(\rho)$ ;

28

29 end for

return $r_{1}$, $\zeta_{1}\coloneqq r_{1}(i_{0})$

Algorithm 1 Adaptive random walk (ARW)

##  IV Measuring the entropy production rate

In Fig. 2MT we compare the spreading performances of ARW, MERW, and URW by
calculating the entropy production rates of each process while they explore
the network. Clearly, there is some freedom in the choice of what to consider
as explored up to a certain moment in the evolution of each process. The most
conservative choice, which is implemented in the calculations appearing in
Fig. 2MT, is to consider the graph visited composed only by the links (and
connecting nodes) that have actually been crossed by the process in
consideration. This is the case we face in the reality when, for instance, a
person takes choices leading them to unknown situations that, in turn, are
connected with other choices that, however, they do not know until the moment
they face that particular situation. In such a scenario, the stochastic matrix
entering in Eq. (19)MT at a certain time $n$ for each process will have to
have non-zero entries only for the edges effectively crossed. This means that
a redesign of the transition matrix—cancelling non-zero entries for unvisited
links—for all the process while they run on the network needs to be done.

Another possibility is to consider the graph discovered composed by all the
links that connect visited nodes. In the real world, this case happens every
time one does not need to cross a link between two nodes to visit them if
these have already been visited via other paths. Furthermore, a more
theoretical reason pushing for this attempt comes from the implementation of
the ARW as in the algorithm 1: the right eigenvector is built on the nodes,
not on the links, and is updated every time a node is visited, independently
on the link used. Similarly, the transition matrix of the estimated driven
process is updated for all the nodes connected to the one the random walk has
just visited. Similarly to the previous case, the transition matrices of ARW,
MERW and URW need to be redesigned—cancelling non-zero entries for links that
connect at least one unvisited node—every time a new node in the network is
discovered. Therefore, in this scenario, the entropy production rate is
calculated as

$h(M)=-\sum_{i,j\in V(n)}\rho_{i}(M)\pi_{ij}(M)\ln\pi_{ij}(M)\ ,$ (32)

where $V(n)$ is the set of visited nodes up to time $n$, and with a little
abuse of notation $\pi_{ij}(M)=\pi_{ij}/(\sum_{i,j^{\prime}\in
V(n)}\pi_{ij^{\prime}})$ if $(i,j)\in V(n)$, $0$ otherwise, and $\rho(M)$ is
calculated as the dominant left eigenvector of $\Pi(M)$ as in the case treated
in the MT.

In Fig. 4 we re-propose Fig. 2MT for the entropy production rate calculated
with Eq. (32).

Figure 4: (a-c) The entropy production rate $h$, calculated as in Eq. (32), of
single trajectories of ARW, MERW and URW (solid lines) is compared to the
corresponding optimal entropy production rate $\bar{h}$ on the discovered
graph (dashed lines). (d-f) Median (solid lines) and first and third quartiles
(shaded area) of the normalized differences $(\bar{h}-h)/\bar{h}$ over an
ensemble of $1000$ trajectories of ARW, MERW, and URW. Results are shown for
random walks running on the giant connected component of an Erdös–Rényi random
graph with $1000$ nodes and average degree $3$ in (a,d), a Barabasi–Albert
network with $1000$ nodes and $m=4$ in (b,e), and an air transportation
network with $3618$ nodes and $14142$ links in (c,f).

Noticeably, the ARW, obtained with a learning rate $a(n)=1/((n+1)^{0.1})$ in
Eq. 16MT or Eq. (30) for $s=1$, is only marginally better than MERW and URW on
the Barabasi–Albert and on the air transportation network (the same analyzed
in the MT), although still much better over the Erdös–Rényi random graph. This
is evidence of the fact that the ARW is effectively optimized only over the
links it actually crosses. For this reason, we argue that the entropy
production rate needs to be calculated over the graph of visited links as done
in Fig. 1MT.

##  V Link-exploration time

As mentioned in the MT, the ARW is characterized by a long warm-up time that
prolongs the time to explore other links in the network. During the warm up,
the ARW is localized in very few links (1 to 3) and by continuously crossing
them it finely tunes the eigenvector centrality and the transition matrix
estimate of the driven process in the visited portion of the graph, so as to
set off an optimal exploration of the network. Indeed, once the warm up is
over, the coverage time scaling exponent of the ARW is very similar to that of
the MERW (see Fig. 3MT). In Fig. 5 we compare the link exploration time of
MERW, URW, and ARW over the same networks considered in Fig. 2MT and in Fig.
4. For each simulation, we save every time step in which the random walk
explores a new link for the first time, and we plot them for a random instance
of each process in Fig. 5(a-c), while in Fig. 5(d-f) we plot the medians,
first and third quartiles, taken over $1000$ trajectories. Differently from
Fig. (3)MT, here we can analyze how the discovery rate of the process
decreases or increases while the network is discovered highlighting the
temporal extension of the warm-up time as well as the number of links
involved.

Figure 5: (a-c) Link exploration time $T$ (solid line) for single trajectories
of ARW, MERW, and URW. (d-f) Median (solid line), with first and third
quartiles (shaded area), of the link exploration time for $1000$ trajectories.
Results are shown for random walks running on the giant connected component of
an Erdös–Rényi random graph of $1000$ nodes and average degree $3$ in (a) and
(d), a Barabasi–Albert network of $1000$ nodes and $m=4$ in (b) and (e), and
an air transportation network with $3618$ nodes and $14142$ links in (c) and
(f).

As mentioned, this warm-up time is responsible for the longer—of about
$10^{2}$ to $10^{4}$ steps—link coverage time of the ARW with respect to MERW
and URW.

##  VI Ideas on how to optimize the Adaptive Random Walk

Such a long warm-up time seems to be necessary for the ARW to set off optimal
initial conditions (in terms of right eigenvector and driven process
transition matrix) for the exploration of the entire network. However, it also
seems that the warm-up time can be shortened at the expense of the optimal
diffusivity of the process. Indeed, the initial condition for the right
eigenvector and the entropy production rate play as two conjugated quantities:
the choice of the initial right eigenvector $r^{(0)}_{1}$ influences both the
spreading properties and the overall coverage time of the process. In the MT
we considered $r^{(0)}_{1}$ to be an $L^{1}$-normalized random vector. In Fig.
6 instead, we study the new $\overline{\text{ARW}}$, which is essentially
equivalent to the ARW so far considered, apart from $r^{(0)}_{1}$ that, in
this case, is a non-normalized random vector with components sampled from the
uniform distribution over $[0,1]$.

Figure 6: (a) Median (solid line), and first and third quartiles (shaded
area), of the link exploration time for $1000$ trajectories of URW, MERW, ARW,
and $\overline{\text{ARW}}$ on the giant connected component of an Erdös–Rényi
random graph of $100$ nodes and average degree $3$. (b) Coverage time $T$ as a
function of $L$, the number of links in the giant connected component of
different Erdös-Rényi graphs with average degree $3$ for URW, MERW, ARW, and
$\overline{\text{ARW}}$. (c) Comparison between entropy production rate $h$
(solid line) and optimal entropy production rate $\bar{h}$ (dashed line) on
the discovered graph, among single trajectories of URW, MERW, ARW, and
$\overline{\text{ARW}}$. (d) Median (solid line), and first and third
quartiles (shaded area), of the entropy production rate relative difference
for $1000$ trajectories of URW, MERW, ARW, and $\overline{\text{ARW}}$.

Interestingly, $\overline{\text{ARW}}$ is much quicker than ARW in exploring
the network—it has a smaller scaling exponent
$\alpha_{\overline{\text{A}}}$—see Fig. 6 (a) and (b) and it has about the
same overall link coverage time of MERW. Furthermore, similarly to the MERW,
it drastically decelerates in the last phase (Fig. 6(a)). Notwithstanding
this, its spreading performance are, on average, no better than MERW (Fig.
6(c)(d)) and certainly not as good as ARW.

This is only a preliminary study that shows an interesting interplay between
initial conditions and overall spreading performances of the adaptive random
walk. Further studies in this direction may help to find the best process that
shows short link exploring/coverage time and maximal entropy production rate
while exploring the network.

## References

  * [1] S. Albeverio, V. Jentsch, and H. Kantz, eds., Extreme Events in Nature and Society.  The Frontiers Collection, Berlin, Heidelberg: Springer Berlin Heidelberg, 2006.
  * [2] V. Kishore, M. S. Santhanam, and R. E. Amritkar, “Extreme events on complex networks,” Physical Review Letters, vol. 106, no. 13, p. 188701, 2011\.
  * [3] V. Kishore, M. S. Santhanam, and R. E. Amritkar, “Extreme events and event size fluctuations in biased random walks on networks,” Physical Review E - Statistical, Nonlinear, and Soft Matter Physics, vol. 85, no. 5, p. 056120, 2012.
  * [4] Y. Z. Chen, Z. G. Huang, and Y. C. Lai, “Controlling extreme events on complex networks,” Scientific Reports, vol. 4, no. 1, pp. 1–10, 2014\.
  * [5] Y. Z. Chen, Z. G. Huang, H. F. Zhang, D. Eisenberg, T. P. Seager, and Y. C. Lai, “Extreme events in multilayer, interdependent complex networks and control,” Scientific Reports, vol. 5, no. 1, pp. 1–13, 2015.
  * [6] J. Hindes and I. B. Schwartz, “Epidemic extinction and control in heterogeneous networks,” Physical Review Letters, vol. 117, no. 2, p. 028302, 2016.
  * [7] C. De Bacco, A. Guggiola, R. Kühn, and P. Paga, “Rare events statistics of random walks on networks: localisation and other dynamical phase transitions,” Journal of Physics A: Mathematical and Theoretical, vol. 49, no. 18, p. 184003, 2016.
  * [8] J. Hindes and I. B. Schwartz, “Large order fluctuations, switching, and control in complex networks,” Scientific Reports, vol. 7, no. 1, pp. 1–9, 2017.
  * [9] G. Bianconi, “Rare events and discontinuous percolation transitions,” Physical Review E, vol. 97, no. 2, p. 022314, 2018.
  * [10] F. Coghi, F. Radicchi, and G. Bianconi, “Controlling the uncertain response of real multiplex networks to random damage,” Physical Review E, vol. 98, no. 6, p. 062317, 2018.
  * [11] W. Staffeldt and A. K. Hartmann, “Rare-event properties of the Nagel-Schreckenberg model,” Physical Review E, vol. 100, no. 6, p. 062301, 2019.
  * [12] F. Coghi, J. Morand, and H. Touchette, “Large deviations of random walks on random graphs,” Physical Review E, vol. 99, no. 2, p. 022137, 2019.
  * [13] R. Gutierrez and C. Perez-Espigares, “Generalized optimal paths and weight distributions revealed through the large deviations of random walks on networks,” Physical Review E, vol. 103, no. 2, p. 022319, 2021.
  * [14] G. Carugno, P. Vivo, and F. Coghi, “Delocalisation-localisation dynamical phase transition of random walks on graphs,” In preparation, 2021.
  * [15] A. Kumar, S. Kulkarni, and M. S. Santhanam, “Extreme events in stochastic transport on networks,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 30, no. 4, p. 043111, 2020.
  * [16] K. Gupta and M. S. Santhanam, “Extreme events in Nagel–Schreckenberg model of traffic flow on complex networks,” The European Physical Journal Special Topics, vol. 230, no. 16-17, pp. 3201–3209, 2021.
  * [17] G. Gandhi and M. S. Santhanam, “Biased random walkers and extreme events on the edges of complex networks,” arXiv:2110.14949v1.
  * [18] J. D. Noh and H. Rieger, “Random walks on complex networks,” Physical Review Letters, vol. 92, no. 11, p. 118701, 2004.
  * [19] J. Gómez-Gardeñes and V. Latora, “Entropy rate of diffusion processes on complex networks,” Physical Review E, vol. 78, no. 6, p. 065102, 2008.
  * [20] Z. Burda, J. Duda, J. M. Luck, and B. Waclaw, “Localization of the maximal entropy random walk,” Physical Review Letters, vol. 102, no. 16, p. 160602, 2009.
  * [21] L. Demetrius and T. Manke, “Robustness and network evolution - an entropic principle,” Physica A: Statistical Mechanics and its Applications, vol. 346, no. 3-4, pp. 682–696, 2005.
  * [22] J.-C. Delvenne and A.-S. Libert, “Centrality measures and thermodynamic formalism for complex networks,” Physical Review E, vol. 83, no. 4, p. 046117, 2011.
  * [23] Y. Lin and Z. Zhang, “Mean first-passage time for maximal-entropy random walks in complex networks,” Scientific Reports, vol. 4, no. 1, pp. 1–7, 2014.
  * [24] F. Battiston, V. Nicosia, and V. Latora, “Efficient exploration of multiplex networks,” New Journal of Physics, vol. 18, no. 4, p. 043035, 2016.
  * [25] E. Estrada, “Informational cost and networks navigability,” Applied Mathematics and Computation, vol. 397, p. 125914, 2021.
  * [26] Y. Wang, X. Cao, T. Weng, H. Yang, and C. Gu, “Lowest-degree preference random walks on complex networks,” Physica A, vol. 577, p. 126075, 2021\.
  * [27] J. K. Ochab and Z. Burda, “Maximal entropy random walk in community finding,” The European Physical Journal Special Topics, vol. 216, no. 1, pp. 73–81, 2012.
  * [28] J. K. Ochab, “Maximal-entropy random walk unifies centrality measures,” Physical Review E, vol. 86, no. 6, p. 066109, 2012.
  * [29] Z. Xu, C. Pu, and J. Yang, “Link prediction based on path entropy,” Physica A: Statistical Mechanics and its Applications, vol. 456, pp. 294–301, 2016.
  * [30] I. Adam, D. Fanelli, T. Carletti, and G. Innocenti, “Reactive explorers to unravel network topology,” The European Physical Journal B, vol. 92, no. 5, pp. 1–8, 2019.
  * [31] Y. Chen, T. Georgiou, M. Pavon, and A. Tannenbaum, “Robust transport over networks,” IEEE transactions on automatic control, vol. 62, no. 9, p. 4675, 2017.
  * [32] Y. W. Niu, H. Liu, G. H. Wang, and G. Y. Yan, “Maximal entropy random walk on heterogenous network for MIRNA-disease association prediction,” Mathematical Biosciences, vol. 306, pp. 1–9, 2018.
  * [33] X. Peng and Z. Zhang, “Maximal entropy random walk improves efficiency of trapping in dendrimers,” The Journal of Chemical Physics, vol. 140, no. 23, p. 234104, 2014.
  * [34] M. Smerlak, “Neutral quasispecies evolution and the maximal entropy random walk,” Science Advances, vol. 7, no. 16, p. eabb2376, 2021.
  * [35] R. Sinatra, J. Gómez-Gardeñes, R. Lambiotte, V. Nicosia, and V. Latora, “Maximal-entropy random walks in complex networks with limited information,” Physical Review E, vol. 83, no. 3, p. 030103, 2011.
  * [36] F. Coghi, Large deviation methods for the study of nonequilibrium systems: variational and spectral approaches.  PhD thesis, Queen Mary, University of London, 2021.
  * [37] T. M. Cover and J. A. Thomas, Elements of Information Theory, vol. 2.  Wiley, 66.
  * [38] U. Seifert, “Entropy production along a stochastic trajectory and an integral fluctuation theorem,” Physical Review Letters, vol. 95, no. 4, p. 040602, 2005.
  * [39] H. Touchette, “The large deviation approach to statistical mechanics,” Physics Reports, vol. 478, no. 1-3, pp. 1–69, 2009.
  * [40] F. den Hollander, Large Deviations.  American Mathematical Society, 2000.
  * [41] A. Dembo and O. Zeitouni, Large Deviations Techniques and Applications, vol. 38 of Stochastic Modelling and Applied Probability.  Springer Berlin Heidelberg, 2010.
  * [42] R. Chetrite and H. Touchette, “Nonequilibrium microcanonical and canonical ensembles and their equivalence,” Physical Review Letters, vol. 111, no. 12, p. 120601, 2013.
  * [43] R. Chetrite and H. Touchette, “Nonequilibrium Markov processes conditioned on large deviations,” Annales Henri Poincare, vol. 16, no. 9, pp. 2005–2057, 2015.
  * [44] R. Chetrite and H. Touchette, “Variational and optimal control representations of conditioned and driven processes,” Journal of Statistical Mechanics: Theory and Experiment, vol. 2015, no. 12, p. P12001, 2015\.
  * [45] H. Touchette, “Legendre–Fenchel transforms in a nutshell,” Report Online, 2005.
  * [46] V. S. Borkar, S. Juneja, and A. A. Kherani, “Peformance analysis conditioned on rare events: an adaptive simulation scheme,” Communications in Information and Systems, vol. 3, no. 4, pp. 259–278, 2003.
  * [47] T. P. Ahamed, V. S. Borkar, and S. Juneja, “Adaptive importance sampling technique for Markov chains using stochastic approximation,” Operations Research, vol. 54, no. 3, pp. 489–504, 2006.
  * [48] G. Ferré and H. Touchette, “Adaptive sampling of large deviations,” Journal of Statistical Physics, vol. 172, no. 6, pp. 1525–1544, 2018.
  * [49] G. Ferré and H. Touchette, “Correction to: Adaptive sampling of large deviations,” Journal of Statistical Physics, vol. 173, no. 5, p. 1547, 2018.
  * [50] J. A. Bucklew, Introduction to Rare Event Simulation.  Springer Series in Statistics, Springer New York, 2004.
  * [51] V. S. Borkar, “Asynchronous stochastic approximations,” SIAM Journal on Control and Optimization, vol. 36, no. 3, pp. 840–851, 1998.
  * [52] R. Guimera, S. Mossa, A. Turtschi, and L. N. Amaral, “The worldwide air transportation network: Anomalous centrality, community structure, and cities’ global roles,” Proceedings of the National Academy of Sciences, vol. 102, no. 22, pp. 7794–7799, 2005.
  * [53] L. Gao, W. Wang, L. Pan, M. Tang, and H. F. Zhang, “Effective information spreading based on local information in correlated networks,” Scientific Reports 2016 6:1, vol. 6, no. 1, pp. 1–10, 2016.
  * [54] T. P. Peixoto, “Network Reconstruction and Community Detection from Dynamics,” Physical Review Letters, vol. 123, no. 12, p. 128301, 2019\.
  * [55] P. Holme and J. Saramäki, “Temporal networks,” Physics Reports, vol. 519, no. 3, pp. 97–125, 2012.
  * [56] N. Masuda and R. Lambiotte, A Guide to Temporal Networks, vol. 4 of Series on Complexity Science.  World Scientific (Europe), 2016.
  * [57] P. Holme and J. Saramäki, eds., Temporal Network Theory.  Computational Social Sciences, Cham: Springer International Publishing, 2019.
  * [58] S. Manfredi, E. Di Tucci, and V. Latora, “Mobility and congestion in dynamical multilayer networks with finite storage capacity,” Physical Review Letters, vol. 120, no. 6, p. 068301, 2018.
  * [59] T. Carletti, M. Asllani, D. Fanelli, and V. Latora, “Nonlinear walkers and efficient exploration of congested networks,” Physical Review Research, vol. 2, no. 3, p. 033012, 2020.
  * [60] V. S. Borkar and S. P. Meyn, “Risk-sensitive optimal control for Markov decision processes with monotone cost,” Mathematics of Operations Research, vol. 27, no. 1, pp. 192–209, 2002.
  * [61] A. Basu, T. Bhattacharyya, and V. S. Borkar, “A learning algorithm for risk-sensitive cost,” Mathematics of Operations Research, vol. 33, no. 4, pp. 880–898, 2008.
  * [62] M. Benaïm, Dynamics of stochastic approximation algorithms.  1999\.