# Indeterminate Stieltjes moment problems revisited

Christian Berg

(11.7.2024)

###### Abstract

We consider a normalized indeterminate Hamburger moment sequence $s$ which is
supposed to be Stieltjes. We revisit old results about
determinacy/indeterminacy in the sense of Stieltjes for $s$ and we prove some
new results about the concepts involved.

Dedicated to Mourad Ismail at the occasion of his 80’th birthday

Mathematics Subject Classification: Primary 44A60; Secondary 30B50

Keywords. Stieltjes and Hamburger moment problems. Indeterminacy. Friedrichs
solution to a Stieltjes problem.

##  1 Introduction

In a recent paper Sokal and Walrad [13] gave a continued fraction
characterization of Stieltjes moment sequences for which there exists a
representing measure with support in $[\xi,\infty)$. In Remark 3 of their
paper the authors referred to a private communication by the present author.
The purpose of the present paper is to make this communication available.

For a normalized Stieltjes moment sequence $s$ which is indeterminate as a
Hamburger moment sequence, there is a quantity $\alpha=\alpha(s)\leq 0$, which
can be used to characterize Stieltjes indeterminacy of $s$. In fact there is
exactly one measure (resp. several measures) supported by $[0,\infty)$ with
moment sequence $s$ if $\alpha(s)=0$ (resp. $\alpha(s)<0$). Chihara [7, Lemma
1] proved an equivalent of

$\alpha(s)=\lim_{n\to\infty}\frac{P_{n}(0)}{Q_{n}(0)}$ (1)

in terms of chain sequences, where $P_{n}$ and $Q_{n}$ are the orthonormal
polynomials and those of the second kind affiliated with $s$. In [5] there was
given another characterization of $\alpha(s)$ as

$\alpha(s)=\lim_{x\to-\infty}\frac{D(x)}{B(x)},$ (2)

where $B,D$ are two of the Nevanlinna functions for the indeterminate
Hamburger problem, see Section 2 for definitions in connection with the
Nevanlinna parametrization (25) of the solutions to the indeterminate
Hamburger moment problem. The formula (2) for $\alpha(s)$ together with the
graph of the function $j(x)=D(x)/B(x)$ in Section 3 makes it transparent that
the zeros of $B(z)t-D(z)$ are bounded below, hence of the form
$x_{1}(t)<x_{2}(t)<\ldots$ for a sequence $(x_{n}(t))$ tending to infinity.
When $t=\infty$ then $B(z)t-D(z)$ shall be interpreted as $B(z)$. Therefore,
the N-extremal solutions $\nu_{t},t\in\mathbb{R}\cup\\{\infty\\}$ of the
Hamburger moment problem are of the form

$\nu_{t}=\sum_{n=1}^{\infty}\rho(x_{n}(t))\varepsilon_{x_{n}(t)},$ (3)

where

$\rho(x)=\left(\sum_{n=0}^{\infty}P_{n}^{2}(x)\right)^{-1},\quad
x\in\mathbb{R}.$ (4)

By $\varepsilon_{a}$ we denote the measure with mass 1 concentrated in the
point $a\in\mathbb{R}$.

For any solution $\mu$ of the Hamburger moment problem we define

$\xi(\mu):=\inf\operatorname{supp}(\mu).$ (5)

We therefore get

$\operatorname{supp}(\nu_{t})=\\{x_{n}(t),n\geq
1\\},\quad\xi(\nu_{t})=x_{1}(t).$

In particular, the graph of $j(x)$ shows that $\operatorname{supp}(\nu_{t})$
is contained in $[0,\infty)$ if and only if $\alpha(s)\leq t\leq 0$. The last
result can also be deduced from [8, p. 340]. Notice that there are solutions
$\mu$ with $\operatorname{supp}(\mu)=\mathbb{R}$ and hence $\xi(\mu)=-\infty$.
This is true for the solution

$\mu=\frac{1}{\pi}\left(B^{2}(x)+D^{2}(x)\right)^{-1}dx,$

cf. formula (2.15) in [5].

The standard orthonormal basis of $\ell^{2}$ is denoted $e_{0},e_{1},\ldots$.
We consider the Jacobi matrix (14) as an operator in $\ell^{2}$ with the
domain $\mathcal{F}$ of sequences with only finitely many non-zero elements.
It is a symmetric operator of deficiency indices $(1,1)$. Its self-adjoint
extensions in $\ell^{2}$ can be parametrized as
$(T_{t},D(T_{t})),t\in\mathbb{R}\cup\\{\infty\\}$ such that

$\nu_{t}=\langle E_{t}e_{0},e_{0}\rangle,t\in\mathbb{R}\cup\\{\infty\\},$ (6)

where $E_{t}$ is the spectral measure of $T_{t}$. The operators $T_{t}$ are
bounded below with the lower bound

$\xi(\nu_{t})=x_{1}(t)=\inf\\{\langle T_{t}e,e\rangle\mid e\in
D(T_{t}),||e||=1\\}.$ (7)

Later the two N-extremal solutions $\nu_{0}$ and $\nu_{\alpha(s)}$ were
identified with respectively the Krein and the Friedrichs self-adjoint
extensions of the Jacobi matrix. As far as we know the measure
$\nu_{\alpha(s)}$ was first associated with the Friedrichs extension in
Pedersen’s paper [11].

For information about the Friedrichs extension in general see [9]. It has the
same lower bound as the Jacobi matrix, viz.

$x_{1}(\alpha(s))=\inf\\{\langle Je,e\rangle\mid e\in\mathcal{F},||e||=1\\}.$

######  Remark 1.1.

The self-adjoint operators $(T_{t},D(T_{t}))$ are parametrized differently in
[14, Theorem 6.23] because of a different form of the Nevanlinna
parametrization. See Remark 7.14 in [14].

We prove that $\alpha(s)$ can also be calculated as

$\alpha(s)=\lim_{x\to-\infty}\frac{C(x)}{A(x)},$

where $A,C$ are the two other Nevanlinna functions.

Let $x_{n,k},k=1,\dots,n$ be the zeros of $P_{n}$ in increasing order. By the
interlacing of the zeros of $P_{n}$ and $P_{n+1}$ we have that
$(x_{n,k})_{n\geq k}$ is decreasing in $n$ for each $k\geq 1$, and therefore
the following limit exists for each $k\geq 1$:

$\xi_{k}:=\lim_{n\to\infty}x_{n,k},$ (8)

and clearly $0\leq\xi_{1}\leq\xi_{2}\leq\ldots$. These numbers appear on the
graph in Section 3.

In [7, Lemma 2] Chihara proved the existence of an N-extremal measure
$\psi_{*}$, which is concentrated on the sequence $(\xi_{k})$. For the benefit
of the reader we give a complete proof in Theorem 3.4. Chihara also proved
that $\psi_{*}$ has the extremality property, that any other solution
$\mu\neq\psi_{*}$ of the indeterminate Hamburger moment problem satisfies

$\xi(\mu)<\xi(\psi_{*})=\xi_{1}.$ (9)

Chihara does not mention the Friedrichs extension, but (9) is essential in
Pedersen’s proof that $\psi_{*}$ is the Friedrichs extension, because it has
the greatest lower bound among the self-adjoint extensions of the Jacobi
matrix.

Chihara’s proofs are largely based on chain sequences. The result is given in
Theorem 3.5 with a proof independent of chain sequences. The measure
$\psi_{*}$ equals $\nu_{\alpha(s)}$, so
$\xi_{k}=x_{k}(\alpha(s)),k=1,2,\ldots$ with the notation of (3).

In Section 4 we consider the moment problem corresponding to the shifted
Jacobi matrix $J^{\dagger}$ formed by the shifted Jacobi parameters
$(a_{k+1}),(b_{k+1})$. This means that $J^{\dagger}$ is obtained from $J$ by
deleting the first row and column of $J$. The corresponding normalized
Stieltjes moment sequence $s^{\dagger}$ is always indeterminate as Stieltjes
moment sequence because of (45). This has also been obtained in [2, Theorem
4.3] by another proof. In Theorem 4.1 we finally prove the interlacing of the
sequences $(\xi_{k})$ and $(\xi^{\dagger}_{k})$ associated with $s$ and
$s^{\dagger}$ via (8).

##  2 Preliminaries about indeterminate Hamburger moment problems

The starting point is a normalized indeterminate Hamburger moment sequence
$s=(s_{n})_{n\geq 0}$ as in Section 2.1 of [5], i.e.,

$s_{n}=\int x^{n}\,d\mu(x),\quad\mu\in V,n=0,1,\ldots,$ (10)

where $V$ is the infinite convex set of probability measures on $\mathbb{R}$
with moments of any order and having the moments $s_{n}$. As in [5] we denote
by $P_{n}(x)$ and $Q_{n}(x)$ the orthonormal polynomials and the corresponding
polynomials of the second kind. Both polynomial sequences satisfy the
difference equation

$xy_{n}=b_{n}y_{n+1}+a_{n}y_{n}+b_{n-1}y_{n-1},\quad n\geq 1$ (11)

together with the initial conditions

$P_{0}(x)=1,P_{1}(x)=\frac{1}{b_{0}}(x-a_{0}),\quad
Q_{0}(x)=0,Q_{1}(x)=\frac{1}{b_{0}}.$ (12)

Here

$a_{n}=\int xP_{n}^{2}(x)\,d\mu(x),\quad b_{n}=\int
xP_{n}(x)P_{n+1}(x)\,d\mu(x),\quad\mu\in V.$ (13)

The infinite matrix

$J=\begin{pmatrix}a_{0}&b_{0}&0&\ldots\\\ b_{0}&a_{1}&b_{1}&\ldots\\\
0&b_{1}&a_{2}&\ldots\\\ \vdots&\vdots&\vdots&\ddots\end{pmatrix},$ (14)

is called the Jacobi matrix of the moment problem. The matrix $J$ defines an
operator in $\ell^{2}$ with the domain $\mathcal{F}$ of sequences with only
finitely many non-zero elements. The closure of $J$ is a symmetric operator
$(T,D(T))$ in $\ell^{2}$ with deficiency indices $(1,1)$ called the Jacobi
operator of $s$. For results about the domain $D(T)$ of the Jacobi operator,
and the domains $D(T_{t})$ of its self-adjoint extensions see [4].

We combine the polynomials $P_{n},Q_{n}$ to form polynomials of two variables,
and the following formulas hold:

######  Proposition 2.1.

[14, Proposition 5.24] For $u,v\in\mathbb{C}$ and $n\geq 0$ we have

$\displaystyle A_{n}(u,v)$ $\displaystyle:=$
$\displaystyle(u-v)\sum_{k=0}^{n}Q_{k}(u)Q_{k}(v)=b_{n}\left|\begin{array}[]{cc}Q_{n+1}(u)&\;Q_{n+1}(v)\\\
Q_{n}(u)&\;Q_{n}(v)\end{array}\right|$ $\displaystyle B_{n}(u,v)$
$\displaystyle:=$
$\displaystyle-1+(u-v)\sum_{k=0}^{n}P_{k}(u)Q_{k}(v)=b_{n}\left|\begin{array}[]{cc}P_{n+1}(u)&\;Q_{n+1}(v)\\\
P_{n}(u)&\;Q_{n}(v)\end{array}\right|$ $\displaystyle C_{n}(u,v)$
$\displaystyle:=$ $\displaystyle
1+(u-v)\sum_{k=0}^{n}Q_{k}(u)P_{k}(v)=b_{n}\left|\begin{array}[]{cc}Q_{n+1}(u)&\;P_{n+1}(v)\\\
Q_{n}(u)&\;P_{n}(v)\end{array}\right|$ $\displaystyle D_{n}(u,v)$
$\displaystyle:=$
$\displaystyle(u-v)\sum_{k=0}^{n}P_{k}(u)P_{k}(v)=b_{n}\left|\begin{array}[]{cc}P_{n+1}(u)&\;P_{n+1}(v)\\\
P_{n}(u)&\;P_{n}(v)\end{array}\right|.$

It is important to notice that

$\left|\begin{array}[]{cc}\ A_{n}(u,v)&\;B_{n}(u,v)\\\
C_{n}(u,v)&\;D_{n}(u,v)\end{array}\right|=1\textrm{ for
}(u,v)\in\mathbb{C}^{2},$

cf. [14, Equation(5.57)]. These polynomials were also introduced in [1, p.
123].

Since we assume $(s_{n})$ to be indeterminate , we know that
$(P_{n}(u)),(Q_{n}(u))\in\ell^{2}$ for all $u\in\mathbb{C}$, so we can let $n$
tend to infinity in the expressions in the middle of Proposition 2.1, and we
get four entire functions of two complex variables called the Nevanlinna
functions of the indeterminate Hamburger moment problem:

$\displaystyle A(u,v)$ $\displaystyle=$
$\displaystyle(u-v)\sum_{k=0}^{\infty}Q_{k}(u)Q_{k}(v)$ (19) $\displaystyle
B(u,v)$ $\displaystyle=$
$\displaystyle-1+(u-v)\sum_{k=0}^{\infty}P_{k}(u)Q_{k}(v)$ (20) $\displaystyle
C(u,v)$ $\displaystyle=$ $\displaystyle
1+(u-v)\sum_{k=0}^{\infty}Q_{k}(u)P_{k}(v)$ (21) $\displaystyle D(u,v)$
$\displaystyle=$ $\displaystyle(u-v)\sum_{k=0}^{\infty}P_{k}(u)P_{k}(v),$ (22)

see Section 7.1 in [14]. The two-variable functions were introduced in [6] in
a slightly different form, which was subsequently used in [3], [10]. If the
functions of [6] are marked with a $*$, we have

$\displaystyle A^{*}(u,v)$ $\displaystyle=$
$\displaystyle-A(u,v),\;B^{*}(u,v)=-C(u,v),$ $\displaystyle C^{*}(u,v)$
$\displaystyle=$ $\displaystyle-B(u,v),\;D^{*}(u,v)=-D(u,v).$

We clearly have the determinant equation

$\left|\begin{array}[]{cc}\ A(u,v)&\;B(u,v)\\\
C(u,v)&\;D(u,v)\end{array}\right|=1\textrm{ for }(u,v)\in\mathbb{C}^{2}.$

We define the following polynomials and entire functions of one variable

$A_{n}(u)=A_{n}(u,0),\;B_{n}(u)=B_{n}(u,0),\;C_{n}(u)=C_{n}(u,0),\;D_{n}(u)=D_{n}(u,0),$
(23) $A(u)=A(u,0),\;B(u)=B(u,0),\;C(u)=C(u,0),\;D(u)=D(u,0).$ (24)

The latter are also called the Nevanlinna functions of the moment problem, and
they enter in the Nevanlinna parametrization of the indeterminate Hamburger
moment problem,
$V=\\{\nu_{\varphi}\mid\varphi\in\mathcal{N}\cup\\{\infty\\}\\}$:

$\int\frac{d\nu_{\varphi}(x)}{x-z}=-\frac{A(z)\varphi(z)-C(z)}{B(z)\varphi(z)-D(z)},\quad
z\in\mathbb{C}\setminus\mathbb{R},$ (25)

where $\mathcal{N}$ is the set of Pick functions, i.e., the holomorphic
functions $\varphi:\mathbb{C}\setminus\mathbb{R}\to\mathbb{C}$ satisfying
$\Im(\varphi(z))/\Im(z)\geq 0$. See [1] or [5] for details. When the Pick
function $\varphi\in\mathcal{N}$ is a constant
$t\in\mathbb{R}\cup\\{\infty\\}$, the solution $\nu_{t}$ is called N-extremal
or a von Neumann solution. It is well-known that the entire functions
$A,B,C,D$ have infinitely many zeros which are all real and simple.

The two-variable polynomials of Proposition 2.1 can be expressed in terms of
the one variable polynomials (23) as follows, see [1, p. 123], or [4]:

######  Proposition 2.2.

For $u,v\in\mathbb{C}$ we have

$\displaystyle A_{n}(u,v)$ $\displaystyle=$
$\displaystyle\left|\begin{array}[]{cc}A_{n}(u)&\;A_{n}(v)\\\
C_{n}(u)&\;C_{n}(v)\end{array}\right|$ $\displaystyle B_{n}(u,v)$
$\displaystyle=$ $\displaystyle\left|\begin{array}[]{cc}B_{n}(u)&\;A_{n}(v)\\\
D_{n}(u)&\;C_{n}(v)\end{array}\right|$ $\displaystyle C_{n}(u,v)$
$\displaystyle=$ $\displaystyle\left|\begin{array}[]{cc}A_{n}(u)&\;B_{n}(v)\\\
C_{n}(u)&\;D_{n}(v)\end{array}\right|$ $\displaystyle D_{n}(u,v)$
$\displaystyle=$ $\displaystyle\left|\begin{array}[]{cc}B_{n}(u)&\;B_{n}(v)\\\
D_{n}(u)&\;D_{n}(v)\end{array}\right|.$

There are similar expressions for the Nevanlinna functions of one and two
variables just by letting $n$ tend to infinity.

Differentiating the formula for $A(u,v)$ with respect to $u$ and setting
$u=v=x$ we get

$A^{\prime}(x)C(x)-C^{\prime}(x)A(x)=\sum_{k=0}^{\infty}Q_{k}^{2}(x),$ (30)

and the same procedure for $D(u,v)$ leads to

$B^{\prime}(x)D(x)-D^{\prime}(x)B(x)=\sum_{k=0}^{\infty}P_{k}^{2}(x).$ (31)

For zeros $\alpha$ of $A$ and $\gamma$ of $C$ we get from (30)

$A^{\prime}(\alpha)C(\alpha)>0,\quad C^{\prime}(\gamma)A(\gamma)<0,$

which shows that $C$ has a zero between any two consecutive zeros of $A$, and
that $A$ has a zero between any two consecutive zeros of $C$. In other words
the zeros of $A$ and $C$ are interlacing. In the same way (31) shows that the
zeros of $B$ and $D$ are interlacing.

##  3 The Stieltjes case

From now on we assume that $s=(s_{n})$ in addition to the previous assumptions
is a Stieltjes moment sequence, so there exists at least one measure $\mu\in
V$ supported by the half-line $[0,\infty)$.

There are two possibilities for a Stieltjes moment problem which is
indeterminate as Hamburger moment problem. There can be precisely one solution
$\mu\in V$ supported by $[0,\infty)$, or there can be more than one such
solution and then infinitely many, since they form a convex set. We say that
the Stieltjes moment problem is determinate, respectively indeterminate in the
sense of Stieltjes in the two cases. In short the problem is called det(S) or
indet(S).

In [5, Section 4] there is a detailed study of the Al-Salam–Carlitz
polynomials $V_{n}^{(a)}(x;q)$, which depends on parameters $a>0$ and $0<q<1$.
There is given a complete classification of when the corresponding Hamburger
moment problem is determinate or indeterminate, and when the indeterminate
Hamburger problem is det(S) or indet(S).

Since the zeros of $P_{n}$ and $Q_{n}$ are located in $(0,\infty)$ and the
leading coefficients of $P_{n}$ and $Q_{n}$ are positive, we see that for
$k\geq 1,x\leq 0$

$P_{k}(x)P_{k}(0)>0,\;Q_{k}(x)Q_{k}(0)>0,\;P_{k}(x)Q_{k}(0)<0,\;Q_{k}(x)P_{k}(0)<0,$
(32)

while for $k=0$ the first quantity equals 1 and the three others equal 0. It
follows that $B_{n}(x)<B_{n+1}(x)<B(x),D_{n}(x)>D_{n+1}(x)>D(x)$ for $x<0$,
and in particular for $x\leq 0$:

$B(x)\geq B_{1}(x)=-1+xP_{1}(x)Q_{1}(0)=\frac{1}{b_{0}^{2}}x(x-a_{0})-1,\quad
D(x)\leq D_{1}(x)=x.$ (33)

From this equation we see that $B$ has a smallest zero $\beta_{1}<0$ and that
$D$ has a smallest zero $\delta_{1}=0$. This means that we can write the zeros
$\beta_{n},\delta_{n},n\geq 1$ of $B$ and $D$ in increasing order, and by the
interlacing property mentioned at the end of Section 2 we have

$\beta_{1}<\delta_{1}=0<\beta_{2}<\delta_{2}<\beta_{3}<\ldots,$

and $\beta_{n},\delta_{n}$ tend to infinity.

From (32) we similarly get $A_{n}(x)>A_{n+1}(x)>A(x),C_{n}(x)<C_{n+1}(x)<C(x)$
for $x<0$, and in particular for $x\leq 0$:

$A(x)\leq A_{1}(x)=\frac{x}{b_{0}^{2}},\quad C(x)\geq
C_{1}(x)=1-\frac{a_{0}}{b_{0}^{2}}x\geq 1.$ (34)

From this equation we see that $A$ has a smallest zero $\alpha_{1}=0$ and that
$C$ has a smallest zero $\gamma_{1}>0$. This means that we can write the zeros
$\alpha_{n},\gamma_{n},n\geq 1$ of $A$ and $C$ in increasing order, and by the
interlacing property we have

$\alpha_{1}=0<\gamma_{1}<\alpha_{2}<\gamma_{2}<\alpha_{3}<\ldots,$

and $\alpha_{n},\gamma_{n}$ tend to infinity.

Based on previous work of Chihara, see [7],[8], the following result was
established in [2] and [5, section 2.2]:

######  Proposition 3.1.

Let $s=(s_{n})$ be a normalized indeterminate Hamburger moment sequence which
is Stieltjes. The function $j(x)=D(x)/B(x)$ is strictly decreasing in each of
the intervals
$(-\infty,\beta_{1}),(\beta_{1},\beta_{2}),(\beta_{2},\beta_{3}),\ldots$ with
limits

$j(\beta_{k}-)=-\infty,\quad j(\beta_{k}+)=\infty,\quad k=1,2,\ldots.$

Furthermore, the quantity $\alpha=\alpha(s)$ given by

$\alpha=\lim_{x\to-\infty}\frac{D(x)}{B(x)}=\lim_{n\to\infty}\frac{P_{n}(0)}{Q_{n}(0)}$
(35)

belongs to $(-\infty,0]$, and the sequence $P_{n}(0)/Q_{n}(0)$ is strictly
increasing.

The Stieltjes problem is det(S) if $\alpha(s)=0$ and indet(S) if
$\alpha(s)<0$.

The graph of $j(x)=D(x)/B(x)$

The quantity $\alpha(s)$ can also be determined from the functions $A,C$:

######  Proposition 3.2.

Let $s=(s_{n})$ be a normalized indeterminate Hamburger moment sequence which
is Stieltjes. The function $l(x)=C(x)/A(x)$ is strictly decreasing in each of
the intervals
$(-\infty,\alpha_{1}),(\alpha_{1},\alpha_{2}),(\alpha_{2},\alpha_{3}),\ldots$
with limits

$l(\alpha_{k}-)=-\infty,\quad l(\alpha_{k}+)=\infty,\quad k=1,2,\ldots.$

We have

$\alpha(s)=\lim_{x\to-\infty}\frac{C(x)}{A(x)}.$ (36)

###### Proof.

We know from Berg–Szwarc [4] that $A(z)/C(z)$ is a Pick function so
$C(x)/A(x)$ is strictly decreasing in the intervals
$(-\infty,\alpha_{1}),(\alpha_{k},\alpha_{k+1}),k=1,2,\ldots$. (That
$A(z)/C(z)$ is a Pick function is also deduced in Section 4 of the present
paper). We have $A(x)<0,B(x)>0$ for $x<\beta_{1}$, so for those $x$ we have
$D(x)/B(x)-C(x)/A(x)=1/(A(x)B(x))\to 0$ for $x\to-\infty$ because
$A(x)\to-\infty$ and $B(x)\to\infty$ by (34) and (33). This shows (36). ∎

If the Stieltjes problem is det(S), i.e., $\alpha(s)=0$ then $\nu_{0}$ is the
only solution to the Hamburger problem which is concentrated on $[0,\infty[$.
For all other solutions $\mu$ we have $\xi(\mu)<0$.

Let us now assume that the Stieltjes problem is indet(S), i.e., that
$\alpha(s)<0$.

For $t\in\mathbb{R}\cup\\{\infty\\}$ the N-extremal solutions $\nu_{t}\in V$
are the discrete measures satisfying

$\int\frac{d\nu_{t}(x)}{x-z}=-\frac{A(z)t-C(z)}{B(z)t-D(z)},\quad
z\in\mathbb{C}\setminus\mathbb{R},$ (37)

so $\nu_{t}$ is concentrated in the zero-set

$\Lambda_{t}:=\\{x\in\mathbb{R}\mid B(x)t-D(x)=0\\}=\\{x\in\mathbb{R}\mid
j(x)=t\\},$

which can be described as a strictly increasing sequence $(x_{n}(t))_{n\geq
1}$. From Proposition 3.1 we clearly have:

(i) For $t>0$: $x_{n}(t)\in(\beta_{n},\delta_{n})$, so $\nu_{t}$ has one
negative mass-point $x_{1}(t)$ and the other mass-points are positive. When
$t$ increases from 0 to $\infty$, then $x_{1}(t)$ decreases from
$\delta_{1}=0$ to $\beta_{1}$ with $x_{1}(\infty)=\beta_{1}$.

(ii) For $t<\alpha(s)$:
$x_{1}(t)\in(-\infty,\beta_{1}),\;x_{n}(t)\in(\delta_{n-1},\beta_{n})$ for
$n\geq 2$, so $\nu_{t}$ has one negative mass-point and the other mass-points
are positive. When $t$ increases from $-\infty$ to $\alpha(s)$, then
$x_{1}(t)$ decreases from $\beta_{1}$ to $-\infty$.

(iii) For $t=\infty$: $x_{n}(t)=\beta_{n}$, so $\nu_{\infty}$ is concentrated
in the zero-set of $B$ and has one negative mass-point.

(iv) For $\alpha(s)<t<0$: $x_{n}(t)\in(\delta_{n},\beta_{n+1})$, so $\nu_{t}$
is supported by a sequence belonging to $(0,\infty)$. When $t$ increases from
$\alpha(s)$ to 0, then $x_{1}(t)$ decreases from $\xi_{1}$ to 0.

(v) For $t=0$: $x_{n}(t)=\delta_{n}$, so $\nu_{0}$ is supported by the zero-
set of $D$. It is called the Krein solution of the moment problem.

(vi) For $t=\alpha(s)$: $x_{n}(t)>\delta_{n}$, and $x_{1}(\alpha(s))>x_{1}(t)$
for all $t\neq\alpha(s)$ (including $t=\infty$). The measure $\nu_{\alpha(s)}$
is called the Friedrichs solution of the moment problem.

That the Friedrichs solution $\nu_{\alpha(s)}$ corresponds to the Friedrichs
extension of the Jacobi matrix $J$ or the Jacobi operator $(T,D(T))$ was
established in [11]. In the paper [12] Pedersen gave a very satisfactory
characterization of the solutions $\nu_{\varphi}\in V$ which are supported by
$[0,\infty)$, in terms of the corresponding Pick function $\varphi$. For the
benefit of the reader we repeat it here.

######  Theorem 3.3.

[12] The solution $\nu_{\varphi}\in V$ corresponding to the Pick function
$\varphi$ is supported by $[0,\infty)$ if an only if $\varphi$ has a
holomorphic extension from $\mathbb{C}\setminus\mathbb{R}$ to
$\mathbb{C}\setminus[0,\infty)$ such that

$\alpha(s)\leq\varphi(x)\leq 0,\quad x<0.$

We shall now describe the Friedrichs solution $\nu_{\alpha(s)}$ in more
detail.

######  Theorem 3.4.

The sequence $(\xi_{k})$ defined in (8) satisfies
$0<\xi_{1}<\xi_{2}<\ldots<\xi_{k}\to\infty$ and the $\xi_{k}$ are the zeros of
$B(z)\alpha(s)-D(z)$. The support of the Friedrichs solution is $(\xi_{k})$.

###### Proof.

It is elementary to see that

$-\frac{P_{n}(z)}{Q_{n}(0)}\to B(z)\alpha(s)-D(z),\quad z\in\mathbb{C},$ (38)

locally uniformly in $z$, cf. [11, Equation (8)], and hence also the
derivatives $-P^{\prime}_{n}(z)/Q_{n}(0)$ converge to the derivative
$B^{\prime}(z)\alpha(s)-D^{\prime}(z)$. Define $y_{n,k}$ such that
$x_{n,k}<y_{n,k}<x_{n,k+1}$ be such that $P_{n}^{\prime}(y_{n,k})=0$ for
$k=1,\ldots,n-1$.

From (38) it is clear that each $\xi_{k}$ is a zero of $B(z)\alpha(s)-D(z)$.
We claim that $\xi_{k}<\xi_{k+1}$, for if there was equality we would have
$\lim_{n\to\infty}y_{n,k}=\xi_{k}$, and hence
$B^{\prime}(\xi_{k})\alpha(s)-D^{\prime}(\xi_{k})=0$, contradicting that the
zeros of $B(z)\alpha(s)-D(z)$ are simple. In addition $\xi_{k}\to\infty$ for
otherwise the zeros accumulate.

Finally $B(z)\alpha(s)-D(z)$ has no other zeros than $\xi_{k}$, for if $\xi$
is a real zero different from all $\xi_{k}$, then there exists an
$\varepsilon>0$ such that

$[\xi-\varepsilon,\xi+\varepsilon]\cap[\xi_{k}-\varepsilon,\xi_{k}+\varepsilon]=\emptyset,\forall
k.$ (39)

By the Theorem of Rouché $P_{n}$ must have zeros in
$(\xi-\varepsilon,\xi+\varepsilon)$ for $n$ sufficiently large, but this is
impossible because of (39).

Clearly $\xi_{1}\geq 0$, but if $\xi_{1}=0$ then $D(\xi_{1})=0$, hence
$B(\xi_{1})\alpha(s)=0$, which gives the contradiction $B(\xi_{1})=0$ since we
assume $\alpha(s)<0$. ∎

We know that the Friedrich measure $\nu_{\alpha(s)}$ satisfies

$\xi_{1}=\xi(\nu_{\alpha(s)}),\quad\alpha(s)=D(\xi_{1})/B(\xi_{1}).$ (40)

We have the following result about $\xi_{1}$, which is the same as [7, Lemma
2], except that there is no mentioning of the Friedrichs extension.

######  Theorem 3.5.

Let $s$ be a normalized Stieltjes moment sequence which is indet(S). For any
solution $\nu\in V$ to the indeterminate Hamburger problem we have
$\xi(\nu)<\xi(\nu_{\alpha(s)})$ unless $\nu$ is the Friedrichs solution
$\nu_{\alpha(s)}$.

###### Proof.

Suppose that $\nu$ is a solution satisfying
$\operatorname{supp}(\nu)\subseteq[\xi_{1},\infty)$, and we want to show that
$\nu=\nu_{\alpha(s)}$.

Considering the translation $\tau(x)=x-\xi_{1}$, the image measures
$\tau(\nu_{\alpha(s)}),\tau(\nu)$ are supported on $[0,\infty)$ and belong to
the Stieltjes moment problem with moments

$\tilde{s}_{n}=\int(x-\xi_{1})^{n}\,d\nu(x)=\sum_{j=0}^{n}\binom{n}{j}(-\xi_{1})^{j}s_{n-j}.$

It is clearly indeterminate as Hamburger moment problem and the full set of
solutions is $\\{\tau(\mu)\mid\mu\in V\\}$. The orthonormal polynomials and
those of the second kind are $P_{n}(x+\xi_{1}),Q_{n}(x+\xi_{1})$. The
corresponding Nevanlinna functions marked with $\tilde{}$ are

$\tilde{A}(u,v)=A(u+\xi_{1},v+\xi_{1}),\ldots,\tilde{D}(u,v)=D(u+\xi_{1},v+\xi_{1}),\quad
u,v\in\mathbb{C}$

From the formulas for the Nevanlinna functions of two variables expressed by
the Nevanlinna functions of one variable, see Proposition 2.2 (where
$n\to\infty$), we get

$\displaystyle\tilde{A}(z)$ $\displaystyle=$ $\displaystyle
A(z+\xi_{1})C(\xi_{1})-C(z+\xi_{1})A(\xi_{1})$ $\displaystyle\tilde{B}(z)$
$\displaystyle=$ $\displaystyle B(z+\xi_{1})C(\xi_{1})-D(z+\xi_{1})A(\xi_{1})$
$\displaystyle\tilde{C}(z)$ $\displaystyle=$ $\displaystyle
A(z+\xi_{1})D(\xi_{1})-C(z+\xi_{1})B(\xi_{1})$ $\displaystyle\tilde{D}(z)$
$\displaystyle=$ $\displaystyle
B(z+\xi_{1})D(\xi_{1})-D(z+\xi_{1})B(\xi_{1}).$

These formulas were found in [11, Proposition 3.3]. According to these
formulas we get using (40)

$\frac{\tilde{D}(x)}{\tilde{B}(x)}=\frac{D(\xi_{1})-B(\xi_{1})\frac{D(x+\xi_{1})}{B(x+\xi_{1})}}{C(\xi_{1})-A(\xi_{1})\frac{D(x+\xi_{1})}{B(x+\xi_{1})}}\to\frac{D(\xi_{1})-B(\xi_{1})\alpha(s)}{C(\xi_{1})-A(\xi_{1})\alpha(s)}=0$

for $x\to-\infty$, and this shows that $\tilde{\alpha}:=\alpha(\tilde{s})=0$
for the new Stieltjes sequence $\tilde{s}$, so $\tilde{s}$ is det(S). We
conclude that $\tau(\nu_{\alpha(s)})=\tau(\nu)$ and hence
$\nu_{\alpha(s)}=\nu$. ∎

######  Remark 3.6.

Let $s$ be a normalized indeterminate Hamburger moment sequence and assume
that there exists a solution $\mu$ with $\xi(\mu)\in\mathbb{R}$. Then there
exists a solution $\mu_{*}$ such that $\xi(\mu)<\xi(\mu_{*})$ for all
solutions $\mu\neq\mu_{*}$.

This follows by using the translation $\tau(x)=x-\xi(\mu)$ to obtain a
Stieltjes moment sequence.

##  4 The shifted problem

If the first row and column of the Jacobi matrix (14) is deleted, we obtain
the matrix

$J^{\dagger}=\begin{pmatrix}a_{1}&b_{1}&0&\ldots\\\
b_{1}&a_{2}&b_{2}&\ldots\\\ 0&b_{2}&a_{3}&\ldots\\\
\vdots&\vdots&\vdots&\ddots\end{pmatrix},$ (41)

which by Favard’s Theorem, cf. [14, Theorem 5.14], is the Jacobi matrix of a
uniquely determined normalized moment sequence , which we denote
$s^{\dagger}=(s_{n}^{\dagger})$ and call the shifted moment sequence. The
shifted problem has been studied by many authors, and we shall use results
from [2] and [10].

The original moment problem is Stieltjes if and only if $J$ is positive in the
sense that $\langle Jc,c\rangle\geq 0$ for all sequences
$c=(c_{n})\in\mathcal{F}$. Therefore, if $s$ is Stieltjes, then so is
$s^{\dagger}$.

The orthonormal polynomials and the associated polynomials for $s^{\dagger}$
are given as

$P^{\dagger}_{n}(x)=b_{0}Q_{n+1}(x),\quad
Q^{\dagger}_{n}(x)=P_{1}(x)Q_{n+1}(x)-\frac{1}{b_{0}}P_{n+1}(x),$ (42)

and these formulas show that $s$ is indeterminate as a Hamburger problem if
and only if $s^{\dagger}$ is so. They are in [10] as well as expressions for
the Nevanlinna functions for the shifted problem in terms of those for the
original problem. In particular Pedersen found the formulas

$B^{\dagger}(z)=-C(z)-a_{0}A(z),\quad D^{\dagger}(z)=b_{0}^{2}A(z).$ (43)

(From this we get $b_{0}^{2}B^{\dagger}(z)/D^{\dagger}(z)=-a_{0}-C(z)/A(z)$,
and from this we see that $A(z)/C(z)$ is a Pick function, a property used in
Proposition 3.2.)

In the following we assume that $s$ is a Stieltjes sequence for an
indeterminate Hamburger problem, and we shall relate $\alpha(s)$ and
$\alpha(s^{\dagger})$.

First of all we know from Proposition 3.1 that the sequence
$P_{n}(0)/Q_{n}(0)$ is strictly increasing to $\alpha(s)$. In particular

$P_{1}(0)/Q_{1}(0)=\frac{-a_{0}/b_{0}}{1/b_{0}}=-a_{0}<\alpha(s),$

and hence

$\alpha(s)+a_{0}>0.$ (44)

Since $P^{\dagger}_{n}(x)=b_{0}Q_{n+1}(x)$, we know that the zeros $x_{n,k}$
of $P_{n}$ and the zeros $x^{\dagger}_{n,k}$ of $P^{\dagger}_{n}$ satisfy

$x_{n,1}<x^{\dagger}_{n-1,1}<x_{n,2}<x^{\dagger}_{n-1,2}\ldots<x^{\dagger}_{n-1,n-1}<x_{n,n},$

and therefore with the notation from Theorem 3.4

$\xi_{1}\leq\xi^{\dagger}_{1}\leq\xi_{2}\leq\xi^{\dagger}_{2}\leq\ldots.$

We claim that there are strict inequalities everywhere.

We notice that

$\alpha(s^{\dagger})=-\frac{b_{0}^{2}}{\alpha(s)+a_{0}}<0,$ (45)

so the shifted problem is always indet(S), an observation also done in [2,
Remark 4.5].

In fact, by (43) we find

$\alpha(s^{\dagger})=\lim_{x\to-\infty}\frac{D^{\dagger}(x)}{B^{\dagger}(x)}=\lim_{x\to-\infty}\frac{b_{0}^{2}A(x)}{-C(x)-a_{0}A(x)},$

and (45) follows because $C(x)/A(x)\to\alpha(s)$ for $x\to-\infty$ by (36).

We know that $\xi^{\dagger}_{k}$ are the zeros of

$B^{\dagger}(z)\alpha(s^{\dagger})-D^{\dagger}(z)=-\frac{\alpha(s)b_{0}^{2}}{\alpha(s)+a_{0}}A(z)+\frac{b_{0}^{2}}{\alpha(s)+a_{0}}C(z).$

If $\alpha(s)=0$ this function is proportional to $C(z)$ so $\xi_{k}$ are the
zeros of $D$ and $\xi^{\dagger}_{k}$ are the zeros of $C$ and hence
$\xi_{k}\neq\xi^{\dagger}_{k}$. If $\alpha(s)<0$ then
$B^{\dagger}(z)\alpha(s^{\dagger})-D^{\dagger}(z)$ is proportional to
$A(z)\alpha(s)-C(z)$, which then has the zeros $\xi^{\dagger}_{k}$, while
$B(z)\alpha(s)-D(z)$ has the zeros $\xi_{k}$, and again
$\xi_{k}\neq\xi^{\dagger}_{k}$.

Summing up we have proved:

######  Theorem 4.1.

Let $s$ be a normalized indeterminate Hamburger moment sequence which is
Stieltjes. The zeros $x_{n,k}$ of $P_{n}$ and the zeros $x^{\dagger}_{n,k}$ of
$Q_{n+1}$ have limits

$\xi_{k}:=\lim_{n\to\infty}x_{n,k},\quad\xi^{\dagger}_{k}:=\lim_{n\to\infty}x^{\dagger}_{n,k},\quad
k=1,2,\ldots$

satisfying

$\xi_{1}<\xi^{\dagger}_{1}<\xi_{2}<\xi^{\dagger}_{2}<\ldots.$

The zeros of $B(z)\alpha(s)-D(z)$ are $(\xi_{k})$ and the zeros of
$A(z)\alpha(s)-C(z)$ are $(\xi^{\dagger}_{k})$.

######  Remark 4.2.

It is a general result about indeterminate Hamburger problems that for each
$t\in\mathbb{R}$ the entire functions $A(z)t-C(z)$ and $B(z)t-D(z)$ have
infinitely many zeros which are all real and simple and they interlace.

## References

  * [1] N. I. Akhiezer, The classical moment problem. Oliver & Boyd, Edinburgh, 1965.
  * [2] C. Berg, Markov’s Theorem revisited, J. Approx. Theory 78 (1994), 260–275.
  * [3] C. Berg, Indeterminate moment problems and the theory of entire functions, J. Comput. Appl. Math. 65 (1995), 27–55.
  * [4] C. Berg, R. Szwarc, Indeterminate Jacobi operators. To appear in J. Operator Theory. ArXiv:2301.00586.
  * [5] C. Berg, G. Valent, The Nevanlinna parametrization for some indeterminate Stieltjes moment problems associated with birth and death processes. Methods and Applications of Analysis 1 (2) (1994), 169–209.
  * [6] H. Buchwalter and G. Cassier, La paramétrisation de Nevanlinna dans le problème des moments de Hamburger, Expo Math. 2 (1984), 155–178.
  * [7] T. S. Chihara, On indeterminate Hamburger moment problems, Pacific J. Math. 27 (1968), 475–484.
  * [8] T. S. Chihara, Indeterminate symmetric moment problems, J. Math. Anal. appl. 85 (1982), 331–346.
  * [9] G. K. Pedersen, Analysis Now, Graduate Texts in Mathematics Vol. 118. Springer Verlag 1989.
  * [10] H. L. Pedersen, The Nevanlinna matrix of entire functions associated with the shifted indeterminate Hamburger moment problem, Math. Scand. 74 (1994), 152–160.
  * [11] H. L. Pedersen, Stieltjes Moment Problems and the Friedrichs Extension of a Positive Definite Operator, J. Approx. Theory 83 (1995), 289–307.
  * [12] H. L. Pedersen, La paramétrisation de Nevanlinna et le problème des moments de Stieltjes indéterminé, Expo Math. 15 (1997), 273–278.
  * [13] A. D. Sokal, J. Walrad, Continued-fraction characterization of Stieltjes moment sequences with support in $[\xi,\infty)$. ArXiv:2404.12131.
  * [14] K. Schmüdgen, The Moment Problem, Graduate Texts in Mathematics Vol. 277. Springer International Publishing AG 2017.

Christian Berg
Department of Mathematical Sciences, University of Copenhagen
Universitetsparken 5, DK-2100 Copenhagen, Denmark
email<EMAIL_ADDRESS>


# Hierarchical Dynamic Routing in Complex Networks via Topologically-decoupled
and Cooperative Reinforcement Learning Agents

Shiyuan Hu<EMAIL_ADDRESS>Network Technology Lab, Huawei Technologies
Co., Ltd., Beijing 100095, China  Shihan Xiao<EMAIL_ADDRESS>Network
Technology Lab, Huawei Technologies Co., Ltd., Beijing 100095, China

###### Abstract

The transport capacity of a communication network can be characterized by the
transition from a free-flow state to a congested state. Here, we propose a
dynamic routing strategy in complex networks based on hierarchical bypass
selections. The routing decisions are made by the reinforcement learning
agents implemented at selected nodes with high betweenness centrality. The
learning processes of the agents are decoupled from each other due to the
degeneracy of their bypasses. Through interactions mediated by the underlying
traffic dynamics, the agents act cooperatively, and coherent actions arise
spontaneously. With only a small number of agents, the transport capacities
are significantly improved, including in real-world Internet networks at the
router level and the autonomous system level. Our strategy is also resilient
to link removals.

##  I Introduction

Efficient and congestion-free transport is critical in many networks, such as
the Internet, the power grid, and airport networks. These networks often show
complex structures with important effects on various dynamic processes [1, 2,
3]. Typically in communication networks, there exists a critical packet
generation rate, beyond which congestion occurs [4, 5, 6, 7, 8], i.e., the
number of buffered packets increases at some intermediate nodes, leading to
performance degradation or even collapse. The onset of congestion defines a
transport capacity for the network, which depends not only on the network
structure but also on the routing strategies.

The widely adopted shortest-path (SP) routing in many systems is to send
packets from source to destination through paths with the fewest number of
links. However, SP routing is inefficient and may lead to severe congestion
due to jamming at nodes with high usage, especially in networks with
heterogeneous degree distributions [9, 10]. By following static paths with the
least sum of node degrees, the transport capacity is significantly improved
[9]. A variety of other elements have also been considered in developing
efficient routing strategies, including local [11, 12, 13, 14] and global [15,
16, 17] traffic conditions, packet priorities [18, 19], and network resource
allocation [20].

With the rapid expansion of scale and explosive growth of traffic volumes in
modern communication networks, traffic control becomes increasingly difficult
and complex. Instead of pushing the limit of transport capacity, we focus on
two other perspectives: dynamic adaptation and distributed control. Dynamic
routing can achieve better load balance than static routing by self-adapting
to diverse network conditions, and distributed control is more feasible and
reactive than centralized control in large-scale systems [21, 22]. As an
unsupervised learning and control approach, reinforcement learning (RL) has
been applied to many domains [23, 24, 25, 26] and also to dynamic routing
problems (e.g., [27] and see reviews in [28, 29]). Recent studies of
distributed dynamic routing with RL mainly focus on hop-by-hop routing [30, 8,
31, 32]. To guarantee fast convergence and packet delivery without loops,
various techniques are needed, such as node identification using some encoding
schemes [30, 32], link reversal [33, 34], and limiting routing changes to only
a few hops (otherwise following SP) [31].

Although the dynamic routing strategies considered in previous studies can
adapt to time-changing traffic conditions, several questions that are both
fundamentally and practically important remain largely unexplored. Instead of
selecting the next hop at each step, how to limit the number of routing
changes and utilize topological properties to bypass congested nodes? If the
routing decisions are distributed and made at each node, how to achieve a
global optimum coherently? Another desired feature is low reliance on the
communications between nodes, since frequent communications may add additional
loads to the network.

In this work, we propose a novel hierarchical dynamic (HD) routing strategy in
heterogeneous networks based on sets of bypasses around selected nodes with
high betweenness centrality (BC). The bypass decisions are made online
according to the real-time traffic conditions by the RL agents implemented at
each selected node. Due to the degeneracy of their action spaces, the agents
are decoupled from each other, avoiding complex competitive multi-agent
settings [35]. Although seemingly independent, different RL nodes interact
with each other through the underlying traffic dynamics, and coherent actions
arise spontaneously across the agents. Our strategy outperforms the least-
degree (LD) routing [9] with larger capacity and lower travel time. Quite
remarkably, even with a small number of agents, the transport capacity is
increased significantly. Different from the routing strategy based on global
traffic conditions [15, 16, 17], the agents can converge and make optimal
decisions based only on the local traffic conditions, which eliminates the
cost of global communications across the network. We further demonstrate that
HD routing is resilient to link removals.

##  II Simulation setup

Figure 1: Schematic of bypass in a BA network with 60 nodes. The nodes are
colored according to their BC. Red line is the SP from node 45 to 58. Green
and blue lines show two bypasses around node 0 computed using Eq. (1) with
$\beta=0.4$ and $1.0$, respectively. 11footnotetext: In Ref. [9], the static
least-degree routing is computed with node degrees by setting
$\beta=1$.22footnotetext: We have also considered to include a bounded history
of the queue lengths and found that its effect on the convergence of the
learning process is negligible.

Traffic model.—Consider the transport of information packets on a complex
network with $N$ nodes. Each node acts as both router and host [7, 9]. The
traffic dynamics evolves forward with discrete time steps. At each time step,
$R$ packets are generated on the network with random sources and destinations.
Each node has a finite buffer (maximum queue length) and can forward at most
one packet following first-in-first-out rule at each step. Additional arriving
packets are dropped when the buffer is full. Delivered packets and dropped
packets are both removed from the network.

We select nodes with high BC to implement the RL agents (hereinafter referred
to as RL nodes), since they are more exposed to traffic under SP routing. The
BC of node $\nu$ is defined as $b(\nu)=\sum_{s\neq
d}\sigma_{s,d}(\nu)/\sigma_{s,d}$ [36], where $\sigma_{s,d}$ is the number of
SPs from node $s$ to $d$ and $\sigma_{s,d}(\nu)$ is how many of them passing
through $\nu$. Denote the set of RL nodes as
$\mathcal{N}^{\alpha}=\\{n_{1}^{\alpha},\cdots,n_{K}^{\alpha}\\}$,
corresponding to the first $K$ ($K\ll N$) nodes with the highest BC. A packet
is transported along one of the SPs after generation. If the next hop is one
of the RL nodes, the packet may take a bypass around that RL node from its
current location to its destination. We compute the bypass as follows.
Consider a SP between source $s$ and destination $d$, $\\{s,\cdots
x,n_{k}^{\alpha},\cdots d\\}$ with $n_{k}^{\alpha}\in\mathcal{N}^{\alpha}$,
and denote the set of all paths between $x$ and $d$ as $\mathcal{P}$, the
bypass around $n_{k}^{\alpha}$ parameterized by $\beta$ is computed as

$p(x\to
d;n_{k}^{\alpha},\beta)=\mathop{\arg\min}_{p\in\mathcal{P}}\left(\sum_{\nu\in
p}b(\nu)^{\beta}\right),$ (1)

where the summation is over nodes along the path $p\in\mathcal{P}$
[[TheBCinEq.(1)maybereplacedwithnodedegreesforfastercomputation,
duetothecorrelationsbetweendegreesandBC, see][]Holme02, 38]. A typical example
of bypass in a Barabási-Albert (BA) network [2] is depicted in Fig. 1.
Equation (1) returns the SP for $\beta=0$ (no bypass). Given different values
of $\beta$, a set of hierarchical bypasses can be obtained: as $\beta$
increases, the bypass becomes longer, but the average BC of nodes on the
bypasses becomes smaller.

Since Eq. (1) only considers topological information, the bypasses between any
pair of nodes only need to form once. If a SP contains no RL node, the SP is
followed without bypass; if a SP contains multiple RL nodes, the one with
higher BC is bypassed, i.e., a packet bypasses at most once.

RL node.—In RL, the intelligent agents interact with the environment and
discover optimal policies through trial and error. At each RL node, we
implement a Q-network [23], a popular RL algorithm that takes continuous state
input and approximates the action-value function corresponding to the optimal
policy. Based on the state of the local traffic condition, the agent performs
an $\epsilon-$greedy action selection at each time step and receives a reward,
the traffic condition then moves to a new state. The experience is then
stored, and the Q-network is updated by replaying the stored experiences.

We use a simple fully connected neural network of three layers with a
rectifier nonlinearity [39]. The received state of a RL node at time $t$
consists of its queue length and the queue lengths of other RL nodes [40]. The
latter is considered since a packet may be directed to another RL node, and
the traffic conditions there may affect the decision-making. The actions are
the bypasses corresponding to a set of $\beta$ values. For all the RL nodes,
we use the same set of $\beta$ values, $\\{0,0.2,0.4,0.6,0.8,1.0,2.0\\}$. The
last ingredient is the reward function, which is computed as
$\text{reward}=-\langle\text{travel time}\rangle-\langle\text{drop
rate}\rangle$, where the first term is the travel time scaled by the maximum
travel time following SP (the product of SP length and the buffer size) and
the second term is the packet drop rate. Both $\langle\cdots\rangle$ denote
average across packets. These two terms represent two competing effects: small
travel time favors shorter paths, while low drop rate favors longer but less-
congested paths with a lower risk of being dropped at severely congested
nodes. Because both statistics are not immediately available after each
action, we divide time into consecutive monitor intervals (MI) [41, 42]. The
RL nodes take actions at the beginning of each MI and perform the same actions
throughout the MI. A RL node only records the packets directed by itself and
computes the reward at the end of the MI.

The training process breaks down into episodes. Each episode restarts the
network traffic and ends after a fixed number of MIs. Throughout the
simulations, the buffer size is 40, $\text{MI}=$ 5–15, and each episode has 50
MIs. The number of neurons in each hidden layer varies from 16–128, depending
on the input size.

Network topologies.—Several practically relevant network topologies are used
in our simulations. We first consider the BA scale-free network, which has a
power-law degree distribution, $\phi(\kappa)\sim\kappa^{-3.0}$. To evaluate
the effect of varying degree heterogeneity, we generate networks using
configuration model [43] with a $J$-tier composite exponential (CE-$J$) degree
distribution, which is given by a sum of exponential distributions with mean
$\lambda^{j+1}$ weighted by $a^{-(j+1)}$ [44]:
$\phi(\kappa)\sim\sum_{j=0}^{J}a^{-(j+1)}\lambda^{j+1}\exp(-\lambda^{j+1}\kappa)$.
When $J\to\infty$, $\phi(\kappa)$ resembles a power law; when $J\to 0$,
$\phi(\kappa)$ degenerates into an exponential distribution. We also consider
two real-world networks, the Internet at AS level [45] and the router-level
topology of an Internet service provider (ISP) [46]. The AS network also has a
power-law degree distribution, $\phi(\kappa)\sim\kappa^{-2.2}$. The basic
topological features of these networks are summarized in Table 1. Both BA and
AS networks have high degree heterogeneity.

| $N$ | $\langle\kappa\rangle$ | $\langle l\rangle$ | $\text{RSD}(\kappa)$ | $H$
---|---|---|---|---|---
BA | $10^{3}$ | 6.0 | 3.5 | 1.2 | 3.5
CE-7 | $10^{3}$ | 5.8 | 3.8 | 0.8 | 2.4
CE-3 | $10^{3}$ | 5.4 | 4.2 | 0.6 | 1.8
ISP | 624 | 17.0 | 3.4 | 1.2 | 2.5
AS | 6474 | 3.9 | 3.7 | 6.4 | 42.4
Table 1: Number of nodes $N$, mean degree $\langle\kappa\rangle$, average
length of shortest paths $\langle l\rangle$, degree relative standard
deviation $\text{RSD}(\kappa)$, and degree heterogeneity,
$H=\langle\kappa^{2}\rangle/\langle\kappa\rangle^{2}$ of three synthetic
networks and two real-world Internet networks used in our simulations. As a
comparison, $\text{RSD}(\kappa)$ and $H$ for relatively homogeneous Erdös-
Rényi network [47] with $N=1000$ and $\langle\kappa\rangle=6.0$ are 0.4 and
1.6, respectively. The statistics of the synthetic networks are averaged from
3 random generations.

##  III Results

To evaluate the effect of nonlocal traffic conditions, we implement 20 RL
nodes in BA network. Figure 2(a) shows that the communications between RL
nodes on their traffic conditions has little effect on the reward, and both
cases converge within 30 episodes, i.e., the learning process of a RL node
only depends on its own traffic condition. This is because the frequency of
redirecting packets to another RL node is around $1/N$ [Fig. 2(a) inset],
which is small in large networks. Similar results are observed in other
networks explored in this study. We further compute the average BC of the
bypasses, $\langle\bar{b}(\beta,n^{\alpha}_{k})\rangle$, where
$\bar{b}(\beta,n^{\alpha}_{k})$ is the average BC of nodes on a bypass of
$\beta$ around a RL node $n^{\alpha}_{k}$ and $\langle\cdots\rangle$ denotes
average obtained by random sampling the source and destination pairs in the
network. Figures 2(b)–(d) demonstrate two folds of hierarchies of bypasses in
BA, ISP, and AS networks. First, for the same RL node
$\langle\bar{b}(\beta,n^{\alpha}_{k})\rangle$ decreases as $\beta$ is
increased and the bypass becomes more marginal. Second, for the same value of
$\beta$, $\langle\bar{b}(\beta,n^{\alpha}_{k})\rangle$ of a node with lower BC
is, in general, degenerated to that of a node with higher BC, i.e., the action
spaces of different RL nodes with different BC are decoupled from each other.
Therefore, we speculate that the inter-node communications on their traffic
conditions may not be important even when $K\sim N$.

Figure 2: The decoupling of the RL nodes. (a) Training reward averaged with a
moving window of 5 episodes for the RL node with the highest BC in BA network
with $K=20$. The green/red curve shows the reward obtained with/without
communications to other RL nodes. The error bars are computed from simulations
on three random network generations. Inset shows the frequency of directing
packets to other RL nodes obtained from the last 30 episodes. (b)–(d) As a
function of $\beta$, $\langle\bar{b}(\beta,n_{k}^{\alpha})\rangle$ of the
first seven nodes with the highest BC (see legend) in (b) BA network, (c) ISP
network, and (d) AS network. We normalize BC by $2/\left[(N-1)(N-2)\right]$.

We show the frequency distributions of actions (different values of $\beta$)
in Fig. 3. During the traffic simulation, even after convergence, the queue
length at a RL node is time-changing, and different values of $\beta$ may be
selected. However, the most frequent actions are different at different packet
generation rates. At small $R/N$, packets travel across the network with
little queuing, and the bypasses of small values of $\beta$ are frequently
selected to minimize the travel time (blue lines in Fig. 3). As $R/N$
increases, the RL nodes become more congested, and the frequencies of larger
values of $\beta$ increase (orange and red lines). Arriving packets that need
to wait in long queues or exceed the finite buffers are gradually
redistributed to nodes with lower BC.

Figure 3: The frequency distributions of actions $P(\beta)$ for three RL nodes
(corresponding to three columns with decreasing BC from left to right) with
different $R/N$ in (a) BA network, (b) CE-3 network, (c) ISP network, and (d)
AS network. Results are obtained from the last 30 episodes. Figure 4: The
order parameter $\eta$ as a function of packet generation rate $R/N$ in semi-
log scale for an increasing $K$ (from light red to dark red circles) in five
different networks. The gray circles show data for SP routing and the blue
circles show data for LD routing. For (a)-(d), $K=0$–15. The results on BA,
CE-7, and CE-3 networks are averaged from independent simulations on three
random network generations. The transport capacity $R_{c}$ is estimated as the
intersection of the linear trend and $\eta=0$ (see green line as an example).
Insets show the relative increase of $R_{c}$ as a function of $K$ compared
with SP routing.

Another important observation from Fig. 3 is that, for different RL nodes in
the same network, their frequency distributions are quite similar with the
most frequent values of $\beta$ close to each other, even though the actions
are taken independently based on their own traffic conditions. The coherent
action selections indicate that different RL nodes act cooperatively and are
non-competitive over network resources, since the bypasses of the same $\beta$
for different RL nodes are degenerated from each other [Figs. 2(b)–(c)]. As an
example, in BA network with $R/N=0.06$, the most frequency action for the RL
node with the highest BC [Fig. 3(a) left] is $\beta=0.6$. If the second RL
node [Fig. 3(a) middle] selects smaller values of $\beta$ more frequently, it
would be discouraged by long queuing and even packet loss, since its bypasses
of $\beta<0.6$ have larger $\langle\bar{b}(\beta,n^{\alpha}_{k})\rangle$ than
the bypasses of $\beta=0.6$ of the first RL node [Fig. 2(b)] and are already
occupied by the packets directed from the first RL node. Actions of larger
values of $\beta$ are also not optimal due to long bypass length. As a result,
the most frequent action at $R/N=0.06$ for the second RL node is also around
$\beta=0.6$. This hierarchical exploitation of network resources continues as
we add more RL nodes. Therefore, the coherent actions depicted in Fig. 3 arise
spontaneously without communications with each other, and the global optimum
is achieved after each RL node behaves optimally.

To characterize the transition from free-flow state to congested state, we
compute the order parameter [6],

$\eta=\frac{1}{R}\frac{\langle W(t+\Delta t)-W(t)\rangle}{\Delta t}.$ (2)

Here, $\langle\cdots\rangle$ indicates a moving average over a time window of
size $\Delta t$, and $W(t)$ is the sum of the number of in-transit and
cumulative dropped packets at time $t$. At small $R/N$, packets are delivered
with little queuing on the network and $\eta=0$. The transport capacity is
characterized by a critical generation rate $R_{c}/N$, beyond which packets
start to accumulate in queues and $\eta>0$. Figure 4 shows $\eta$ for five
different networks with different number of RL nodes $K$. The increase of
$R_{c}$ relative to SP routing is computed as, $\Delta
R_{c}(K)=[R_{c}(K)-R_{c}(0)]/R_{c}(0)$ (Fig. 4 insets). We observe significant
increases in $R_{c}$, especially in networks with high degree heterogeneity
(BA and AS networks).

Figure 5: (a)–(c) The frequency distributions of packet travel time $P(T)$ in
semi-log scale for SP (gray circles), HD (red circles) with $K=5$, and LD
routing (blue circles) at three different packet generation rates in BA
network, (a) $R/N=0.002$, (b) $R/N=0.02$, and (c) $R/N=0.05$. (d) The
distribution of cumulative packet loss on BC at $t=500$.

In BA network, $R_{c}$ is increased by more than 10 times with $K/N\sim 0.01$;
in AS network, $R_{c}$ is increased by around 6 times with $K/N\sim 0.003$. In
these networks, as $K$ is increased, $R_{c}$ increases sharply first and then
gradually saturates, since the routings of the majority of packets are
controlled by a few RL nodes due to the fat-tailed $\phi(\kappa)$. As
$\phi(\kappa)$ becomes more homogeneous from BA to CE-7, and then to CE-3, the
increase of $R_{c}$ slows down. In CE-3 network, with $\langle H\rangle$ close
to that of Erdös-Rényi network [Table 1], $R_{c}$ is still doubled with only
1% RL nodes. Our strategy outperforms LD routing in BA network [Fig. 4(a) blue
line]. This is because in LD routing packets are distributed to nodes with
small degrees and travel long distances, regardless of the congestion level.
But in our strategy, a set of hierarchical bypasses of different lengths are
dynamically selected depending on the traffic conditions [Fig. 3].

To further demonstrate that our strategy achieves a balance between routing
through hub nodes at the network center and nodes at the network periphery, we
show in Fig. 5 the statistics of packet travel time and dropped packets. At
small $R/N$, the frequency distribution of packet travel time of HD routing is
close to that of SP routing, and both are much smaller than LD routing [Fig.
5(a)]. As $R/N$ increases, the distribution of HD routing gradually shifts
towards LD routing as more and more packets take longer bypasses. The peaks in
SP routing in Figs. 5(b) and 5(c) around $T=40$ and 80 are due to overfilled
buffers: packets at the end of a queue have to wait for the whole queue to be
forwarded. This also indicates severe packet loss. Figure 5(d) shows the
distribution of cumulative packet loss on BC. As two extremes, there are
considerable packet losses at nodes with high BC in SP routing and nodes with
low BC in LD routing. The introduction of 3 RL nodes first significantly
reduces packet losses by about two orders of magnitude at nodes with high BC
[green circles in Figure 5(d)]. As $K$ is increased from 3 to 9, the packet
losses at nodes with intermediate BC around 0.05–0.1 are also reduced since
more packets are redistributed to nodes with low BC.

Figure 6: Effect of link removals with $R/N=0.001$ and $K=3$ in BA network.
(a) Training reward averaged with a moving window of 3 episodes for no removal
(gray), random removal (green), and BC removal (red). (b) The frequency
distributions of $\beta$ for the RL node with the highest BC. The error bars
are computed from independent simulations on three random network generations.

Finally, we evaluate the effect of link removals. Two removal strategies are
considered, random removal and removal with probabilities proportional to the
average BC of the node pairs connecting to the links (BC removal). The former
may be caused by attacks with little information about the network; the latter
may be caused by attacks with partial information. Different from the previous
approaches to increase the transport capacity by deliberate link removal or
addition [48, 49], neither the RL nodes nor the packets are aware of the
topological changes: the paths are not updated on the network after link
removal. If the link between a packet’s current node and its next-hop node is
missing, the packet is dropped. We remove 2% of the links at the end of
episode 20. In the free-flow state ($R/N=0.001$), Fig. 6(a) shows that both
removal strategies reduce the reward immediately caused by packet losses at
missing links. Surprisingly, the recovery from BC removal is faster than
random removal. This is due to the build-in hierarchies of the bypasses. Since
at small $R/N$ most packets follow bypasses of small values of $\beta$, more
penalty is assigned to these actions than actions of large values of $\beta$.
Indeed, the frequencies of $\beta=0$ decrease, and the frequencies of larger
values of $\beta$ increase for both removal strategies [Fig. 6(b)]. However,
for the BC removal, links in bypasses of small values of $\beta$ are removed
with larger probabilities. Therefore, actions of small values of $\beta$ are
penalized more compared with random removal, providing clearer feedback to the
RL nodes that bypasses of small values of $\beta$ should be avoided. Indeed,
the frequency distribution shifts towards larger $\beta$ values for the BC
removal.

##  IV Discussion

We have studied a hierarchical dynamic routing strategy in heterogeneous
complex networks with the routing decision-making distributed at
topologically-decoupled RL nodes. The cooperative behaviors of the RL nodes
and their coherent actions do not require explicit coordination through inter-
node communications. It arises from the degeneracy of their actions spaces and
indirect interactions mediated by the traffic dynamics. Most importantly, our
results suggest that the transport capacity can be significantly increased by
implementing only a small number of RL nodes, much smaller than the total
number of nodes of the network. Our results may also be useful for the design
of distributed intelligent agents in complex systems and provide insights into
the understanding of their collective behaviors.

A balance of traffic between hub nodes and peripheral nodes has also been
realized by computing an indicator of the traffic conditions at hub nodes, but
the queue lengths at all nodes are required [17]. Unlike hop-by-hop routing,
our strategy is based on selections of bypasses between source and
destination, and therefore loop-free routing is guaranteed. The conventional
congestion control mechanism on the Internet along fixed paths typically only
involves the source and destination nodes with the network structure
functioning only to transport the packets. Based on the feedback information
from the destination, the source can decrease its sending rate when congestion
has likely occurred [50]. In contrast to bypassing congested nodes
dynamically, this end-to-end approach often limits network utilization.

Our strategy may be further improved from several perspectives. For large-
population agents in large-scale networks, we can speed up training by sharing
the parameters of a single Q-network among all RL nodes [51, 30]. The RL nodes
can still make different decisions with their node identifications as
additional inputs to the Q-network. We may also assign different weights to
the two terms in the reward function and tune the preference over low travel
time or low drop rate. In this study, we have only considered discrete action
spaces (discrete $\beta$ values), continuous actions can be achieved with
policy gradient algorithm [52], but may require large networks to
differentiate bypasses with small difference in $\beta$.

###### Acknowledgements.

We thank Xiaofei Xu for fruitful discussions and Gang Yan for insightful
comments and suggestions.

## References

  * Watts and Strogatz [1998] D. J. Watts and S. H. Strogatz, Nature 393, 440 (1998).
  * Barabási and Albert [1999] A.-L. Barabási and R. Albert, Science 286, 509 (1999).
  * Boccaletti  _et al._ [2006] S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, and D.-U. Hwang, Phys. Rep. 424, 175 (2006).
  * Takayasu  _et al._ [1996] M. Takayasu, H. Takayasu, and T. Sato, Phys. A: Stat. Mech. Appl. 233, 824 (1996).
  * Ohira and Sawatari [1998] T. Ohira and R. Sawatari, Phys. Rev. E 58, 193 (1998).
  * Arenas  _et al._ [2001] A. Arenas, A. Díaz-Guilera, and R. Guimerà, Phys. Rev. Lett. 86, 3196 (2001).
  * Zhao  _et al._ [2005] L. Zhao, Y.-C. Lai, K. Park, and N. Ye, Phys. Rev. E 71, 026125 (2005).
  * Ding  _et al._ [2020] R. Ding, Y. Yang, J. Liu, H. Li, and F. Gao, in  _2020 International Conference on Computing, Networking and Communications (ICNC)_ (2020) pp. 932–937.
  * Yan  _et al._ [2006] G. Yan, T. Zhou, B. Hu, Z.-Q. Fu, and B.-H. Wang, Phys. Rev. E 73, 046108 (2006).
  * De Martino  _et al._ [2009] D. De Martino, L. Dall’Asta, G. Bianconi, and M. Marsili, Phys. Rev. E 79, 015101 (2009).
  * Echenique  _et al._ [2004] P. Echenique, J. Gómez-Gardeñes, and Y. Moreno, Phys. Rev. E 70, 056105 (2004).
  * Echenique  _et al._ [2005] P. Echenique, J. Gómez-Gardeñes, and Y. Moreno, EPL 71, 325 (2005).
  * Wang  _et al._ [2006] W.-X. Wang, C.-Y. Yin, G. Yan, and B.-H. Wang, Phys. Rev. E 74, 016101 (2006).
  * Tang and Zhou [2011] M. Tang and T. Zhou, Phys. Rev. E 84, 026116 (2011).
  * Zhang  _et al._ [2007a] H. Zhang, Z. Liu, M. Tang, and P. Hui, Phys. Lett. A 364, 177 (2007a).
  * Ling  _et al._ [2010] X. Ling, M.-B. Hu, R. Jiang, and Q.-S. Wu, Phys. Rev. E 81, 016113 (2010).
  * Tan and Xia [2013] F. Tan and Y. Xia, Phys. A: Stat. Mech. Appl. 392, 4146 (2013).
  * Kim  _et al._ [2009] K. Kim, B. Kahng, and D. Kim, EPL 86, 58002 (2009).
  * Du  _et al._ [2013] W.-B. Du, Z.-X. Wu, and K.-Q. Cai, Phys. A: Stat. Mech. Appl. 392, 3505 (2013).
  * Wu  _et al._ [2013] J. Wu, C. K. Tse, F. C. M. Lau, and I. W. H. Ho, IEEE Trans. Circuits Syst. I Regul. Pap. 60, 3303 (2013).
  * Antonelli [2013] G. Antonelli, IEEE Control Syst. Mag. 33, 76 (2013).
  * Marden and Shamma [2015] J. R. Marden and J. S. Shamma, Game theory and distributed control, in  _Handbook of Game Theory with Economic Applications_, Vol. 4 (Elsevier, 2015) pp. 861–899.
  * Mnih  _et al._ [2015] V. Mnih _et al._ , Nature 518, 529 (2015).
  * Boutaba  _et al._ [2018] R. Boutaba _et al._ , J. Internet Serv. Appl. 9, 16 (2018).
  * Garnier  _et al._ [2021] P. Garnier, J. Viquerat, J. Rabault, A. Larcher, A. Kuhnle, and E. Hachem, Comput. Fluids 225, 104973 (2021).
  * Degrave  _et al._ [2022] J. Degrave _et al._ , Nature 602, 414 (2022).
  * Boyan and Littman [1993] J. Boyan and M. Littman, in  _Advances in Neural Information Processing Systems_, Vol. 6 (Morgan-Kaufmann, 1993).
  * Al-Rawi  _et al._ [2015] H. A. A. Al-Rawi, M. A. Ng, and K.-L. A. Yau, Artif. Intell. Rev. 43, 381 (2015).
  * Dai  _et al._ [2021] B. Dai, Y. Cao, Z. Wu, Z. Dai, R. Yao, and Y. Xu, Neurocomputing 459, 44 (2021).
  * Mukhutdinov  _et al._ [2019] D. Mukhutdinov, A. Filchenkov, A. Shalyto, and V. Vyatkin, Future Gener. Comput. Syst. 94, 587 (2019).
  * Kang  _et al._ [2021] Y. Kang, X. Wang, and Z. Lan, in  _Proceedings of the 30th International Symposium on High-Performance Parallel and Distributed Computing_ (2021) pp. 189–200.
  * You  _et al._ [2022] X. You, X. Li, Y. Xu, H. Feng, J. Zhao, and H. Yan, IEEE Trans. Syst. Man Cybern.: Syst. 52, 855 (2022).
  * Gafni and Bertsekas [1981] E. Gafni and D. Bertsekas, IEEE Trans. Commun. 29, 11 (1981).
  * Xiao  _et al._ [2020] S. Xiao, H. Mao, B. Wu, W. Liu, and F. Li, in  _Proceedings of the Workshop on Network Meets AI & ML_ (2020) pp. 28–34.
  * Zhang  _et al._ [2021] K. Zhang, Z. Yang, and T. Başar, Multi-agent reinforcement learning: A selective overview of theories and algorithms, in  _Handbook of Reinforcement Learning and Control_ (Springer International Publishing, Cham, 2021) pp. 321–384.
  * Freeman [1977] L. C. Freeman, Sociometry 40, 35 (1977).
  * Holme  _et al._ [2002] P. Holme, B. J. Kim, C. N. Yoon, and S. K. Han, Phys. Rev. E 65, 056109 (2002).
  * Note [1] In Ref. [9], the static least-degree routing is computed with node degrees by setting $\beta=1$.
  * Nair and Hinton [2010] V. Nair and G. E. Hinton, in  _ICML_ (2010) pp. 807–814.
  * Note [2] We have also considered to include a bounded history of the queue lengths and found that its effect on the convergence of the learning process is negligible.
  * Dong  _et al._ [2015] M. Dong, Q. Li, D. Zarchy, P. B. Godfrey, and M. Schapira, in  _12th USENIX Symposium on Networked Systems Design and Implementation (NSDI 15)_ (2015) pp. 395–408.
  * Jay  _et al._ [2019] N. Jay, N. Rotman, B. Godfrey, M. Schapira, and A. Tamar, in  _Proceedings of the 36th International Conference on Machine Learning_, Vol. 97 (2019) pp. 3050–3059.
  * Newman [2003] M. E. J. Newman, SIAM Rev. 45, 167 (2003).
  * Reynolds [2014] A. M. Reynolds, Sci. Rep. 4, 4409 (2014).
  * Leskovec and Krevl [2014] J. Leskovec and A. Krevl, SNAP Datasets: Stanford large network dataset collection, http://snap.stanford.edu/data (2014).
  * Spring  _et al._ [2002] N. Spring, R. Mahajan, and D. Wetherall, SIGCOMM Comput. Commun. Rev. 32, 133 (2002).
  * Erdös and Rényi [1960] P. Erdös and A. Rényi, Publ. Math. Inst. Hung. Acad. Sci. 5, 17 (1960).
  * Zhang  _et al._ [2007b] G.-Q. Zhang, D. Wang, and G.-J. Li, Phys. Rev. E 76, 017101 (2007b).
  * Jiang  _et al._ [2011] Z. Jiang, M. Liang, and D. Guo, Int. J. Mod. Phys. 22, 1211 (2011).
  * Jacobson [1988] V. Jacobson, SIGCOMM Comput. Commun. Rev. 18, 314 (1988).
  * Zheng  _et al._ [2018] L. Zheng, J. Yang, H. Cai, M. Zhou, W. Zhang, J. Wang, and Y. Yu, in  _Proceedings of the AAAI Conference on Artificial Intelligence_, Vol. 32 (2018).
  * Lillicrap  _et al._ [2015] T. P. Lillicrap _et al._ , arXiv preprint arXiv:1509.02971 (2015).