Konstantin Avrachenkov
email: k.avrachenkov@sophia.inria.fr

Alexey Piunovskiy

Yi Zhang
email: yi.zhang@liv.ac.uk

Hitting Times in Markov Chains with Restart and their Application to Network Centrality

Keywords: Discrete-time Markov process with restart, Expected hitting time, Network centrality Mathematics Subject Classification (2000) 60J05

Motivated by applications in telecommunications, computer science and physics, we consider a discrete-time Markov process with restart. At each step the process either with a positive probability restarts from a given distribution, or with the complementary probability continues according to a Markov transition kernel. The main contribution of the present work is that we obtain an explicit expression for the expectation of the hitting time (to a given target set) of the process with restart. The formula is convenient when considering the problem of optimization of the expected hitting time with respect to the restart probability. We illustrate our results with two examples in uncountable and countable state spaces and with an application to network centrality.

policy is also used to speedup the Las Vegas type randomized algorithms [START_REF] Alt | A method for obtaining randomized algorithms with small tail probabilities[END_REF], [START_REF] Luby | Optimal speedup of Las Vegas algorithms[END_REF]. Finally, human and animal mobility patterns can be modeled by Markov processes that restart from some locations [START_REF] González | Understanding individual human mobility patterns[END_REF], [START_REF] Walsh | Monkey and cell-phone-user mobilities scale similarly[END_REF].

The main focus of the present work is the expectation of the hitting time (to a given target set) of the process with restart, for which we obtain simple explicit expressions in terms of the expected discounted hitting time of the original process without restart; see Theorem 1. These formulae are useful in the optimization of the expected hitting time with respect to the restart probability. In addition, the formulae allow us to either refine, or give simple and self-contained proofs of the stability results for the process with restart in terms of hitting times. Finally, in Section 3, we illustrate the main results with two examples in uncountable and countable state spaces and with an application to network centrality. In particular, we show that the hitting time based network centrality can be more discerning than PageRank.

Let us mention some related work to the present one in the current literature. The continuous-time Markov process with restart was considered in [START_REF] Avrachenkov | Markov processes with restart[END_REF][START_REF] Evans | Diffusion with stochastic resetting[END_REF][START_REF] Evans | Diffusion with optimal resetting[END_REF][START_REF] Montero | Monotonic continuous-time random walks with drift and stochastic reset events[END_REF]. According to Theorem 2.2 in [START_REF] Avrachenkov | Markov processes with restart[END_REF], the continuous-time Markov process with restart is positive Harris recurrent in case the original process is honest. At the same time, the process with restart is not positive Harris recurrent if the original process is not honest (i.e., the transition kernel is substochastic; in case the state space is countable, that means the accumulation of jumps). The objective of [START_REF] Avrachenkov | Markov processes with restart[END_REF] does not lie in the expected hitting time, but in the representation of the transition probability function of the (continuous-time) process with restart in terms of the one of the original (continuous-time) process without restart. This is trivial in the present discrete-time setup. Here our focus is on the characterization of the expected hitting time. We also would like to mention the two works [START_REF] Dumitriu | On playing golf with two balls[END_REF] and [START_REF] Janson | Hitting times for random walks with restarts[END_REF], dealing with the control theoretic formulation, where the controller decides (dynamically) whether it is beneficial or not to perform a restart at the current state. That line of research can be considered complementary to ours.

The rest of this paper is organized as follows. The description of the process with restart and the main statements are presented in Section 2, which are illustrated by two examples and network centrality application in Section 3. The paper is finished with a conclusion and future research directions in Section 4.

Main statements

Let us introduce the model formally. As in [START_REF] Meyn | Markov chains and stochastic stability[END_REF], let E be a nonempty locally compact Borel state space endowed with its Borel σ-algebra B(E). Consider a discrete-time Markov chain X = { Xt , t = 0, 1, . . . } in the state space E with the transition probability function P (x, dy) being defined by

P (x, Γ ) := pν(Γ ) + (1 -p)P (x, Γ ), (1)
for each Γ ∈ B(E), where p ∈ (0, 1) and P (x, dy) is a transition probability function, and ν is a probability measure on B(E). Let X := {X t , t = 0, 1, . . . } denote the Markov chain corresponding to the transition probability P (x, dy). We assume the two processes X and X are defined on the common probability space (Ω, F, P); if we emphasize that the initial state is x ∈ E, then E x denotes the corresponding expectation operator, and the notation P x is similarly understood. The process X is understood as the modified version of the process X, and is obtained by restarting (independently of anything else) the process X after each transition with probability p ∈ (0, 1) and the distribution of the state after each restart being given by ν; whereas if there is no restart after the transition (with probability 1 -p), the distribution of the post-transition state is P (x, dy) (given that the current state is x).

The following notation is used throughout this paper. Let P t (x, dy), t = 0, 1, . . . , be defined iteratively as follows; for each Γ ∈ B(E),

P 0 (x, Γ ) := I{x ∈ Γ }, P t+1 (x, Γ ) := E P t (x, dy)P (y, Γ ).
The power of other kernels is understood similarly.

Finally, thoughout this paper, the convention 0 • +∞ := 0 is in use.

Known facts

The materials in this subsection are standard and known from [START_REF] Meyn | Markov chains and stochastic stability[END_REF]. The purpose here is to give a short and self-contained presentation. The main result of this paper is postponed to the next subsection.

From (1) it is clear that the process X is Harris recurrent with an irreducibility measure ν. (Recall that the process X is ν-irreducible if for each set Γ ∈ B(E) satisfying ν(Γ ) > 0, it holds that P x (τ Γ < ∞) > 0 for each x ∈ E, where

τ Γ := inf{t = 1, 2, • • • : Xt ∈ Γ }. (2)
As usual, inf ∅ := ∞.) If a Harris recurrent process admits an invariant probability, then it is called positive Harris recurrent; in that case the invariant probability is unique. We verify that X is positive Harris recurrent with the unique invariant probability q given in the next statement.

Proposition 1 The process X is positive Harris recurrent with the unique invariant probability measure q(dy) given by

q(Γ ) = E ∞ t=0 p(1 -p) t P t (y, Γ )ν(dy) (3)
for each Γ ∈ B(E).

Proof. Clearly, q(dy) is a probability measure, and routine calculations verify

q(Γ ) = E q(dx) P (x, Γ ) for each Γ ∈ B(E). 2
We strengthen the above statement in the next corollary. Observe that the process X is aperiodic in the sense of p.118 of [START_REF] Meyn | Markov chains and stochastic stability[END_REF], and the state space E is a petite set since P (x, •) ≥ pν(•) for each x ∈ E. By Theorem 16.0.2 (vi) of [START_REF] Meyn | Markov chains and stochastic stability[END_REF], we see that the process X is uniformly ergodic (see also related arguments in [START_REF] Nummelin | MC's for MCMC'ists[END_REF][START_REF] Nummelin | Geometric ergodicity of Harris recurrent Marcov chains with applications to renewal theory[END_REF]), and satisfies the Doeblin condition, see p.391 of [START_REF] Meyn | Markov chains and stochastic stability[END_REF]. In this connection, the next statement is known as the so called Doeblin's theorem , which is often proved using a coupling argument. However, our model admits a direct simple proof of this fact, and we include it as follows. Below, || • || T V stands for the total variation norm of finite signed measures, i.e., ||µ|| T V := sup A∈B(E) µ(A) -inf A∈B(E) µ(A) for each finite signed measure µ on B(E).

Corollary 1 The process X is uniformly ergodic with the unique invariant probability measure q(•) given in Proposition 1. In particular, we have

|| P n (x, •) -q(•)|| T V ≤ 2(1 -p) n . (4)
Proof. We first note that

P n (x, Γ ) = E n-1 t=0 p(1 -p) t P t (y, Γ )ν(dy) + (1 -p) n P n (x, Γ ),
which can be easily shown by induction. Next, using the above expression and Proposition 1, we can write

|| P n (x, •) -q(•)|| T V = || - E ∞ t=n p(1 -p) t P t (y, •)ν(dy) + (1 -p) n P n (x, •)|| T V ≤ (1 -p) n || - E ∞ t=0 p(1 -p) t P t+n (y, •)ν(dy) + P n (x, •)|| T V ≤ 2(1 -p) n .
Thus, we have established inequality (4) and hence the uniform ergodicity. 2

Hitting times

In this paper we are primarily interested in the expected hitting time of the process X to a given set H ∈ B(E), defined by

η H := inf{t = 0, 1, • • • : Xt ∈ H}.
For the future reference, we put

H P 0 (x, E) := I{x ∈ E \ H}. Denote V (x) := E x [η H ] (5)
and let

H P (x, Γ ) = P (x, Γ \ H)
be the taboo transition kernel with respect to the set H. Then, one can write

V (x) = 1 + p E V (y)ν(dy) + (1 -p) E V (y) H P (x, dy), ∀ x ∈ E \ H, V (x) = 0, ∀ x ∈ H. (6)
Furthermore, it is well known that the function V (x) defined by ( 5) is the minimal nonnegative (measurable) solution to equation ( 6), and can be obtained by iterations

V (n+1) (x) = 1 + p E V (n) (y)ν(dy) + (1 -p) E V (n) (y) H P (x, dy), x ∈ E \ H, n = 0, 1, ... with V (n) (x) = 0 if
x ∈ H for each n = 0, 1, . . . , and V (0) (x) ≡ 0; c.f. e.g., Proposition 9.10 of [START_REF] Bertsekas | Stochastic optimal control: the discrete-time case[END_REF] or Theorem 2 of [START_REF] Silvestrov | Recurrence relations for generalized hitting times for semi-Markov processes[END_REF]. One can actually obtain the minimal nonnegative solution to [START_REF] Avrachenkov | The effect of new links on Google PageRank[END_REF] in the explicit form.

Theorem 1 (a) The minimal nonnegative solution to ( 6) is given by the following explicit form

V (x) = V 1 (x) ∞ t=0 p E V 1 (y)ν(dy) t , ∀ x ∈ E \ H, V (x) = 0, ∀ x ∈ H, (7)
where the function V 1 is given by

V 1 (x) := ∞ t=0 (1 -p) t H P t (x, E), ∀ x ∈ E \ H; V 1 (x) := 0, ∀ x ∈ H. ( 8
)
It coincides with the unique bounded solution to the equation

V 1 (x) = 1 + (1 -p) E V 1 (y) H P (x, dy), ∀ x ∈ E \ H; V 1 (x) = 0, ∀ x ∈ H. ( 9
) (b) If q(H) > 0, then V (x) = V 1 (x) 1 -p E V 1 (y)ν(dy) < ∞, ∀ x ∈ E \ H, V (x) = 0, ∀ x ∈ H. (10)
Proof. (a) Observe that the function V 1 given by ( 8) represents the expected total discounted time before the first hitting of the process X at the set H given the initial state x and the discount factor 1 -p. It thus follows from the standard result about the discounted dynamic programming with a bounded reward that the function V 1 is the unique bounded solution to equation [START_REF] Avrachenkov | Markov processes with restart[END_REF]; see e.g., Theorem 8.3.6 of [START_REF] Hernández-Lerma | Further topics in discrete-time Markov control processes[END_REF]. Now by multiplying both sides of the equation ( 9) by the expression

∞ t=0 p E V 1 (y)ν(dy)
t for all x ∈ E \ H, it can be directly verified that the function V defined in terms of V 1 by ( 7) is a nonnegative solution to [START_REF] Avrachenkov | The effect of new links on Google PageRank[END_REF]. We show that it is indeed the minimal nonnegative solution to (6) as follows.

Let U ≥ 0 be an arbitrarily fixed nonnegative solution to [START_REF] Avrachenkov | The effect of new links on Google PageRank[END_REF]. It will be shown by induction that

U (x) ≥ V 1 (x) n t=0 p E V 1 (y)ν(dy) t ∀ n = 0, 1, . . . , ∀ x ∈ E. (11)
The case when x ∈ H is trivial. Let x ∈ E \ H be arbitrarily fixed. It follows from (6) that

U (x) ≥ 1 + (1 -p) E H P (x, dy) = 1 t=0 (1 -p) t E H P t (x, dy) ≥ 1. ( 12
)
If for some n ≥ 1

U (x) ≥ n t=0 (1 -p) t E H P t (x, dy) , (13)
then by [START_REF] Avrachenkov | The effect of new links on Google PageRank[END_REF],

U (x) ≥ 1 + (1 -p) E U (y) H P (x, dy) ≥ 1 + (1 -p) E n t=0 (1 -p) t E H P t (y, dz) H P (x, dy) = n+1 t=0 (1 -p) t E H P t (x, dy) ,
and so [START_REF] Bonacich | Power and centrality: A family of measures[END_REF] holds for all n ≥ 0 and thus by ( 8)

U (x) ≥ V 1 (x).
Consequently, [START_REF] Bertsekas | Stochastic optimal control: the discrete-time case[END_REF] holds when n = 0. Suppose [START_REF] Bertsekas | Stochastic optimal control: the discrete-time case[END_REF] holds for n, and consider the case of n + 1. Then from (6),

U (x) = 1 + p E U (y)ν(dy) + (1 -p) E U (y) H P (x, dy) ≥ n+1 t=0 p E V 1 (y)ν(dy) t + (1 -p) E U (y) H P (x, dy), (14)
where the inequality follows from the inductive supposition. Define the function W on E by

W (z) = U (z) n+1 t=0 p E V 1 (y)ν(dy) t , ∀ z ∈ E.
Then by [START_REF] Borgatti | Centrality and network flow[END_REF],

W (x) ≥ 1 + (1 -p) E H P (x, dy) = 1 t=0
(1 -p) t E H P t (x, dy) ; c.f. [START_REF] Bianchini | Inside PageRank[END_REF]. Now, based on ( 14), a similar reasoning by induction as to the verification of ( 13) for all n ≥ 0 shows that

W (x) ≥ k t=0 (1 -p) t E H P t (x, dy) , ∀ k = 0, 1, . . . ,
and thus W (x) ≥ V 1 (x). This means U (x) ≥ V 1 (x) n+1 t=0 p E V 1 (y)ν(dy) t .
Thus by induction, [START_REF] Bertsekas | Stochastic optimal control: the discrete-time case[END_REF] holds for all n ≥ 0, and thus

U (x) ≥ V (x)
by ( 7), as desired.

(b) If q(H) > 0, then there exists some T ≥ 0 such that E P T (y, H)ν(dy) > 0, meaning that there exists some T ≤ T such that E H P T (x, E)ν(dx) < 1. Thus,

0 ≤ p E ∞ t=0 (1 -p) t H P t (y, E)ν(dy) < 1,
and so the geometric series

∞ t=0 p E V 1 (y)ν(dy) t converges. The statement follows. 2
The next corollary is immediate.

Corollary 2 If q(H) > 0, then both V 1 (x) and V (x) are bounded with respect to the state x ∈ E. In particular, we have

V 1 (x) ≤ 1 p .
There is a nice probabilistic interpretation of the decomposition presented in Theorem 1. Suppose q(H) > 0, which is case (b) in Theorem 1. Let us make a change of variable

ρ(x) = 1 -pV 1 (x) = 1 -p ∞ t=0 (1 -p) t H P t (x, E) = P x [hit before restart].
Then, equation [START_REF] Baumgärtel | Analytic perturbation theory for matrices and operators[END_REF] takes the form

V (x) = 1 -ρ(x) p E ρ(y)ν(dy) = P x [no hit before restart] pP ν [hit before restart] .
Let us note that the moments of restart determine regenerative cycles [START_REF] Feller | An introduction to probability theory and its applications: Volume I[END_REF]. The expected number of steps in a cycle from one restart to another is 1/p, and the expected number of cycles until we hit the set H is roughly 1/P ν [hit before restart]. Thus, the total expected steps is approximately equal to 1/(pP ν [hit before restart]). This is only an approximation since we did not carefully take into account the impact of the first and the last incomplete cycles. It turns out that the factor P x [no hit before restart] gives the necessary correction. Theorem 1 can be viewed as a generalization of the results of [START_REF] Avrachenkov | The effect of new links on Google PageRank[END_REF][START_REF] Hopcroft | Manipulation-resistant reputations using hitting time[END_REF] to the general state space from the finite state space and to arbitrary initial state. We also have the following chain of equivalent statements.

Theorem 2

The following statements are equivalent. (a) q(H) > 0.

(b) V (x) = E x [η H ] < ∞ for each x ∈ E. (c) V (x) = E x [η H ] < ∞ for almost all x ∈ E with respect to q(dy). (d) sup x∈E V (x) = sup x∈E E x [η H ] < ∞.
The proof of this theorem follows from the next lemma.

Lemma 1 Statements (a) and (c) in Theorem 2 are equivalent.

Proof. Statement (a) implies statement (c) by Theorem 1(b). It remains to show that statement (c) implies statement (a). To this end, suppose for contradiction that statement (a) does not hold, i.e., q(H) = 0. Then it follows from (3) that

E ∞ t=0 P t (y, H)ν(dy) = 0.
Thus, there exists a measurable subset Γ of E \ H such that

ν(Γ ) = 1 (15)
and

P x [η X H = ∞] = 1 ∀ x ∈ Γ, (16)
where

η X H := inf{t = 0, 1, • • • : X t ∈ H}. (17)
Let x ∈ Γ be fixed. Then by ( 8) and ( 16)

V 1 (x) = E x   η X H -1 t=0 (1 -p) t   = E x ∞ t=0 (1 -p) t = 1 p .
Since x ∈ Γ is arbitrarily fixed, and ν(Γ ) = 1 by ( 15), we see from [START_REF] Avrachenkov | A singular perturbation approach for choosing the PageRank damping factor[END_REF] and the fact that . However, we gave its self-contained proof in the above without referring to the result concerning the general Markov chains in [START_REF] Meyn | Markov chains and stochastic stability[END_REF].

V 1 (y) ≥ 1 if y ∈ E \ H that V (y) = ∞ for each y ∈ E \ H. Since q(E \ H) = 1,

Optimization problem

Next let us consider the dependance of V (x) on the restart probability p for each fixed x ∈ E. When we emphasize the dependance on p, we explicitly write V (x, p) and V 1 (x, p). In the above, p ∈ (0, 1) was fixed. Now we formally put

V (x, 0) := ∞ t=0 H P t (x, E) ∈ [1, ∞], ∀ x ∈ E \ H; V (x, 0) := 0, ∀ x ∈ H,
which represents the expected hitting time of the process without restart, and

V (x, 1) := 1 ν(H) ∈ [1, ∞], ∀ x ∈ E \ H; V (x, 1) := 0, ∀ x ∈ H,
which represents the expected hitting time of the process that restarts with full probability at each transition.

Here and below, c 0 := ∞ for any c > 0. As usual, the continuity of V (x, p) at p = a means lim p→a

V (x, p) = V (x, a) ∈ [-∞, ∞]. Theorem 3 Suppose E ∞ t=0 P t (y, H)ν(dy) > 0. ( 18
)
Then the function V (x, p) is infinitely many times differentiable in p ∈ (0, 1) and is continuous in p ∈ [0, 1]. As a consequence, the problem

Miniminize V (x, p) with respect to p ∈ [0, 1] (19)
is solvable.

Remark 2

The condition [START_REF] Engström | PageRank for networks, graphs and Markov chains[END_REF] in the above statement is equivalent to q(H) > 0 for some and then all p ∈ (0, 1) by ( 3).

Proof of Theorem 3. If x ∈ H, the statement holds trivially since V (x, p) = 0 for each p ∈ [0, 1]. Consider now

x ∈ E \ H. One can see that V 1 (x, p) (as given by ( 8)) and (1 -p E V 1 (y, p)ν(dy)) -1 are both infinitely many times differentiable in p ∈ (0, 1). It follows from this and (10) that V (x, p) is infinitely many times differentiable in p ∈ (0, 1). For the continuity of V (x, p) at p = 0, it holds that

1 -p E V 1 (y, p)ν(dy) = 1 -p E E y   η X H -1 t=0 (1 -p) t   ν(dy) = 1 - E E y 1 -(1 -p) η X H ν(dy) = E E y (1 -p) η X H ν(dy) → 1,
as p → 0, by the monotone convergence theorem; recall [START_REF] Dumitriu | On playing golf with two balls[END_REF] for the definition of η X H . It follows from this fact, ( 8), [START_REF] Baumgärtel | Analytic perturbation theory for matrices and operators[END_REF] and the monotone convergence theorem that

lim p→0 V (x, p) = lim p→0 V 1 (x, p) lim p→0 1 -p E V 1 (y, p)ν(dy) = V (x, 0)
as desired. The continuity of the function V (x, p) at p = 1 can be similarly established. The last assertion is a well known fact; see e.g., [START_REF] Bertsekas | Stochastic optimal control: the discrete-time case[END_REF].

2
We shall illustrate the optimization problem [START_REF] Evans | Diffusion with stochastic resetting[END_REF] by two examples in the next section.

3 Numerical examples and application

Uni-directional random walk on the line

Let E = R, H = [a, b]
with a < b being two real numbers. The process X only moves to the right, and the increments of each of the transitions are i.i.d. exponential random variables with the common mean 1 µ > 0. The restart probability is denoted as p ∈ (0, 1) as usual, and the restart distribution ν is arbitrary. Below by using Theorem 1 we provide the explicit formula for the expected hitting time at H of the restarted process X. (Clearly, if the initial state is outside H, then the expected hitting time of the process X at the set H is infinite.)

We can give the following informal description of this example. There is a treasure hidden in the interval [a, b] and one tries to find the treasure. Once the searcher checks one point in the interval [a, b], he finds the treasure. The searcher has the means only to stop and to check points between the exponentially distributed steps. This models the cost of checking frequently. It is also natural to restart the search from some base. Intuitively, by restarting too frequently, the searcher spends most of the time near the base and does not explore the area sufficiently. On the other hand, restarting too seldom leads the search to very far locations where the searcher spends most of the time for nothing. Hence, intuitively there should be an optimal value for restarting probability.

One can verify that in this example the unique bounded solution to ( 9) is given by

V 1 (x) = 0 for each x ∈ [a, b], V 1 (x) = 1 p (20)
for each x > b,

V 1 (x) = 1 p - 1 -p p 1 -e -µ(b-a) e -µ(a-x)p (21)
for each x < a. In fact, this can be conveniently established using the following probabilistic argument. Recall that V 1 (x) represents the expected total discounted time up to the hitting of the set H by the process X; see [START_REF] Avrachenkov | Beta current flow centrality for weighted networks[END_REF]. So for [START_REF] Evans | Diffusion with optimal resetting[END_REF] one merely notes that with the initial state x > b, η X H = ∞, where η X H is defined by [START_REF] Dumitriu | On playing golf with two balls[END_REF]. For [START_REF] Feller | An introduction to probability theory and its applications: Volume I[END_REF], one can write for each x < a that

V 1 (x) = E x   η X H -1 t=0 (1 -p) t   = E x   E x   η X H -1 t=0 (1 -p) t η X H     .
Now the expected hitting time of the restarted process X to the set H = [a, b] is given by

V (x) = 1 p - 1 -p p 1 -e -µ(b-a) e -µ(a-x)p ∞ t=0 p E V 1 (y)ν(dy) t ,
with the initial state x < a, and by

V (x) = 1 p ∞ t=0 p E V 1 (y)ν(dy) t ,
with the initial state x > b, recall [START_REF] Avrachenkov | A singular perturbation approach for choosing the PageRank damping factor[END_REF]. If the restart distribution ν is not concentrated on (b, ∞), then q(H) > 0, and by Theorem 1(b) we have

V (x) = 1 p -1-p p 1 -e -µ(b-a) e -µ(a-x)p 1 -p E V 1 (y)ν(dy) ,
with the initial state x < a and by

V (x) = 1 p 1 -p E V 1 (y)ν(dy) ,
for each x > b. In particular, if the process restarts from a single point r < a, the above expressions can be specified to

V (x) = 1 -(1 -p)(1 -e -µ(b-a) )e -µ(a-x)p p(1 -p)(1 -e -µ(b-a) )e -µ(a-r)p ,
for the initial state x < a and to

V (x) = 1 1 -e -µ(b-a) 1 p(1 -p)e -µ(a-r)p ,
for each x > b. For the latter case (x > b), by standard analysis of derivatives, one can find the optimal value of the restart probability minimizing the expected hitting time of the process with restart in a closed form, as given by

p opt = 2 2 + µ(a -r) + 4 + µ 2 (a -r) 2 .
Now we can make several observations: the first somewhat interesting observation is that in the case when the initial state is to the right of the interval [a, b], the value of the optimal restart probability does not depend on the length of the interval but only on the average step size and on the restart position. The second observation is that when µ(a -r) is small, i.e., when either the average step size is large or the restart position is close to H, the optimal restart probability is close to 1/2. Thirdly, when µ(a -r) is large, the optimal restart probability is small and reads

p opt = 1 1 + µ(a -r)
+ o 1 µ(a -r) .

Random walk on the one dimensional lattice

Let us consider a symmetric random walk on the one dimensional lattice which aims to hit H = {0} with restart at some node r. Assume without loss of generality that the restart state r is on the positive half-line, i.e., r > 0.

From Theorem 1 we conclude that it is sufficient to solve the following equations

V 1 (k) = 1 + 1 -p 2 [V 1 (k -1) + V 1 (k + 1)], k = 0, V 1 (0) = 0.
Following the standard approach for solution of difference equations, we obtain

V 1 (k) = cα k 1 + 1 p ,
where α 1 < 1 is the minimal solution to the characteristic equation

α = 1 -p 2 [1 + α 2 ],
and the constant c = -1 p comes from the condition V 1 (0) = 0. Consequently,

V (k) = V 1 (k) 1 -pV 1 (r) = 1 -α k 1 pα r 1 .
An elegant analysis can be done for the limiting case when the initial position k goes large, and hence we now minimize lim k→∞ V (k) = 1/(pα r 1 ), or equivalently, maximize pα r 1 with respect to p ∈ (0, 1). This leads to the following equation for the optimal restart probability p α 1

dα 1 dp = - 1 r . (22)
Indeed, we note that

lim p→0 p α 1 dα 1 dp = 0, lim p→1 p α 1 dα 1 dp = -∞, and
d dp p α 1 dα 1 dp = - 1 1 -(1 -p) 2 1 -(1 -p 2 )(1 -p) (1 -p) 2 (1 -(1 -p) 2 ) < 0.
Thus, the left hand side of ( 22) is a monotone function decreasing from zero to minus infinity. Consequently, the unique solution of equation ( 22) is the global minimizer of 1/(pα r 1 ). The equation ( 22) can be transformed to the polynomial equation 1

r 2 (1 -p) 2 (2 -p) = p
Consider the case of large r. This is a so-called case of singular perturbation, as the small parameter 1/r 2 is in front of the largest degree term, [START_REF] Avrachenkov | Analytic perturbation theory and its applications[END_REF], [START_REF] Baumgärtel | Analytic perturbation theory for matrices and operators[END_REF]. It is not difficult to see that for large values of r, the equation has one real root that can be expanded as

p opt = c 1 r 2 + c 2 r 4 + ... (23)
and two complex roots that move to infinity as r → ∞. By substituting the series (23) into the polynomial equation, we can identify the terms c i , i = 1, 2, .... Thus, we obtain p opt = 2 r 2 -10 r 4 + o(r -4 ).

Application to network centrality

One of the main tasks in network analysis is to determine which nodes are more "central" than the others. Node degree and PageRank [START_REF] Brin | The anatomy of a large-scale hypertextual Web search engine[END_REF] are examples of widely used centrality measures. An interested reader can find comprehensive discussions on various aspects of network centrality in [2, 12-15, 23, 30, 31]. We note that PageRank is a stationary distribution of the random walk with restart and in our setting it is just measure q. Both node degree and PageRank are prone to manipulation or so-called "sybil attack". To mitigate this problem, the authors of [START_REF] Hopcroft | Manipulation-resistant reputations using hitting time[END_REF][START_REF] Liu | Personalized hitting time for informative trust mechanisms despite sybils[END_REF] proposed hitting time based centrality measures. Here we show that the hitting time based centrality can be more discerning. Let us first consider a simple 6-node network with weighted edges, see Figure 1. The weights are depicted near the edges. In Table 1 we give centrality values of the nodes with respect to node degree, PageRank and Expected Hitting Times starting and restarting both from the uniform distribution. The last row in Table 1 lists down the expected hitting times to the target state A, B, and so on. It is intuitively clear that nodes A and D are more central in this network. However, both node degree and PageRank indicate equal importance for the nodes. In contrast, the hitting time based centrality clearly indicates that nodes A and D are more central than the other nodes.

Let us now consider an example of a real social network. The example was taken from online social network VKontakte and represents a principal component of the interest group about Game Theory [START_REF] Avrachenkov | Beta current flow centrality for weighted networks[END_REF]. The example has 71 nodes and 116 weighted edges, see Figure 2 (taken from [START_REF] Avrachenkov | Beta current flow centrality for weighted networks[END_REF]). The edge weight is equal to the number of common friends. Only the edges with a weight more than two have been kept.

In Table 2 we provide top-10 lists of nodes according to node degree, PageRank and the expected hitting time. We observe that the top-10 list by PageRank has 9 nodes from the top-10 list by node degree. The top-10 list by the expected hitting time has only 5 nodes from the top-10 list by node degree. Note that nodes 22 and 13, which intuitively look quite central, are not in the top-10 list by PageRank.

Conclusion

In conclusion, in this paper we present a self-contained study of a discrete-time Markov process with restart. Our primary interest is in the expected hitting time of the process with restart to a target set. We obtained the formula of the expected hitting time of the restarted process to a target set, and considered the optimization problem of the expected hitting time with respect to the restart probability. We illustrated our results with two examples in uncountable and countable state spaces and one application to network centrality. In particular, we show that the network centrality based on hitting times is more selective. We think that it would be interesting to study perturbation and optimization of the network centrality based on hitting times along similar analysis for PageRank [START_REF] Avrachenkov | A singular perturbation approach for choosing the PageRank damping factor[END_REF][START_REF] Bianchini | Inside PageRank[END_REF][START_REF] Engström | PageRank for networks, graphs and Markov chains[END_REF][START_REF] Fercoq | Ergodic control and polyhedral approaches to PageRank optimization[END_REF][START_REF] Langville | Deeper inside PageRank[END_REF]. Our general results may also have potential application to network community analysis, which we intend to explore in the future.

2 Proof of Theorem 2 . 2 Remark 1

 2221 we see that statement (c) does not hold. From Theorem 1(b), we see that (a) implies (b), which implies (c). From Lemma 1, (c) implies (a). Clearly, (d) implies (b). Finally, (a) implies (d) because V 1 (x) is bounded; see Corollary 2. Corollary 1 and Theorem 16.0.2 of [35] immediately give the relation (a) ⇒ (d) ⇒ (b) ⇒ (c) in Theorem 2

Fig. 1

 1 Fig. 1 Example of 6-node network with weighted edges.

  time, p = 0.15 6.28 8.28 8.28 6.28 8.28 8.28

Fig. 2

 2 Fig. 2 Example of a social network.

FinallyFig. 3

 3 Fig.3The expected hitting time from A to B.

Table 1

 1 Centralities for 6-node network.

Table 2

 2 Top-10 lists for the social network example.

	Node degree	1 8 4 20 6 56 7 28 44 32
	PageRank, p = 0.15	1 8 56 28 44 4 32 20 63 6
	Hitting time, p = 0.15 1 8 56 28 63 22 13 33 69 4

Acknowledgements

This work was partially supported by the European Commission within the framework of the CONGAS project FP7-ICT-2011-8-317672. Y.Zhang's work was carried out with a financial grant from the Research Fund for Coal and Steel of the European Commission, within the INDUSE-2-SAFETY project (Grant No. RFSR-CT-2014-00025). This is the author version of the paper published in Methodology and Computing in Applied Probability, v.20(4), pp.1173-1188, 2018.


A Stamm

O Commowick

A Menafoglio

S K Warfield

A Bayes Hilbert Space for Compartment Model Computing in Diffusion MRI

The single diffusion tensor model for mapping the brain white matter microstructure has long been criticized as providing sensitive yet non-specific clinical biomarkers for neurodegenerative diseases because (i) voxels in diffusion images actually contain more than one homogeneous tissue population and (ii) diffusion in a single homogeneous tissue can be non-Gaussian. Analytic models for compartmental diffusion signals have thus naturally emerged but there is surprisingly little for processing such images (estimation, smoothing, registration, atlasing, statistical analysis). We propose to embed these signals into a Bayes Hilbert space that we properly define and motivate. This provides a unified framework for compartment diffusion image computing. Experiments show that (i) interpolation in Bayes space features improved robustness to noise compared to the widely used log-Euclidean space for tensors and (ii) it is possible to trace complex key pathways such as the pyramidal tract using basic deterministic tractography thanks to the combined use of Bayes interpolation and multi-compartment diffusion models.

Introduction

Diffusion MRI is a unique in-vivo and non-invasive imaging technique that probes the microstructure of tissues by tracking water diffusion [START_REF] Le Bihan | Looking into the functional architecture of the brain with diffusion MRI[END_REF]. The diffusion process is a probability distribution that describes the random 3-dimensional displacements of water molecules due to thermal agitation. The signal recovered in diffusion MRI is the Fourier transform of the diffusion density evaluated at the magnetic field spatial gradient to which the imaged brain was subjected.

A practical approach to access the diffusion density pertains to resorting to parametric models that lead to an analytic expression of the observed diffusion signals. The first parametric model was introduced two decades ago in [START_REF] Basser | Estimation of the effective self-diffusion tensor from the NMR spin echo[END_REF] and is called single tensor (ST) model. The ST model essentially assumes that the diffusion in the voxel is well characterized by a zero-mean Gaussian distribution with covariance matrix proportional to the so-called diffusion tensor. It has become widely used in clinical practice, which raised the demand for a solid mathematical framework for processing tensor images (estimation, smoothing, registration, atlasing, statistical analysis). As a response, several frameworks have been proposed, including a Riemannian affine-invariant framework [START_REF] Pennec | A riemannian framework for tensor computing[END_REF][START_REF] Lenglet | Statistics on the manifold of multivariate normal distributions: Theory and application to diffusion tensor MRI processing[END_REF] and the popular log-Euclidean space for tensors, embedding them in a Lie group structure [START_REF] Arsigny | Log-Euclidean metrics for fast and simple calculus on diffusion tensors[END_REF].

However, because of the low spatial resolution in diffusion MRI, it has been shown that the ST model is an average of several diffusion processes arising from multiple populations of tissues in the voxel, ultimately leading to an inaccurate description of the microstructure in most parts of the brain white matter. In addition, even within a single homogeneous tissue population, non-Gaussian diffusion has been observed [START_REF] Wang | When Brownian diffusion is not Gaussian[END_REF]. This has led to the extension of the ST model to mixture models, often called multi-compartment models (MCM), where the voxelwise diffusion signal is modeled as a linear combination of compartmental diffusion signals arising from underlying homogeneous diffusion processes [START_REF] Panagiotaki | Compartment models of the diffusion MR signal in brain white matter: a taxonomy and comparison[END_REF]. To the best of our knowledge, the only mathematical framework available for multicompartment image processing [START_REF] Taquet | A mathematical framework for the registration and analysis of multi-fascicle models for population studies of the brain microstructure[END_REF] is however limited to multi-tensor images (mixture of single-tensor signals, see [START_REF] Scherrer | Parametric representation of multiple white matter fascicles from cube and sphere diffusion MRI[END_REF]) as it relies on log-Euclidean geometry.

In this work, we present an alternative to the log-Euclidean space on tensors called the Bayes Hilbert space [START_REF] Van Den Boogaart | Bayes Hilbert Spaces[END_REF] for processing compartmental diffusion signals instead of tensors. It provides a unified framework that can accommodate analytic compartment diffusion signals of any type (Gaussian and non-Gaussian). Section 2 gives a general introduction to Bayes Hilbert spaces with motivation and setup in diffusion MRI. Section 3 describes a simulation study and a tractography application on real data to compare Euclidean, log-Euclidean and Bayes interpolation. Results commented in Section 4 show that Bayes interpolation features improved robustness to noise and yields better pyramidal tract reconstructions when combined with MCM-based deterministic tractography to account for streamline atlas priors when diffusion models cannot be locally trusted.

Theory

Bayes Hilbert Spaces and Compartmental Diffusion

The signal in diffusion MRI, hereafter diffusion signal, observed after the application of a diffusion gradient q, where q 2 = b is the well known b-value, undergoes a signal decay starting from a baseline signal. The form of this decay depends on how water molecules diffuse in tissues in the vicinity of the spatial location where the signal is observed. In general, in a voxel composed of K different tissue structures, the compound diffusion signal reads:

S(q) = K j=1 S j (q) = K j=1 S 0j A j (q), S 0j > 0, A j ∈ [0, 1], (1)
in which all the information pertaining to the microstructure (set of surrounding tissue structures) is grasped by the signal attenuations A j 's. This is the signal decomposition at the foundation of diffusion compartment imaging [START_REF] Panagiotaki | Compartment models of the diffusion MR signal in brain white matter: a taxonomy and comparison[END_REF]. Hence, up to some multiplicative constant, the compartmental diffusion signal S j can be interpreted as the density of a probability measure ν j on the space of diffusion gradients, which is absolutely continuous w.r.t. the Lebesgue measure λ such that dνj dλ (Radon-Nykodym derivative w.r.t. the Lebesgue measure λ) is proportional to S j . Furthermore, this multiplicative constant is related to the T 2 relaxation time of tissue j and thus not relevant for microstructure mapping.

To account for this, we propose to embed the probability measures induced by compartmental diffusion signals into the Bayes space B 2 (P ) of (classes of equivalence of B(R 3 )-finite) measures on the Borel space (Ω, B(Ω)), with squareintegrable log-densities w.r.t. a reference measure P [START_REF] Van Den Boogaart | Bayes Hilbert Spaces[END_REF], which reads:

B 2 (P ) = ν : R 3 log 2 dν dP dP < ∞, ν > 0 . (2)
In this view, if a diffusion signal S j1 (or, equivalently a probability measure ν 1 induced by S j1 ) carries a piece of information about the structure of tissue j, then a diffusion signal S j2 proportional to S j1 does not carry additional information about the structure of tissue j. Hence, S j1 and S j2 are regarded as equivalent for microstructure mapping. This key property, known as scale invariance, is accounted for in the Bayes space B 2 (P ) by the induced classes of equivalence. Not every Bayes space B 2 (P ) however can accommodate compartmental diffusion signals because, depending on the choice of the reference measure P , the logarithm of compartmental diffusion signals might not be square-integrable. The choice of the reference measure of the Bayes space is thus critical. Provided that an appropriate reference measure P exists, embedding diffusion signals in B 2 (P ) provides a unified framework for processing any type of analytic microstructure compartment models [START_REF] Panagiotaki | Compartment models of the diffusion MR signal in brain white matter: a taxonomy and comparison[END_REF] at the cost of numerical integrations. The space B 2 (P ) is a vector space when endowed with the perturbation and powering operations (⊕, ), defined for ν 1 , ν 2 ∈ B 2 (P ), A ∈ B(R 3 ), α ∈ R as:

(ν 1 ⊕ ν 2 )(A) = B(P ) A dν 1 dP dν 2 dP dP, (α ν 1 )(A) = B(P ) A dν 1 dP α dP,
and becomes a separable Hilbert space when equipped with a proper inner product [START_REF] Egozcue | Bayes spaces: use of improper distributions and exponential families[END_REF]. In the following subsection, we focus on the special case of Gaussian compartmental diffusion for which analytic expressions can be obtained for the perturbation and powering operations as well as for the distance.

Gaussian Compartmental Diffusion

The diffusion signal arising from Gaussian compartmental diffusion reads [START_REF] Basser | Estimation of the effective self-diffusion tensor from the NMR spin echo[END_REF]:

S(q) = S 0 e -q Dq , ( 3
)
where D is the so-called diffusion tensor, i.e. a 3 × 3 symmetric positive definite (SPD) matrix. Hence, the probability measure ν induced by S is a Gaussian distribution with mean 0 and covariance matrix D -1 /2 and we can write, for the purpose of notation, S ∼ N (0, D -1 /2). As a result, Gaussian diffusion signals cannot be embedded in B 2 (λ) because choosing the Lebesgue measure as reference makes the integral of the squared logarithm of S diverge. However since all moments of the Gaussian distribution are finite, if we define the reference measure P as a Gaussian distribution with mean 0 and covariance matrix D ref SPD, i.e. P ∼ N (0, D ref ), then Gaussian diffusion signals can be embedded in B 2 (P ). Furthermore, it is possible to obtain analytic expression for the operations ⊕ and . The perturbation of measure ν 1 by measure ν 2 reads:

(ν 1 ⊕ ν 2 )(A) = B A dν 1 dλ dν 2 dλ dλ dP dλ = B A e -q (D1+D2-D ref )q dλ(q).
Hence, (ν

1 ⊕ ν 2 ) ∼ N (0, D ref + (D 1 -D ref ) + (D 2 -D ref )
). The result is another Gaussian distribution whose covariance matrix is centered around the one of the reference Gaussian distribution but perturbed by the covariance information in D 1 and D 2 that are not already present in the reference covariance structure.

The multiplication of measure ν by a scalar α ∈ R reads:

(α ν)(A) = B A dν dλ α dλ dP α-1 dλ = B A e -q (αD+(1-α)D ref )q dλ(q).
Hence, (α ν) ∼ N (0, αD

+ (1 -α)D ref ).
The result is another Gaussian distribution whose covariance matrix is a linear combination between the covariances of ν and P . This offers the possibility of giving more or less weight to the information content of a diffusion signal in terms of microstructure. If one think of α ∈ [0, 1], then (α ν) means that less weight is put on the information content of ν (which could come from a poor model estimate, low SNR, etc) in favor of the information content of the reference P . In practice, it shrinks the diffusion tensor of the signal towards the one of the reference. Note that writing the results of the operations of ⊕ and as Gaussian distributions is an abuse of notation. In effect, the space of diffusion signals is a subspace of B 2 (P ) that is not closed in the Bayes geometry. This is because operations on diffusion signals in B 2 (P ) might yield non positive "covariance" matrices. However, most applications in which analytic diffusion model computing is required involve weighted average operations, i.e. linear combinations with positive weights whose sum is less than one, which always produce SDP tensors in B 2 (P ). In effect, it is easy to show that:

N i=1 w i ν i ∼ N 0, N i=1 w i D i + 1 - N i=1 w i D ref . (4)
In the rare events that operations in Bayes space produce non-positive matrices, one can perform an orthogonal projection of the non-positive matrix back into the space of diffusion signals by solving the following minimization problem:

min D∈B(P ) s.t. D>0 ad 2 B(P ) (M, D) + (1 -a)d 2 B(P ) (D ref , D), a ∈ [0, 1], ( 5
)
where M is the symmetric matrix assumed to have some negative eigenvalues and d B(P ) is the distance on B 2 (P ) given by:

ν 1 ν 2 2 B(P ) := Tr (D 1 -D 2 )D -1 ref (D 1 -D 2 )D -1 ref = (D 1 -D 2 )D -1 ref 2 F .
3 Experimental Setup

Simulations: Robustness to Noise

The goal of the simulated data is to assess the robustness of Bayes space interpolation to MRI-induced noise, compared to more traditional approaches that focus directly on the tensor and embed it either in Euclidean or log-Euclidean space.

In Bayes space, as shown in Eq. ( 4), the tensor associated with the interpolated signal is a convex combination of the tensor interpolated in Euclidean space and the tensor of the reference signal, where the weight w data associated with the tensor interpolated in Euclidean space shall indicate how much we are willing to trust its information content. For this purpose, we set w data := 1 -e -SNR/β and include in the comparison three Bayes spaces (Bayes10, Bayes20, Bayes30), using β = 10, 20, 30 respectively. We generate (i) a set of 31 noise-free diffusion signals according to Eq. ( 3) with a baseline signal S 0 = 1 and 31 diffusion gradients q uniformly distributed on the hemisphere of radius √ b, with b = 1000s/mm 2 and (ii) a set of n = 8 normalized weights to define spatial weights of the 8 neighbors. Next, for a given SNR and a given method, we obtain a Monte-Carlo estimate of the mean-square error (MSE) of the interpolated tensor by averaging the squared distances between the ground truth tensor and R = 1000 replicates of interpolated tensors from noisy neighboring tensors produced as follows:

(a) Add Rician noise to the noise-free signals using σ = S 0 /SNR; (b) Get a realistic noisy tensor field by estimating a diffusion tensor in each neighboring voxel via maximum likelihood estimation [START_REF] Stamm | Comprehensive Maximum Likelihood Estimation of Diffusion Compartment Models Towards Reliable Mapping of Brain Microstructure[END_REF];

(c) Interpolate the tensor in the central voxel from the tensors in the 8 neighboring voxels to which we associate the initially simulated spatial weights;

(d) Compute distance between the interpolated and ground truth tensor. We used 8 SNR ∈ [6.25, 50] and three metrics to compare interpolations with ground truth: the angle between principal orientations (direction distance), the Euclidean distance between axial diffusivities (diffusivity distance) and the Euclidean distance between radial diffusivities (radius distance as it is often used as a proxy to measure axon radii). We chose these metrics as they focus on microstructural parameters of direct clinical relevance and are independent from the compared spaces. For interpolation in Bayes space, for each SNR, we set the reference tensor as a noisy version of the ground truth tensor using the same procedure as above using SNR √ N . This is because, often, the available data to use as reference comes from atlases generated on similar data but averaged on multiple subjects (N ). Hence, there is still uncertainty in the reference tensor but its variance is likely to be divided by N . In this simulation, we set N = 20.

Case Study: Tractography of the Pyramidal Tract

The pyramidal tract (PT) is of primary importance as it handles volontary motion. Neurons initiate in the primary motor cortex (R0), goes successively through the corona radiata (R1), genu and posterior limb of the internal capsule (R2) and cerebral peduncles (R3) to eventually enter the spinal cord (R4). The PT is difficult to reconstruct from tractography due to the large spreading of neurons on the cortex. We aim at showing that interpolation in Bayes space makes deterministic FACT tractography [START_REF] Mori | Three-dimensional tracking of axonal projections in the brain by magnetic resonance imaging[END_REF] feasible for PT reconstruction. In details, we compare four approaches obtained combining single-or multi-tensor FACT tractography (ST or MT) with log-Euclidean or Bayes interpolation (ST-LogEuclidean, ST-Bayes, MT-LogEuclidean and MT-Bayes).

We scanned a healthy subject on a 3T Siemens Verio magnet for a T 1 MPRAGE image at 1mm 3 resolution and diffusion data at 2mm 3 resolution using the same gradient table as in the simulations. We used an available PT atlas in MNI space based on 20 healthy subjects that underwent the same acquisition protocol. We estimated ST and MT models from the diffusion data [START_REF] Stamm | Comprehensive Maximum Likelihood Estimation of Diffusion Compartment Models Towards Reliable Mapping of Brain Microstructure[END_REF] and brought the resulting images into MNI space where FACT tractography was performed. In the case of MT-FACT, at a given point, a single tensor was picked from each neighboring mixture according to the highest orientation similarity w.r.t. arrival direction so that we could proceed as in ST-FACT. We performed tractography of the PT following regions of interest (ROI). We used R0 as seeding mask, stopped a streamline generation when FA fell below 0.1 (for ST-FACT) or when linearly interpolated fraction of free water exceeded 0.8 (for MT-FACT) and filtered the resulting streamlines progressively through R1, R2, R3 and R4. We defined the reference weight at position x for Bayes interpolation as:

w ref (x, FA, SNR, x ref ; β, δ) := 1 -FA(1 -e -SNR/β ) e -x-x ref /δ ,
where FA is the fractional anisotropy of the tensor interpolated in Euclidean space, x ref is the position of the closest point on the PT atlas and (β, δ) are user-defined parameters that control the decay velocity of SNR and distance weights. In essence, we put more trust in the data when both SNR and orientation coherence between neighbors are high or when the reference is too far.

Results & Discussion

Simulations: Robustness to Noise. Figure 1 shows the MSE between the ground truth tensor in central voxel and the interpolated tensor from noisy neighbors, with increasing amount of noise. For all metrics, interpolation in Bayes space performs uniformly better in terms of robustness to MRI-induced noise for recovering microstructural parameters. The MSE curves for the three Bayes spaces that give increasing importance to the reference measure (from Bayes10 to Bayes30) reveals that even when we trust mostly the data (Bayes10), interpolation in Bayes space is preferable.

Case Study: Tractography of the Pyramidal Tract. Figure 2 shows PT reconstruction with increasing number of filtering ROIs (from left to right) for all four methods. First, observe that only methods based on Bayes interpolation (rows 2 and 4) successfully manage to reconstruct streamlines that go through all four ROIs and thus are more likely to belong to the PT. Also, in general, deterministic FACT tractography fails to recover PT streamlines with the traditional log-Euclidean interpolation (rows 1 and 3) . This is well documented in the literature on ST-FACT. We hypothesize that MT-FACT has a lot more directional information for tracking but, without prior anatomical knowledge, the FACT algorithm (which stepwise follows the most collinear direction) has no mechanism to channel streamlines into following the PT shape. This explanation is supported by the sequential filtering that shows high streamline variability when filtering only by the corona radiata (1st column, 3rd row) and almost no remaining streamlines at the end of the entire filtering process (4th column, 3rd row). The ST-Bayes method seems to mainly follow the shape of atlas streamlines. This is due to an inherent model mis-specification in most parts of the white matter where the ST model provides an insufficient description of the microstructure and therefore presents an articially low FA which uniformly inflates the weight of the reference measure. Conversely, the MT-Bayes version nicely preserves the PT streamlines after complete filtering without heavy influence of the atlas since the MT model provides an accurate description of the microstructure and thus the interpolated tensor has low FA only when neighboring tensors have heterogeneous orientations.

Fig. 1 .

 1 Fig. 1. MSE of interpolators as a function of noise. Axial diffusivity recovery (left), principal eigenvector recovery (middle), radial diffusivity recovery (right). Ground truth tensor: orientation ( √ 2/2, √ 2/2, 0) ; diffusivities 10 -3 (1.71, 0.3, 0.1) mm 2 /s.

Fig. 2 .

 2 Fig. 2. Coronal View of Reconstructed Pyramidal Tracts (with overlaid PT atlas). Columns correspond to increasing number of filtering ROIs: R1 (1st column), R1+R2 (2nd column), R1+R2+R3 (3rd column), R1+R2+R3+R4 (4th column). Rows correspond to the methods (from top to bottom): single-tensor FACT -log-Euclidean, and Bayes, multi-tensor FACT -log-Euclidean, and Bayes.


Konstantin E Avrachenkov
email: k.avrachenkov@inria.fr

Vivek S Borkar
email: borkar.vs@gmail.com

Whittle Index Policy for Crawling Ephemeral Content

We consider the task of scheduling a crawler to retrieve from several sites their ephemeral content. This is content, such as news or posts at social network groups, for which a user typically loses interest after some days or hours. Thus development of a timely crawling policy for ephemeral information sources is very important. We first formulate this problem as an optimal control problem with average reward. The reward can be measured in terms of the number of clicks or relevant search requests. The problem in its exact formulation suffers from the curse of dimensionality and quickly becomes intractable even with moderate number of information sources. Fortunately, this problem admits a Whittle index, a celebrated heuristics which leads to problem decomposition and to a very simple and efficient crawling policy. We derive the Whittle index for a simple deterministic model and provide its theoretical justification. We also outline an extension to a fully stochastic model.

I. INTRODUCTION

Nowadays an overwhelming majority of people find new information on the web at news sites, blogs, forums and social networking groups. Moreover, most information consumed is ephemeral in nature, that is, people tend to lose their interest in the content in several days or hours. The interest in a content can be measured in terms of clicks or number of relevant search requests. It has been demonstrated that the interest decreases exponentially over time [START_REF] Goyal | Learning Influence Probabilities in Social Networks[END_REF], [START_REF] Lefortier | Timely Crawling of High-quality Ephemeral New Content[END_REF], [START_REF] Moon | Refining Recency Search Results with User Click Feedback[END_REF].

In a series of works (see e.g., [START_REF] Cho | Synchronizing a Database to Improve Freshness[END_REF], [START_REF] Cho | Effective Change Detection Using Sampling[END_REF], [START_REF] Cho | Effective Page Refresh Policies for Web Crawlers[END_REF], [START_REF] Olston | Web Crawling[END_REF] and references therein) the authors address the problem of refreshing documents in a database. However, these works do not consider the ephemeral nature of the information. Motivated by this challenge, the authors of [START_REF] Lefortier | Timely Crawling of High-quality Ephemeral New Content[END_REF] suggest a procedure for optimal crawling of ephemeral content. Specifically, they formulate an optimization problem for finding optimal frequencies of crawling for various information sources.

The approach presented in [START_REF] Lefortier | Timely Crawling of High-quality Ephemeral New Content[END_REF] is based on static optimization. In particular, it does not depend on any 'state variable(s)' evolving with time. With a dynamic policy, for instance, if there is not much new material on the principal information sources, the crawler could spend some time to crawl the sources with less popular content but which nevertheless bring noticeable rate of clicks or increase information diversity. In fact, in the numerical results section we demonstrate the superior performance of the Whittle index policy in comparison with some commonly used static policies. Therefore, in the present work we investigate a dynamic formulation of the problem as an optimal control problem with average reward. The direct application of dynamic programming quickly becomes intractable even with moderate number of information sources, due to the so-called curse of dimensionality. Fortunately, the problem admits a Whittle index, which leads to problem decomposition and to a very simple and efficient crawling policy. We derive the Whittle index and provide its theoretical justification first for the deterministic model. We also propose an online algorithm for the relaxed control, which can be used in the case of soft constraints and can provide bounds on the performance of the Whittle index policy. Finally, we study the Whittle index approach in the fully stochastic case when the index does not have an explicit form but can be estimated online.

In [START_REF] Avrachenkov | Optimal Threshold Control by the Robots of Web Search Engines with Obsolescence of Documents[END_REF], [START_REF] Liu | Optimization Issues in Web Search Engines[END_REF], [START_REF] Talim | Controlling the Robots of Web Search Engines[END_REF] the authors study the interaction between the crawler and the indexing engine by means of optimization and control theoretic approaches. One of interesting future research directions is to take into account the indexing engine dynamics in the present context.

The general concept of the Whittle index was introduced by P. Whittle in [START_REF] Whittle | Restless Bandits: Activity Allocation in a Changing World[END_REF]. This has been a very successful heuristic for restless bandits, which, while suboptimal in general, is provably optimal in an asymptotic sense [START_REF] Verloop | Asymptotically Optimal Priority Policies for Indexable and Non-indexable Restless Bandits[END_REF], [START_REF] Weber | On an Index Policy for Restless Bandits[END_REF] and has good empirical performance. It and its variants have been used extensively in logistical and engineering applications, some recent instances of the latter in communications and control being for sensor scheduling [START_REF] Nino-Mora | Sensor Scheduling for Hunting Elusive Hiding Targets via Whittle's Restless Bandit Index Policy[END_REF], multi-UAV coordination [START_REF] Ny | Multi-UAV Dynamic Routing with Partial Observations Using Restless Bandit Allocation Indices[END_REF], congestion control [START_REF] Avrachenkov | Congestion Control of TCP Flows in Internet Routers by Means of Index Policy[END_REF], [START_REF] Avrachenkov | Infinite Horizon Optimal Impulsive Control with Applications to Internet Congestion Control[END_REF], [START_REF] Jacko | Congestion Avoidance with Future-Path Information[END_REF], channel allocation in wireless networks [START_REF] Larranaga | Stochastic and Fluid Index Policies for Resource Allocation Problems[END_REF], cognitive radio [START_REF] Liu | Indexability of Restless Bandit Problems and Optimality of Whittle Index for Dynamic Multichannel Access[END_REF] and real-time wireless multicast [START_REF] Raghunathan | Index Policies for Real-time Multicast Scheduling for Wireless Bradcast Systems[END_REF]. Book length treatments of indexable restless bandits appear in [START_REF] Jacko | Dynamic Priority Allocation in Restless Bandit Models[END_REF], [START_REF] Ruiz-Hernandez | Indexable Restless Bandits[END_REF]. This is a significantly extended version of the conference paper [START_REF] Avrachenkov | Whittle Index Policy for Crawling Ephemeral Content[END_REF]. We have supplied the detailed proofs omitted in ibid. and filled in a gap in the Whittle idexability argument. We also take this opportunity to correct the flawed treatment of the stochastic case in ibid. (see new section VIII below). That section has a large part of new material on the Whitttle index in the fully stochastic case, in particular, on its online estimation. We have also added a completely new section on the relaxed control problem (section VII). We have also significantly reworked the numerical results section adding comparison with the other heuristic policies.

II. MODEL

There are N sources of ephemeral content. A content at source i ∈ {1, ..., N } is published with an initial utility modelled by a nonnegative random variable ξ i and decreasing exponentially over time with a deterministic rate µ i . This model of exponential evolution of the content value was introduced and justified by measurements in [START_REF] Lefortier | Timely Crawling of High-quality Ephemeral New Content[END_REF]. In practical terms, the value of content can be measured in the number of clicks or in the number of relevant queries. Clearly, this sort of information is available to search engines. We further assume that the new content arrives at source i ∈ {1, ..., N } according to a time-homogeneous Poisson process with rate Λ i . Thus, if source i's content is crawled τ time units after its creation, its utility is given by ξ i exp(-µ i τ ). The base utility ξ i is assumed independent identically distributed across contents at a given source, with a finite mean ξi . It is also assumed independent across sources. We assume that the crawler crawls periodically at multiples of time T > 0 and has to choose at each such instant which sources to crawl, subject to a constraint we shall soon specify. When the crawler crawls a content source, we assume that the crawling is done in an exhaustive manner. In such a case, the crawler obtains the following expected reward from crawling the content of source i:

u i = Λ i E[ξ i exp(-µ i τ )] = Λ i ξi µ i (1 -exp(-µ i T )). (1)
Set α i = exp(-µ i T ). Let us define the state of source i at time t as the total expected utility of its content, denoted by X i (t). Then, if we do not crawl source i at epoch t (formally, the control is v i (t) = 0 -we say the source is 'passive'), we obtain zero reward r i (X i (t))v i (t) = 0 and the state evolves as follows:

X i (t + 1) = α i X i (t) + u i . (2)
On the other hand, if we crawl source i (formally, v i (t) = 1we say the source is 'active'), we obtain the expected reward r i (X i (t))v i (t) = X i (t) and the next state of the source is given by

X i (t + 1) = u i . (3)
This is de facto a deterministic model, even though it captures via equation (1) some effects of stochasticity of the underlying Poisson process. An analogous stochastic model which fully takes into account the randomness coming from the underlying arrival process and varying content value will be discussed in detail in Section VIII.

Our aim is to maximize the long run average reward

lim sup t↑∞ N i=1 1 t t m=0 r(X i (m))v i (m) (4)
subject to the constraint

lim sup t↑∞ N i=1 1 t t m=0 C i v i (m) = M (5)
for a prescribed M > 0. If C i = 1, i = 1, ..., N , this case can be interpreted as a constraint on the average number of crawled sites per crawling period T and corresponds to the original Whittle framework for restless bandits [START_REF] Whittle | Restless Bandits: Activity Allocation in a Changing World[END_REF], where this framework arises as a relaxation of a hard constraint of 'at most M ' crawled sites in each time slot (we discuss this later in detail.) The case C i = 1 is slightly more general and can represent the situation when various sites have different limits on the crawling rates (typically specified in the file 'robots.txt').

This is a constrained average reward control problem [START_REF] Altman | Constrained Markov Decision Processes[END_REF], [START_REF] Piunovskiy | Optimal Control of Random Sequences in Problems with Constraints[END_REF]. We address this problem in the framework of restless bandits and derive a simple index policy for the problem, which may be viewed as a variant of the celebrated Whittle index. In the next section, we recall the theory of Whittle index.

III. WHITTLE INDEX

The original formulation of restless bandits is for discrete state space Markov chains, but we consider here Markov chains with closed domains (i.e., closure of an open set) S i ⊂ R d , d ≥ 1, as state space. The original motivation for the index policy remains valid nevertheless as long as we justify the associated dynamic programming equation, which we do. A deterministic dynamics such as ours is a special case, albeit degenerate. The fully stochastic case can be handled similarly and is detailed in a separate Section VIII. While in this section we introduce the broader framework in a general set up, we use the same notation as above to highlight the correspondences. This should not cause any confusion. Thus consider resp. S i -valued processes X i (t), t ≥ 0, 1 ≤ i ≤ N , each with two possible dynamics, dubbed active and passive, wherein they are governed by transition kernels p i (dy|x), q i (dy|x), resp. These are assumed to be continous as maps x ∈ S i → P(S i ). (:= the space of probability measures on S i with Prohorov topology). The control at time t is an

A := {0, 1} N -valued vector v(t) = [v 1 (t), • • • , v N (t)] ∈ A, with the understanding that v i (t) = 1 ⇐⇒ X i (t) is active.
In the original restless bandit problem, exactly N < N processes are active at any given time. The v i (t) are assumed to be adapted to the history, i.e., the σ-field σ(X i (s), s ≤ t; v i (s), s < t; 1 ≤ i ≤ N ). Let r i : S → R + , 1 ≤ i ≤ N, be reward functions so that a reward of r i (X i (t)) is accrued if process i is active at time t. The objective then is to maximize the long run average reward

β := lim sup t↑∞ N i=1 1 t t s=0 E[r i (X i (s))v i (s)].
This problem has state space × N i=1 S i . Whittle's heuristic among other things reduces the problem to separate control problems on S i . The idea is to relax the constraint of 'exactly N are active' to 'on average, N are active', i.e., to

lim sup t↑∞ 1 t t s=0 E[ N i=1 v i (s)] = N .
This makes it a constrained average reward control problem [START_REF] Altman | Constrained Markov Decision Processes[END_REF], [START_REF] Piunovskiy | Optimal Control of Random Sequences in Problems with Constraints[END_REF] which permits a relaxation to an unconstrained average reward problem by replacing the above reward by

lim sup t↑∞ N i=1 1 t t s=0 E[r i (X i (s))v i (s) + λ(N /N -v i (s))],
where λ ∈ R is the Lagrange multiplier. Motivated by this, Whittle introduced a 'subsidy' λ for passivity, i.e., a virtual reward for a process in passive mode. Replace the above control problem by N control problems with the ith problem for process X i (•) seeking to maximize over admissible v i (t), t ≥ 0, the reward

lim sup t↑∞ 1 t t s=0 E[r i (X i (s))v i (s) + λ(1 -v i (s))]. (6)
The dynamic programming equation for this average reward problem is

V i (x) + β = max λ + q i (dy|x)V i (y), r i (x) + p i (dy|x)V i (y) , (7)
where V i (x) is the value (or bias) function. If this can be rigorously justified (which is not always easy), one defines B(λ) as the set of passive states, i.e., B(λ) :=

x : λ + q i (dy|x)V i (y) ≥ r i (x) + p i (dy|x)V i (y) .
If B(λ) increases monotonically from φ to S i as λ increases from -∞ to ∞, the problem is said to be Whittle indexable.

The Whittle index for the ith process in state x i is then defined as the value of λ for which the active and passive modes are equally desirable, i.e.,

γ i (x i ) := {λ : λ + q i (dy|x i )V (y) = r i (x i ) + p i (dy|x i )V (y)}.
The Lagrange multiplier multiplier λ can be interpreted as a subsidy for passivity. Then, the Whittle index has an interpretation of such subsidy that makes active and passive actions equally attractive at the current state. The so-called 'Whittle index policy' [START_REF] Whittle | Restless Bandits: Activity Allocation in a Changing World[END_REF] then is to set v i (t) = 1 for the i with the top N indices and v j (t) = 0 for the rest.

IV. DYNAMIC PROGRAMMING EQUATION

In view of the above, the first step is to justify the counterpart of (7) in our context. For this, we first note that r

i (x) = x, 1 ≤ i ≤ N . Further, let u * i := ui 1-αi > u i .
We argue that without loss of generality, we may take

S i = [u i , u * i ]. To see this, let X i (0) = x 0 . If x 0 ≤ u i and is never crawled, it is easy to see that X i (t) = α t i x 0 + (1 -α t i )u * i ↑ u * i . On the other hand, if x 0 > u * i , then X i (t) = α t i x 0 + (1 -α t i )u * i ↓ u * i as t ↑ ∞ if never crawled.
(Throughout the paper we use ↑ and ↓ to denote the limits from below and above, respectively.) It is reset to u i in either case if there is even a single crawl. Combining these observations and recalling that we consider the long-run average criterion, we conclude that any

x 0 / ∈ [u i , u * i ]
is a transient state and can be ignored. Thus we set

S i = [u i , u * i ].
Henceforth we focus on the average reward problem for source i. We do not delve into the justification for Lagrange multiplier formulation for constrained average reward problem on a general state space, as this is well understood [START_REF] Altman | Constrained Markov Decision Processes[END_REF], [START_REF] Borkar | Convex Analytic Methods in Markov Decision Processes[END_REF], [START_REF] Piunovskiy | Optimal Control of Random Sequences in Problems with Constraints[END_REF]. (In fact, it follows from standard Lagrange multiplier theory applied to the 'occupation measure' formulation of average reward problem which casts it as an abstract linear program. See section 4.2 of [START_REF] Borkar | Convex Analytic Methods in Markov Decision Processes[END_REF] which carries out this program for discrete state space and section 3.2 of ibid. which describes how to extend the same to general compact Polish state spaces as long as the controlled transition probability kernel is continuous in the initial state and control.) For notational simplicity we drop the index i for the time being. We approach the problem by the standard 'vanishing discount' argument. Thus let 0 < δ < 1 be a discount factor and for k(x, v)

:= xv + Cλ(1 -v), consider the infinite horizon discounted reward ∞ m=0 δ t k(X(t), v(t)).
(Here C > 0 replaces the C i in (5).) Denote the associated value function by

V δ (λ, x) := sup {v(t)},X(0)=x ∞ m=0 δ t k(X(t), v(t)) .
We drop the argument λ when it is understood from the context. Then V δ satisfies the discounted reward dynamic programming equation

V δ (λ, x) = max (Cλ + δV δ (λ, αx + u), x + δV δ (λ, u)) . (8) Lemma 1
The solution of equation ( 8) has the following properties:

(1) Equation ( 8) has a unique bounded continuous solution V δ ;

(2) V δ is Lipschitz uniformly in δ ∈ (0, 1);

(3) V δ is monotone increasing and convex. Proof: Claim (1) is standard (See Theorem 4.2.3 and bullet 1 in 'Notes on §4.2', Section 4.2, [START_REF] Hernández-Lerma | Discrete Time Markov Control Processes: Basic Optimality Criteria[END_REF]). For (2), take x, x ∈ S with x > x. Consider processes X(t), t ≥ 0, and X (t), t ≥ 0, with initial conditions x, x resp., both controlled by control sequence v(t), t ≥ 0, that is optimal for the former. Then

V δ (x ) -V δ (x) ≤ ∞ t=0 δ m (k(X (t), v(t)) -k(X(t), v(t))) = 1 -(αδ) τ 1 -α (x -x),
where τ := the time of first crawl (= ∞ if never crawled).

Interchanging the roles of x , x we get a symmetric inequality, whence it follows that

|V δ (x ) -V δ (x)| ≤ 1 -(αδ) τ 1 -α |x -x|.
For the first part of (3), take x > x as above and let X (t), X(t), t ≥ 0, be processes generated by a common admissible control sequence {v(t)} with initial conditions x , x resp. Then it is easy to check that X (t) ≥ X(t) for all t. Therefore

∞ t=0 δ t k(X (t), v(t)) ≥ ∞ t=0 δ t k(X (t), v(t)). (9)
Taking supremum over all admissible controls on both sides, monotonicity of V δ follows. For convexity, define the finite horizon discounted value function

V n (x) = sup {v(t)},X(0)=x n t=0 δ t k(X(t), v(t)).
Then it satisfies the dynamic programming equation

V n (x) = max (Cλ + δV n-1 (αx + u), x + δV n-1 (u))
for n ≥ 1 with V 0 (x) = x. The convexity of V n for each n then follows by a simple induction. Since

V δ (x) = lim n↑∞ V n (x), V δ is also convex. 2 Define Vδ (x) = V δ (x) -V δ (u), x ∈ S.
Then by the above lemma, Vδ is bounded Lipschitz, monotone and convex with Vδ (u) = 0. Also, (1 -δ)V δ (u) is bounded. Using Arzela-Ascoli and Bolzano-Weirstrass theorems, we may pick a subsequence such that ( Vδ , (1-δ)V δ (u)) converge in C(S)×R to (say) (V, β). From (8), we have

Vδ (x) + (1 -δ)V δ (u) = max Cλ + δ Vδ (αx + u), x .
Passing to the limit along an appropriate subsequence as δ ↑ 1, we have

V (x) + β = max (Cλ + V (αx + u), x) (10)
= max

v∈{0,1} vx + (1 -v)(Cλ + V (αx + u)) . (11)
Then [START_REF] Cho | Effective Change Detection Using Sampling[END_REF] is the desired dynamic programming equation for average reward. We study important structural properties of the value function V in the next section.

V. PROPERTIES OF THE VALUE FUNCTION

We begin with the following result.

Lemma 2

The following statements hold:

(1) V is monotone increasing, convex with V (u) = 0;

(2) The maximizer on the right hand side of ( 11) is the optimal control choice at state x and β is the optimal reward. Proof: Since monotonicity and convexity are preserved in pointwise limits, the first claim is immediate.

For the second point, let v * (x) denote the maximizer on the r.h.s. of [START_REF] Goyal | Learning Influence Probabilities in Social Networks[END_REF], any tie being settled arbitrarily. Then under

{v(t) = v * (X(t)), t ≥ 0}, V (X(t)) + β = k(X(t), v(t)) + V (X(t + 1)). (12)
Summing (12) over t = 1, 2, • • • , T , and dividing by T on both sides, then letting T ↑ ∞, we see that β is equal to the average reward under this control policy. On the other hand, for any other control sequence, the equality in (12) will be replaced by ≥, leading to the conclusion that β is greater than or equal to the corresponding average reward by an argument similar to the above. This implies the second claim. 2 Now define

B := {x ∈ S : Cλ + V (αx + u) > x}, B c := {x ∈ S : Cλ + V (αx + u) ≤ x}.
These are respectively the sets of passive and active states under subsidy λ.

Recall the stopping time τ := the time of first crawl. Suppose τ < ∞. (The case τ = ∞ corresponds to 'never crawl' which we consider separately below.) Under optimal policy, iterating equation [START_REF] Cho | Effective Change Detection Using Sampling[END_REF] τ times leads to

V (x) = (Cλ -β)τ + α τ x + 1 -α τ 1 -α u -β .
Under any other policy, we would likewise obtain

V (x) ≥ (Cλ -β)τ + α τ x + 1 -α τ 1 -α u -β .
Thus we have the explicit representation for V given by

V (x) = max (Cλ -β)τ + α τ x + 1 -α τ 1 -α u -β ,
where the maximum is over all admissible control sequences.

In particular, this implies:

Lemma 3 Equation ( 10) has a unique solution.

We shall now eliminate some irrelevant situations. 1) If u * ∈ B, i.e., the optimal action at u * is 0, then u * is a fixed point of the optimally controlled dynamics and the corresponding reward is Cλ. Then β = Cλ and it is optimal to be passive at all states, i.e., B = [u, u * ], B c = φ, and

λ ≥ λ M := max x∈[u,u * ] (x -V (αx + u))/C. (13)
2) If u ∈ B c , then from [START_REF] Cho | Effective Change Detection Using Sampling[END_REF], 0 + β = u, i.e., β = u and it is optimal to crawl when at u. Then u is a fixed point of the controlled dynamics and it is optimal to be active at all states, i.e., B c = [u, u * ], B = φ, and

λ ≤ λ m := min x∈[u,u * ] (x -V (αx + u))/C. (14)
Note that since constant policies v(t) ≡ 0 and v(t) ≡ 1 lead to costs Cλ and u resp., β ≥ (Cλ) ∨ u always and β > (Cλ) ∨ u for λ ∈ (λ m , λ M ). For each λ in (λ m , λ M ), both B, B c are non-empty and there exists an a ∈ (u, u * ) for which the choice of being active or passive is equally desirable. However, since at u * the optimal action is to crawl, we conclude that u * ∈ B c and only the second possibility can occur.

2

Corollary 1 The map x → x -V (αx + u) is monotone non-decreasing on [u, u * ].
Thus the optimal policy is a threshold policy, viz., it is passive for states below a certain threshold and allows the dynamics x k+1 = αx k + u to evolve uninterrupted till the threshold is crossed, then gets active and resets the state to u. In particular, it is periodic after the first reset and x k takes only finitely many distinct values thereafter.

Lemma 5

The above problem is Whittle indexable. Proof: We shall write the value function as V (λ, x) to render its λ-dependence explicit. It suffices to show that the optimal threshold x * (λ) increases with λ for λ ∈ (λ m , λ M ). Since the optimal policy is threshold, we may restrict attention to threshold policies with a threshold in (u, u * ). In particular, this implies a periodic trajectory with a finite range after the first reset, therefore the 'occupation measure' ν on [u, u * ] given by

ν(A) := lim n↑∞ n k=0 I{x k ∈ A} n for A Borel in [u, u * ]
is well defined as a probability measure on [u, u * ]. Then the corresponding average reward is

Cλν(B) + B c xν(dx), (15)
an affine function of λ. Moreover, its slope is ν(B) < some ∆ < 1, otherwise the optimal reward β(λ) would be equal to Cλ, a situation we have excluded from consideration. Since β(λ) is the supremum of ( 15) over all threshold policies, it is convex increasing and Lipschitz, hence absolutely continuous, with β (λ) := d dλ β(λ) ≤ C∆ a.e. For x ∈ (x * (λ), u * ], x is active. Then for x := αx + u,

x ≤ u * = u 1 -α =⇒ x = αx + u ≥ αx + (1 -α)x =⇒ x ≥ x =⇒ x ∈ B c (by Lemma 4) =⇒ V (x ) = x -β.
Combining this with (10), we have the following expression for the threshold:

x * (λ) = Cλ + u -β(λ) 1 -α . ( 16
)
Since C > C∆ ≥ β (λ), this is increasing in λ. Thus the passive set is monotonically increasing, establishing Whittle indexability. 2 VI. DERIVATION OF WHITTLE INDEX Consider the situation when λ = γ(x) for a prescribed x ∈ (u, u * ). It is clear that after the first crawl when the process is reset to u, the optimal X(t) becomes periodic: not crawling and increasing till it hits B c and then crawling -thereby being reset to u -to repeat the process. Since finite initial patches do not affect the long run average reward, we may then take

X(0) = u. Define η(x) = min{t : X(t) ∈ B c }. Then X(η(x)) = (1 -α η(x) )u * (17) =⇒ η(x) = log + α 1 - x u * , ( 18
)
where log + α x = (log α x)I{x > 0}. Since by the renewal reward argument the long run average reward is equal to the average reward over one period, we can write

β = Cλ(η(x) -1) + X(η(x)) η(x) , (19)
where η(x) is given by ( 18) and X(η(x)) is given by [START_REF] Liu | Optimization Issues in Web Search Engines[END_REF].

We now revert to using the index i to identify the source being referred to. In particular, β i , λ i will refer to the optimal reward, resp. Lagrange multiplier, for the ith decoupled problem. Our main result is: Theorem 1 The Whittle index for our problem is given by

γ i (x) := 1 C i η i (x)((1 -α i )x -u i ) + 1 -α ηi(x) i 1 -α i u i , where η i (x) := log + αi u i -(1 -α i )x u i .
Therefore the index policy is to crawl at time t (= mT for some m ≥ 0) the top M sources according to decreasing values of γ i (X i (t)), or alternatively, choose a number of top sources for the constraint to be reached.

Remark: Note that if an arm (say, ith) is crawled even once, the corresponding state process {X i (t)} takes only discrete values thereafter. These depend on α i and u i alone. In fact this is also true for an arm that is never crawled, except that the discrete values taken will also depend on the initial condition. Therefore we need restrict attention to only these values of x for the argument of γ i (•). This results in a further simplification of the index formula, to

γ i (x) = 1 C i (η i (x)((1 -α i )x -u i ) + x) ,
where η i (x) is as before, but without the and the argument x of both γ i and η i is now restricted to the aforementioned discrete set.

Proof: We drop the subscript i for notational convenience. For x ∈ B c , (10) leads to V (x) = x -β. Recall from the proof of Lemma 5 that

x := αx + u =⇒ V (x ) = x -β.
Combining this with [START_REF] Cho | Effective Change Detection Using Sampling[END_REF] and the definition of Whittle index implies that for our problem, it is

γ i (x) = (1 -α i )x -u i + βi (x) C i , ( 20
)
where by virtue of [START_REF] Milgrom | Envelope Theorems for Arbitrary Choice Sets[END_REF], βi (x) := the optimal cost if one were to set λ i = γ i (x). The latter is given by:

βi (x) := 1 η i (x) C i γ i (x)(η i (x) -1) + 1 -α ηi(x) i u * i ,
where

η i (x) := log + αi u i -(1 -α i )x u i .
Substituting this back into [START_REF] Moon | Refining Recency Search Results with User Click Feedback[END_REF], one gets a linear equation for γ i (x) that can be solved to evaluate γ i (x) as

γ i (x) := 1 C i η i (x)((1 -α i )x -u i ) + 1 -α ηi(x) i 1 -α i u i .
This completes the proof. 2

VII. ONLINE ALGORITHM FOR THE RELAXED CONTROL

PROBLEM

If one is satisfied with the time average constraint on the number of crawled sources rather than a hard integer constraint, we can propose an online algorithm based on dual descent. This relaxation of the control problem is also useful to establish an upper bound on the optimal value, therefore also on the Whittle index policy.

In order to apply the dual descent, we first need to derive an expression for the optimal threshold for a stand alone information source. Let us substitute ( 17) into [START_REF] Milgrom | Envelope Theorems for Arbitrary Choice Sets[END_REF] and differentiate the resulting expression for β with respect to η. Then, equating the derivative to zero, we obtain the following equation with respect to η

α η (1 -log(α)η) = 1 -Cλ/u * .
The solution of this equation can be expressed in terms of the Lambert W function (W (z) exp(W (z)) = z) as follows:

η = 1 log(α) 1 + W (e -1 (Cλ/u * -1)) .
The optimal threshold value is then given by

x(λ) = u * (1 -exp(1 + W (e -1 (Cλ/u * -1)))). (21)
This allows us to learn the optimal Lagrange multiplier by dual stochastic descent given by

λ(k + 1) = λ(n) -b(k)(M/N - i C i v i (k)), k ≥ 0, (22)
where stepsizes {b(k)} satisfy the usual conditions for stochastic approximation stepsize:

n b(n) = ∞, n b(n) 2 < ∞.
The fact that ( 22) is indeed a stochastic gradient descent for the Lagrange multiplier follows from the generalized envelope theorem [START_REF] Milgrom | Envelope Theorems for Arbitrary Choice Sets[END_REF] as in [START_REF] Salodkar | A Stable Online Algorithm for energy-efficient multi-user scheduling[END_REF], Lemma 2, for which it is essential that at each time k, the v i (k)'s be chosen according to the optimal threshold policy corresponding to λ = λ(k) as derived above. Being a stochastic gradient descent for a convex function, its convergence to the minimizer, i.e., the actual Lagrange multiplier, is easy to establish by standard arguments (see, e.g., [START_REF] Borkar | Stochastic Approximation: A Dynamical Systems Viewpoint[END_REF], chapter 10). Note that the optimal scheme for the relaxed problem is precisely to crawl those sites for which the Whittle index is above the Lagrange multiplier. In online version above, one simply uses the estimate λ(k) at time k in lieu of the true Lagrange multiplier.

VIII. STOCHASTIC CASE

We now consider the fully stochastic situation when the content generation at each source is observed as a random process. In fact, one could also consider mixed situations when some sources are observed and others are not. As we shall see, the development closely mimics the foregoing.

The stochastic system dynamics can be described as follows: Let {τ i n } denote the successive arrival times of content at source i, with utilities {ξ i n }, resp. We assume these arrival processes to be Poisson with rate possibly depending on the site i, and {ξ i n } to be independent and independent of the process of arrival epochs. Further, we assume that they are identically distributed for each fixed i. The net utility added to source i during k-th epoch will be

U i (k) := τ i n : (k-1)T ≤τ i n <kT ξ i n e -µi(kT -τ i n ) .
Under our hypotheses, the system state at time (k + 1)T is then

X i (k + 1) = α i X i (k) + U i (k + 1) if no crawl, U i (k + 1) if crawled. ( 23
)
We define the average reward as

lim sup t↑∞ N i=1 1 t t m=0 E[r(X i (t))v i (t)],
which we seek to maximize subject to the constraint

lim sup t↑∞ 1 t N i=1 C i E[v i (t)] = M.
Henceforth we drop the subscript i for notational economy, till such time as it is requited. The discounted value function

V δ (x) := sup {v(t)},X(0)=x E ∞ t=0 δ t k(X(t), v(t))
then satisfies the dynamic programming equation

V δ (x) = (24) max Cλ + δ V δ (αx + u)ϕ(du), x + δ V δ (u)ϕ(du) ,
where ϕ is the law of U (t), ∀t.

Lemma 6

The conclusions of Lemma 1 continue to hold.

Proof: The first claim follows as before from the cited results of [START_REF] Hernández-Lerma | Discrete Time Markov Control Processes: Basic Optimality Criteria[END_REF]. For the second, let X(t), X (t) be defined with identical arrival and control processes with the control process v(•) being optimal for X(•), differing only in their initial conditions resp., x, x . Then

V δ (x ) -V δ (x) ≤ E[ ∞ t=0 δ m (k(X (t), v(t)) -k(X(t), v(t)))] ≤ E 1 -(αδ) τ 1 -αδ (x -x),
where τ := min{t ≥ 0 : X(t) = X (t)} is the coupling time, which is also the time of the first crawl. The Lipschitz property follows as before. Next let X(t), X (t) be as in the proof of Lemma 1 (3). Taking expectations in (9) followed by a supremum over all admissible controls proves monotonicity. Convexity also follows as in the deterministic case. 2 The 'vanishing discount' argument of Section IV can now be used to establish the average cost dynamic programming equation

V (x) + β = max(Cλ + V (αx + u)ϕ(du), x), (25)
where we have rendered V unique by setting V (u)ϕ(du) = 0. Monotonicity and convexity properties of V follow as in Lemma 2, using the additional fact that convexity of a function F also implies that of x → F (ax + y)ν(dy) for any scalar a and probability measure ν on R. Equation ( 12) gets modified to

E[V (X(t))] + β = E[k(X(t), v(t))] + E[V (X(t + 1))],
from which the optimality of

v * (x) ∈ Argmax v vx + (1 -v)(λ + V (αx + u)ϕ(du)) ,
with x ∈ S, follows by arguments analogous to those of Lemma 2. Furthermore, iterating the above, V can be shown to be the unique solution of ( 25) by establishing the explicit representation

V (x) = max E x (Cλ -β(λ))σ + α σ x + σ t=1 α σ-t U (t) -β(λ) = max E x [(Cλ -β(λ))σ + α σ x + (1 -α σ )u * ] -β(λ) (26)
for passive x and

V (x) = x -β(λ) (27)
for active x. Here:

• the maximum is over all admissible control sequences {U (t)}, • σ is the first time to cross the threshold, i.e., first time to crawl,

• u := E[U (m)], u * := u
1-α , and the second equality in [START_REF] Raghunathan | Index Policies for Real-time Multicast Scheduling for Wireless Bradcast Systems[END_REF] follows from the optional sampling theorem. As before, we assume that β(λ) > Cλ ∧ u so as to eliminate uninteresting cases. Then the right hand side of (26) will be < x -β(λ) for x > u * , implying in particular that it is less than the right hand side of [START_REF] Ruiz-Hernandez | Indexable Restless Bandits[END_REF], that is, x could not have been passive. Then by the same arguments as for the deterministic case, we can restrict the state space to [u, u * ].

The definitions of B, B c now change to

B := {x ∈ S : Cλ + V (αx + u)ϕ(du) > x}, B c := {x ∈ S : Cλ + V (αx + u)ϕ(du) ≤ x}.
Making similar arguments as those in Lemma 4, we establish that the optimal control policy is of threshold type.

Lemma 7

The problem is Whittle indexable. Proof: Given that the maximum in ( 26) is over all admissible control policies and that the optimal policy is a threshold policy, it is also the maximum over all threshold policies. Fix a threshold, an initial condition x ∈ [u, u * ] and consider a process corresponding to this threshold policy. Then the random variable σ is independent of λ. Since β (λ) < C, it then follows that the expectation on the r.h.s. of ( 26) is monotone increasing in λ. Then the same is true for its maximum over all threshold policies. Letting x * (λ) denote the optimal threshold under λ, by definition both ( 26), ( 27)

hold for x = x * (λ). Let G(λ, x) := max E x [(Cλ -β(λ))σ + α σ x + (1 -α σ )u * ] , ∀x ∈ B, x, ∀x ∈ B c ,
where the maximum is over all threshold policies. Since V is convex increasing, so is G(λ, •). Also, G is increasing in λ. Then x * (λ) is a fixed point of G(λ, •). Since it is optimal to be passive at x = u (because β(λ) > u), G(λ, u) > u ∀λ. On the other hand, it is optimal to be active at u * and therefore G(λ, u * ) = u * . Thus the convex curve x → G(λ, x) crosses the line y = x exactly once in [u, u * ) and the point of crossing, which by definition is x * (λ), increases with λ (because G does). This proves Whittle indexability. 2 Let η denote the first hitting time ≥ 1 of B c when X(0) = U (0), i.e., the time to the next crawl after 0. Then

X(η) = η t=1 α η-t U (t), m ≥ 1.
From dynamic programming arguments, we know that it suffices to consider stationary Markov controls. Consider λ = γ(x), the Whittle index of x. Let η(x) := E[η] and β(x) := the optimal cost, both corresponding to λ = γ(x). Then by standard renewal-reward arguments,

β(x) = Cλ(η(x) -1) + E[X(η(x))] η(x) . Define Z 0 = 0, Z n := n m=1 α n-m U (m), n ≥ 1, (e.g., Z 2 = U (2) + αU (1)),
From the definition of γ(x), it follows that x ∈ B c , so the optimal decision at x is to crawl. By definition of γ(x),

x = Cγ(x) + E[V (αx + U (1))].
sources per crawling period, i.e., M = 2. It is interesting that now the policy becomes much less regular. Source 1 is always crawled. Sources 2 and 3 are crawled in a non-trivial periodic way and sources 4 is crawled periodically with a rather long period.

Let us now consider the stochastic model. Since the computation of the Whittle index for the stochastic model is very cumbersome, even numerically, we propose to use deterministic Whittle index for the stochastic model as well. We compare the performance of the deterministic Whittle index policy with the Round Robin (RR) policy and the Greedy-type policy, which chooses the source maximizing

Λ i ξi µ i (1 -exp(-µ i × t i )) ,
where t i is the time elapsed since the last crawl of source i. To have good averaging, we increased the time interval of simulations to 10000 steps. It is very interesting to observe that if we have randomness only in the content arrival process (we model the content arrival by the Poisson process), the Whittle index policy produces the average reward 253.1, which is smaller than the reward produced by the Greedy-type policy, 260.2. The RR policy results in just 208.4. However, if we introduce randomness also in the initial value of the content (e.g., we now model it as exponential random variable with the parameter ξ i given in Table I), the Whittle index policy results in the average reward 283.8 and outperforms both the Greedy-type policy, 259.3, and RR policy, 207.3. Thus, we can recommend using the deterministic Whittle policy heuristic when we have high level of randomness in the system. As example, in Figure 4 we present the state evolution of the stochastic model with dynamics [START_REF] Olston | Web Crawling[END_REF], with Poisson arrivals but deterministic value of the content using the deterministic Whittle index policy. One can see that now in the stochastic setting source 1 is crawled from time to time.

X. CONCLUSIONS AND FUTURE WORK

We have formulated the problem of crawling web sites with ephemeral content as an average reward optimal control problem and have shown that it is indexable. We have found that the Whittle index has a very simple form for a deterministic model and in principle computable for a fully stochastic case. In particular, we have proposed online algorithms for the Whittle index in the fully stochastic case as well as for the relaxed deterministic control. The numerical example demonstrates that the Whittle index policies, unlike the policies suggested in [START_REF] Lefortier | Timely Crawling of High-quality Ephemeral New Content[END_REF], do not generally have a trivial periodic structure. We have compared the Whittle index policy with the Round Robin and Greedy-type policies and have demonstrated that the deterministic Whittle index provides good performance even in the stochastic case with high randomness. The proposed approach can also be used in the cases when some states are observed. In such cases, the Whittle index will act as a self-tuning mechanism. In the future research we plan to elaborate an adaptive version when some parameters (e.g., the rate of new information arrival) need to be estimated online. The other interesting future research directions are addition to the model the dynamics of the indexing engine and practical implementation and testing of the Whittle index approach in real crawling system.

Lemma 4

 4 The sets B, B c are of the form [u, a), [a, u * ] for some a ∈ [u, u * ]. Proof: Since V is convex, one of the following two must hold: 1) For some a 2 > a 1 , B = [u, a 1 ) ∩ (a 2 , u * ] and B c = [a 1 , a 2 ], or, 2) for some a, B = [u, a), B c = [a, u * ].

Fig. 1 .

 1 Fig. 1. Content value as a function of time.

Fig. 2 .

 2 Fig. 2. The case of M = 1.

Fig. 3 .

 3 Fig. 3. The case of M = 2.

Fig. 4 .

 4 Fig. 4. The case of M = 2 (stochastic model).

ACKNOWLEDGMENTS

The authors gratefully acknowledge the discussions with Liudmila O. Prokhorenkova and Egor Samosvat from Yandex during the preparation of the manuscript.

Learning Schemes for Network Analytics" and IFC/DST-Inria-2016-01/448 "Machine Learning for Network Analytics".

Hence,

which can be solved for γ(x) because all three expectations on the right can be calculated from the problem specifications.

One possible computational approach is as follows. If we fix the threshold to be x, consider V satisfying

The Whittle index γ(x) should satisfy γ(x) = λ with

Write V (x), which is the unique solution to ( 28)-( 30), as V (λ, x) to make its λ-dependence explicit. We can then 'learn' γ(x) by the stochastic approximation scheme

which incrementally adjusts γ n in the direction of forcing [START_REF] Weber | On an Index Policy for Restless Bandits[END_REF] to hold. Here { Ũn } are i.i.d. random variables with distribution ϕ. By standard stochastic approximation theory ([7], Chapter 2), this will converge to a solution of (31) as desired. Note that the infinite linear system (28)- [START_REF] Verloop | Asymptotically Optimal Priority Policies for Indexable and Non-indexable Restless Bandits[END_REF] has to be solved at each step. This can be done approximately by solving a discretized version of it by a standard linear system solver. Alternatively, one can concurrently solve this discretized version iteratively but on a faster time scale, using the 'two time scale stochastic approximation' framework of [START_REF] Borkar | Stochastic Approximation: A Dynamical Systems Viewpoint[END_REF], Chapter 6. For example, one can replace [START_REF] Whittle | Restless Bandits: Activity Allocation in a Changing World[END_REF] with

and augment it with

This is nothing but the classical relative value iteration algorithm [START_REF] Puterman | Markov Decision Processes[END_REF] but for the time-dependent γ n . The algorithm is non-incremental, i.e., on the 'natural' time scale dictated by the clock n = 0, 1, 2, • • • , which is faster than the time scale of (33) dictated by the decreasing step sizes {a(n)}, i.e., 'algorithmic clock'

) See ( [START_REF] Borkar | Stochastic Approximation: A Dynamical Systems Viewpoint[END_REF], chapter 6) for detailed analysis of two time scale algorithms. We summarize only the main conclusions here. Because {γ n } are updated on a slower time scale, (34)-(35) see γ n as quasi-static, i.e., nearly a constant, whence the pair reduces to exactly the classical relative value iteration that would converge to the solution V of ( 28)-( 29) corresponding to V (u) = β for λ 'frozen' at γ n . What this means is that if V n denotes the solution to ( 28)-( 29)-( 30), then

a.s. Thus asymptotically, (33) mimics

which converges to the Whittle index a.s. Note that this has to be done separately for each x. In practice one may do so for a finite family of x's and interpolate.

IX. NUMERICAL EXAMPLE

Let us illustrate the theoretical results by a small numerical example. The goal of the numerical example is to obtain further insights into the performance of the Whittle index policy rather than to investigate an implementation of the Whittle index policy in a real crawling system. Real implementation is of course a good future direction for research and development. Let us consider four information sources with parameters given in Table I. Without loss of generality, we take the crawling period T = 1. One can see how the user interest decreases over time for each source in Figure 1. The initial interest in the content of sources 1 and 2 is high, whereas the initial interest in the content of sources 3 and 4 is relatively small. The interest in the content of sources 1 and 3 decreases faster than the interest in the content of sources 2 and 4.

Let us first consider the deterministic model. We run the system with parameter specified in Table I for 1000 time steps. In Figure 2 we show an example of the state evolution under the Whittle index policy and with the constraint that the crawler can visit only one site per crawling period T , i.e., M = 1. The application of Whittle index results in periodic crawling of sources 1 and 2, crawling each with period two. Sources 3 and 4 should never be crawled. This resulted in the average reward 260.30. We can also apply the relaxed control described in Section VII. This gives the average reward 260.96 and, since the relaxed control provides bounds on the deviation from the optimum, this indicates that the Whittle index policy in this particular case gives solution which is very close to the optimal one. While using the relaxed control, the constraint violation appeared in less than 1% of time steps. In contrast, if one always greedily crawls only source 1, he obtains the average reward 179.8.

In Figure 3 we show the state evolution of the bandits under the constraint that the crawler can visit two information