# Markov chains and mappings of distributions on compact spaces

David J. Aldous Department of Statistics, 367 Evans Hall # 3860, U.C. Berkeley
CA 94720<EMAIL_ADDRESS>www.stat.berkeley.edu/users/aldous. Shi
Feng Department of Mathematics, Cornell University<EMAIL_ADDRESS>
###### Abstract

Consider a compact metric space $S$ and a pair $(j,k)$ with $k\geq 2$ and
$1\leq j\leq k$. For any probability distribution $\theta\in\mathcal{P}(S)$,
define a Markov chain on $S$ by: from state $s$, take $k$ i.i.d. ($\theta$)
samples, and jump to the $j$’th closest. Such a chain converges in
distribution to a unique stationary distribution, say $\pi_{j,k}(\theta)$. So
this defines a mapping $\pi_{j,k}:\mathcal{P}(S)\to\mathcal{P}(S)$. What
happens when we iterate this mapping? In particular, what are the fixed points
of this mapping? We present a few rigorous results, to complement our
extensive simulation study elsewhere.

Keywords: Coupling, Markov chain, compact metric space, dynamical system.

AMS MSC 2020: 60J05, 37A50.

##  1 Introduction

This article discusses a rather novel topic whose motivation may seem obscure,
so we start with informal background that led to the formulation of the topic.
Write $S=(S,d)$ for a compact metric space.. Then the identity function
$f(s):=s$ makes sense for every $S$. Is there any more interesting explicit
function $S\to S$ whose definition makes sense for every $S$? For example one
might try $f(s):=\arg\max_{y}d(s,y)$, that is the most distant point from $s$;
this works for any space $S$ with the property that the most distant point is
always unique, but not for all $S$. Our introspection suggests that in fact
there is no non-trivial such “general” function.

Instead let us write $\mathcal{P}(S)$ for the space of probability
distributions on $S$, and recall that $\mathcal{P}(S)$ is a compact metric
space under the usual weak topology. The observation above suggests that there
may be no non-trivial function $\mathcal{P}(S)\to\mathcal{P}(S)$ whose
definition makes sense for every $S$. But this is false! This article
investigates a particular family of such functions – the reader may care to
try to invent different examples.

#### The Markov chain.

Given $S$ and $\theta\in\mathcal{P}(S)$, consider the following discrete time
Markov chain on state space $S$: from point $s$ make the step to the nearer of
$2$ random points drawn i.i.d. from $\theta$, breaking ties uniformly at
random. This scheme naturally generalizes as follows: fix $k\geq 2$ and $1\leq
j\leq k$, and step from $s$ to the $j$’th nearest of $k$ random points drawn
i.i.d. from $\theta$, again breaking ties uniformly at random. This defines a
kernel $K^{\theta,j,k}$ on $S$. Write the associated chain as
$\mathbf{X}^{\theta,j,k}=(X^{\theta,j,k}(t),t=0,1,2,\ldots)$. Theorem 1 proves
that this chain always has a unique stationary distribution, which we can call
$\pi_{j,k}(\theta)$. So now we have defined a mapping
$\pi_{j,k}:\mathcal{P}(S)\to\mathcal{P}(S)$ for every $S$. Theorem 1 also
proves that the distributions $\theta$ and $\pi_{j,k}(\theta)$ are mutually
absolutely continuous, so in particular have the same support.

#### Invariant measures for the mappings.

These mappings $\pi_{j,k}$ have apparently not been studied previously, even
for special spaces $S$ and the simplest case $k=2$. Amongst the range of
questions one could ask, we will seek to study it as a dynamical system. Given
a mapping $\pi$ from a space to itself, it is mathematically natural to
consider iterates

$\pi^{n+1}(\theta)=\pi(\pi^{n}(\theta)),n\geq 1.$ (1)

In our setting it seems plausible that (at least for typical initial $\theta$)
the iterates should converge to some limit, that is we expect weak convergence

$\pi^{n}_{j,k}(\theta)\to_{w}\phi\mbox{ as }n\to\infty$ (2)

and then we expect111As observed in [1], the Markov chain is not always a
Feller process, so (3) does not immediately follow from (2). the limit $\phi$
to satisfy the fixed point or invariant distribution condition

$\pi_{j,k}(\phi)=\phi.$ (3)

Some comments about this set-up.

(a) The iterative procedure (1) does not have any simple stochastic process
interpretation, in contrast to the mapping $\theta\to\pi_{j,k}(\theta)$
derived from the Markov chain.

(b) The equation for the stationary distribution $\pi_{j,k}(\theta)$, which
for finite $S$ is the elementary Markov chain relation $\pi=\pi\mathbf{P}$, is
a linear equation, whereas the fixed point equation (3) is decidedly non-
linear.

(c) On any $S$ and for any $(j,k)$, two types of measures are always
invariant: we call these the omnipresent measures.

  * •

The distribution $\delta_{s}$ degenerate at one point $s$;

  * •

The uniform two-point distribution
$\delta_{s_{1},s_{2}}=\frac{1}{2}(\delta_{s_{1}}+\delta_{s_{2}})$.

(d) If an invariant distribution on $S$ has support $S_{0}\subset S$ then we
can regard it as an invariant distribution on $S_{0}$. So the essential
question is: given $S$, what are the invariant distributions with full
support? Note that when $\pi^{n}_{j,k}(\theta)\to_{w}\phi$ the distributions
$(\pi^{n}_{j,k}(\theta),n\geq 1)$ all have the same support (by Theorem 1) but
$\phi$ may (as usually has, it turns out) have smaller support.

(e) If $\theta\in\mathcal{P}(S)$ is invariant under an isometry $\iota$ of $S$
then $\pi_{j,k}(\theta)$ is also invariant under $\iota$. So we expect a limit
of the iterative procedure (1) to also be invariant under $\iota$. We say
these invariant measures “exist by symmetry”, for instance the Haar
probability measure on a compact group S with a metric invariant under the
group action.

#### Motivation.

There is no notion of “uniform distribution” applicable to every compact
metric space $S$. The original motivation for this project was the hope that
our invariant distributions might provide a proxy for uniform distributions on
a general $S$. We attempted to find such distributions via numerically
implementing the iterative procedure on various spaces $S$. What we found was
that, in the absence of some special symmetry property preserved under the
iterative procedure, one almost always obtained a limit supported on only one
or two points, the omnipresent measures mentioned above. This seemed counter-
intuitive, and prompted the further study of invariant measures, even though
the original motivation turned out to be unsuccessful.

.

#### What numerics and simulation suggest.

Our quite extensive study via numerics and simulation is described in a
companion document [1], and suggests the following big picture.

(a) For $k=2$, there are no invariant measures other than the omnipresent
ones, except perhaps for “exist by symmetry” ones; with that exception, for
$j=1,k=2$ the iterates (1) converge to some $\delta_{s}$, and for $j=2,k=2$
the iterates (1) converge to some $\delta_{s_{1},s_{2}}$. The precise limits
$(s,s_{1},s_{2})$ may depend on the initial $\theta$. In the case of
$\delta_{s_{1},s_{2}}$, the pair $(s_{1},s_{2})$ is a local maximum of
$d(\cdot,\cdot)$.

(b) For larger $k$, for some types of space $S$ there are additional sporadic
invariant measures; we don’t see a pattern.

(c) For large $k$, as $j$ increases we see a transition, around $j/k=0.7$,
between convergence to some $\delta_{s}$ and convergence to some
$\delta_{s_{1},s_{2}}$. However there seems no reason to believe that there is
a universal value near $0.7$.

(d) Except for the omnipresent ones, all invariant measures $\phi$ that we
have encountered are unstable, in that from any initial distribution that is
$\phi$ plus a generic (not symmetry-preserving) small perturbation, the
iterates converge to some $\delta_{s}$ or $\delta_{s_{1},s_{2}}$.

#### What can we actually prove?

In short: very little. Here are the results that we will derive in this
article.

  * •

Theorem 1 is the Markov chain convergence result.

  * •

Results in section 3 for $|S|=2$ or $3$ are consistent with general picture
above.

  * •

Theorem 3: For every $S$, the set of invariant distributions for $\pi_{1,2}$
is the same as the set of invariant distributions for $\pi_{2,2}$. This is
surprising, in that apparently (as in (a) above) the iterates almost always
converge to some $\delta_{s}$ for $\pi_{1,2}$, but to some
$\delta_{s_{1},s_{2}}$ for $\pi_{2,2}$.

  * •

Theorem 4: There are no $\pi_{1,2}$ or $\pi_{2,2}$-invariant distributions on
the interval $[0,1]$ other than the omnipresent ones.

  * •

Theorem 7: There are no $\pi_{1,2}$ or $\pi_{2,2}$-invariant distributions on
a space of finite binary tree leaves other than the omnipresent ones.

Of course, for any specific $S$, one can simply write out the fixed point
definition (3) and seek some ad hoc method of finding all solutions. The
results above carry this through (for $\pi_{1,2}$) for $|S|=3$ and for the
interval $[0,1]$, and for leaf-labeled binary trees. But these are essentially
“proofs by contradiction” using specific features of the specific class of
spaces. For general $S$ and $\pi_{1,2}$ one feels there should be some
“contraction” argument for the iterates $\pi^{n}_{1,2}(\theta)$ – the
distributions should become more concentrated as $n$ increases – but we are
unable to formalize that general idea.

##  2 Existence and uniqueness of stationary distributions

###### Theorem 1

Consider a compact metric space $(S,d)$ and a probability distribution
$\theta\in\mathcal{P}(S)$. For each pair $1\leq j\leq k,\ k\geq 2$, the Markov
chain $\mathbf{X}^{\theta,j,k}=(X^{\theta,j,k}(t),t=0,1,2,\ldots)$ has a
unique stationary distribution $\pi_{j,k}(\theta)$. From any initial point,
the variation distance $D(t)$ between $\pi_{j,k}(\theta)$ and the distribution
of $X^{\theta,j,k}(t)$ satisfies

$D(2t)\leq(1-1/k^{k-1})^{t},\quad 1\leq t<\infty$ (4)

and so there is convergence to stationarity in variation distance. Moreover,
for $\pi=\pi_{j,k}(\theta)$

$\theta^{k}(A)\leq\pi(A)\leq k\theta(A),\ A\subseteq S$ (5)

and so $\pi$ and $\theta$ are mutually absolutely continuous.

Note that the bound on variation distance depends only on $k$.

Proof.  First note that for any partition $(B_{i},1\leq i\leq k)$ of $S$ we
have

$\sum_{i}(\theta(B_{i}))^{k}\geq 1/k^{k-1}$ (6)

because by convexity the sum is minimized when $\theta(B_{i})\equiv 1/k$.

We construct the process $X(t)=X^{\theta,j,k}(t)$ in the natural way, by
creating i.i.d. $\theta$-distributed $(\mathbf{Y}(t)=(Y_{i}(t),1\leq i\leq
k),t\geq 1)$ and defining for $t\geq 1$

> $X(t)$ is the element of $(Y_{i}(t),1\leq i\leq k)$ attaining the $j$’th
> smallest value of $(d(X(t-1),Y_{i}(t)),1\leq i\leq k)$.

In defining the re-ordering to determine “$j$’th smallest”, we break ties in
accordance with the original $i$ – that is, if
$d(X(t-1),Y_{i_{1}}(t))=d(X(t-1),Y_{i_{2}}(t))$ for $i_{1}<i_{2}$ then we put
the $i_{1}$ term before the $i_{2}$ term in the reordering. Because the
$Y_{i}$ are i.i.d. this has the same effect as breaking the tie randomly.

We define the natural coupling $(X(t),X^{\prime}(t))$ of two chains started
from arbitrary different states by using the same realizations of $Y_{i}(t)$
for each chain. We first seek to upper bound the coupling time
$T:=\min\\{t:X(t)=X^{\prime}(t)\\}$. Consider a realization
$\mathbf{y}=(y_{i},1\leq i\leq k)$ of $\mathbf{Y}(t+2)$. This $\mathbf{y}$
induces a partition of $S$, say $(B_{i}(\mathbf{y}),1\leq i\leq k)$, where
$B_{i}(\mathbf{y})$ is the set of $s\in S$ such that $d(s,y_{i})$ is the
$j$’th smallest of $(d(s,y_{u}),1\leq u\leq k)$, breaking ties as above. The
central part of the proof is the observation that the event $\\{T\leq t+2\\}$
includes the event

each component of $\mathbf{Y}(t+1)$ is in the same set of the partition
$(B_{i}(\mathbf{Y}(t+2)),1\leq i\leq k)$. (7)

Now $\mathbf{Y}(t+1)$ is independent of $\mathbf{Y}(t+2)$, so we can apply (6)
to show that event (7) has probability $\geq 1/k^{k-1}$. This remains true
conditional on $(X(t),X^{\prime}(t))$, and hence conditional on $\\{X(t)\neq
X^{\prime}(t)\\}$, implying that

${\mathbb{P}}(T\leq t+2|T>t)\geq 1/k^{k-1}.$

So inductively

${\mathbb{P}}(T>2t)\leq(1-1/k^{k-1})^{t},\quad 1\leq t<\infty.$ (8)

This is true for arbitrary initial distributions $\theta$ and
$\theta^{\prime}\in\mathcal{P}(S)$, and so in particular for $\theta$ and
$\theta^{(2)}$, where $\theta^{(t)}$ denotes the distribution of
$X^{\theta,j,k}(t)$. So (8) bounds the variation distance

$||\theta^{(2t+2)}-\theta^{(2t)}||_{VD}\leq(1-1/k^{k-1})^{t},\quad 1\leq
t<\infty$

and similarly

$||\theta^{(2t+1)}-\theta^{(2t)}||_{VD}\leq(1-1/k^{k-1})^{t},\quad 1\leq
t<\infty.$

Now variation distance is a complete metric on $\mathcal{P}(S)$, so
$\theta^{(t)}$ converges in variation distance to a limit $\pi$, and $\pi$ is
a stationary distribution for the kernel $K^{\theta,j,k}$. Then applying (8)
to $\pi$ and an arbitrary other initial distribution establishes (4) and shows
that $\pi$ is the unique stationary distribution. Then (5) follows by
considering the first step $(X(0),X(1))$ of the stationary chain, because for
$A\subset S$

$\cap_{i}\\{Y_{i}(1)\in A\\}\subseteq\\{X(1)\in
A\\}\subseteq\cup_{i}\\{Y_{i}(1)\in A\\}.$

Remarks. The variation distance bound (4) is exponentially decreasing in time,
but it is more natural to consider mixing time in the sense of [2]. The
example of the uniform distribution $\theta$ on a 2-point space with $j=1$
shows that the mixing time as a function of $k$ can be order $2^{k}$.

The proof of Theorem 1 does not say anything about $\pi_{j,k}(\theta)$ except
(5). We do not know if there are informative analytic descriptions of
$\pi_{j,k}(\theta)$ in terms of $\theta$.

##  3 Two or 3 points

###  3.1 Two points – the binomial case

The case of a 2-element space $S=\\{a,b\\}$ and general $(j,k)$ is not
completely trivial. Here is an outline – for more details see [1].

Parametrizing a distribution $\theta$ on $S$ by $p:=\theta(a)$, we view the
mapping $\pi_{j,k}:\mathcal{P}(S)\to\mathcal{P}(S)$ as a mapping
$\pi_{j,k}:[0,1]\to[0,1]$. From the stationary distribution we find, in terms
of binomial variables,

$\pi_{j,k}(p)=\frac{{\mathbb{P}}(\mathrm{Bin}(k,p)>k-j)}{{\mathbb{P}}(\mathrm{Bin}(k,p)>k-j)+{\mathbb{P}}(\mathrm{Bin}(k,p)<j)}.$

So a fixed point is a solution of the equation

$\pi_{j,k}(p)=p.$ (9)

The omnipresent fixed points are $p=0,p=1/2,p=1$; are there others? By
symmetry it is enough to consider $0<p<1/2$.

For given $(j,k)$, we observe in [1] three possible types of qualitative
behavior:

(i) $\pi^{n}_{j,k}(p)\to 0\mbox{ as }n\to\infty,\mbox{ for all }0<p<1/2$.

(ii) $\pi^{n}_{j,k}(p)\to 1/2\mbox{ as }n\to\infty,\mbox{ for all }0<p<1/2$.

(iii) There exists a critical value $p_{crit}\in(0,1/2)$ which is unstable:
that is,

$p_{crit}$ is invariant

and $\pi^{n}_{j,k}(p)\to 0\mbox{ as }n\to\infty,\mbox{ for all }0<p<p_{crit}$

and $\pi^{n}_{j,k}(p)\to 1/2\mbox{ as }n\to\infty,\mbox{ for all
}p_{crit}<p<1/2$.

Case (iii) first arises with $k=5,j=4$, and persists for larger values of $k$.
For instance, with $k=8$ we observe case (i) for $1\leq j\leq 5$, case (iii)
for $j=6$ with $p_{crit}=0.26405$, and case (ii) for $j=7,8$.

Of course the 2-point space is very special. The occurrence of these
“sporadic” case (iii) fixed points seems much rarer in other spaces.

###  3.2 Three elements

Here we consider $S=\\{a,b,c\\}$ where the three distances are distinct, say

$d(a,b)<d(a,c)<d(b,c).$ (10)

###### Theorem 2

If $S$ satisfies (10) then there is no $\pi_{1,2}$-invariant distribution
except the omnipresent ones.

Proof.  It is enough to prove that there is no invariant distribution
$\theta=(\theta_{a},\theta_{b},\theta_{c})$ with each term strictly positive.
So suppose, to get a contradiction, such $\theta$ exists.

Take $Y,Y_{1},Y_{2}$ independent with distribution $\theta$. Invariance says
that the random variable $X$ defined as

$\displaystyle X$ $\displaystyle=$ $\displaystyle Y_{1}\mbox{ if
}d(Y,Y_{1})<d(Y,Y_{2})$ $\displaystyle=$ $\displaystyle Y_{2}\mbox{ if
}d(Y,Y_{2})<d(Y,Y_{1})$

will also have distribution $\theta$. Writing out the ways that $X$ can be $c$
or $b$ or $a$ gives the equations

$\displaystyle\theta_{c}$ $\displaystyle=$
$\displaystyle\theta_{c}(1-(1-\theta_{c})^{2})+\theta_{b}\theta_{c}^{2}+\theta_{a}\theta_{c}^{2}$
$\displaystyle\theta_{b}$ $\displaystyle=$
$\displaystyle\theta_{c}\theta_{b}^{2}+\theta_{b}(1-(1-\theta_{b})^{2})+\theta_{a}(\theta_{b}^{2}+2\theta_{b}\theta_{c})$
$\displaystyle\theta_{a}$ $\displaystyle=$
$\displaystyle\theta_{c}(\theta_{a}^{2}+2\theta_{a}\theta_{b})+\theta_{b}(\theta_{a}^{2}+2\theta_{a}\theta_{c})+\theta_{a}(1-(1-\theta_{a})^{2}).$

Because each term of $\theta$ is strictly positive, we can cancel the common
terms to get

$\displaystyle 1$ $\displaystyle=$ $\displaystyle
1-(1-\theta_{c})^{2}+(1-\theta_{c})\theta_{c}$ (11) $\displaystyle 1$
$\displaystyle=$
$\displaystyle\theta_{c}\theta_{b}+(1-(1-\theta_{b})^{2})+\theta_{a}(\theta_{b}+2\theta_{c})$
(12) $\displaystyle 1$ $\displaystyle=$
$\displaystyle\theta_{c}(\theta_{a}+2\theta_{b})+\theta_{b}(\theta_{a}+2\theta_{c})+1-(1-\theta_{a})^{2}).$

Equation (11) reduces to $2\theta_{c}^{2}-3\theta_{c}+1=0$ with solutions
$\theta_{c}=1$ or $1/2$. The solution with $\theta_{c}=1$ is excluded by
supposition, so we must have $\theta_{c}=1/2$. Now we have
$\theta_{a}=1/2-\theta_{b}$; inserting into (12), the equation reduces to
$2\theta_{b}^{2}-2\theta_{b}+\frac{1}{2}=0$ with solution
$\theta_{b}=\frac{1}{2}$. So $\theta_{a}=0$, contradicting the supposition.

Theorem 7 establishes a more general result, but we have given the simpler
proof above to demonstrate the style of “proof by contradiction” to be used
later.

##  4 The $k=2$ case.

###### Theorem 3

For every compact metric space $S$, the set of invariant distributions for
$\pi_{1,2}$ is the same as the set of invariant distributions for $\pi_{2,2}$.

Proof.  Given $\theta\in\mathcal{P}(S)$, the transition kernel
$K=K^{\theta,1,2}$ for $\pi_{1,2}$ can be written as a Radon-Nikodym density
w.r.t. $\theta$ as follows.

$\displaystyle\frac{dK(x,\cdot)}{d\theta(\cdot)}(y)$ $\displaystyle=$
$\displaystyle 2\theta\\{z:d(x,z)>d(x,y)\\}+\theta\\{z:d(x,z)=d(x,y)\\}$
$\displaystyle=$ $\displaystyle\int\left(2\cdot
1_{\\{z:d(x,z)>d(x,y)\\}}+1_{\\{z:d(x,z)=d(x,y)\\}}\right)\theta(dz).$

So the identity $\theta=\theta K$ characterizing a $\pi_{1,2}$-invariant
distribution $\theta$ can be written in density form as

$\displaystyle 1$ $\displaystyle=$ $\displaystyle\int\theta(dx)\
\frac{dK(x,\cdot)}{d\theta(\cdot)}(y)$ (13) $\displaystyle=$
$\displaystyle\int\int\left(2\cdot
1_{\\{z:d(x,z)>d(x,y)\\}}+1_{\\{z:d(x,z)=d(x,y)\\}}\right)\theta(dz)\theta(dx)$

where the equality holds for $\theta$-a.a. $y$. Because

$1=\int\int 1\ \
\theta(dz)\theta(dx)=\int\int\left(1_{\\{z:d(x,z)>d(x,y)\\}}+1_{\\{z:d(x,z)=d(x,y)\\}}+1_{\\{z:d(x,z)<d(x,y)\\}}\right)\
\ \theta(dz)\theta(dx)$

we have from (13) that

$\int\int 1_{\\{z:d(x,z)>d(x,y)\\}}\ \ \theta(dz)\theta(dx)=\int\int
1_{\\{z:d(x,z)<d(x,y)\\}}\ \ \ \theta(dz)\theta(dx).$ (14)

Analogous to (13), the identity characterizing a $\pi_{2,2}$-invariant
distribution $\phi$ can be written as

$1=\int\int\left(2\cdot
1_{\\{z:d(x,z)<d(x,y)\\}}+1_{\\{z:d(x,z)=d(x,y)\\}}\right)\phi(dz)\phi(dx).$
(15)

By (13) and (14), any $\pi_{1,2}$-invariant distribution $\theta$ satisfies
(15) and is therefore a $\pi_{2,2}$-invariant distribution. The converse holds
via the analog of (14) for $\phi$.

##  5 The case $S=[0,1]$

Numerical study in [1] suggests that there are no invariant distributions on
$[0,1]$ with full support, for any $(j,k)$. Theorem 4 proves a slightly
stronger result in the case $k=2$ (recall that by Theorem 3 the cases $j=1$
and $j=2$ here are identical). The stronger form is not true for general
$(j,k)$, for instance the uniform distribution on the 4 points
$\\{0,0.4,0.6,1\\}$ is invariant for $\pi_{3,4}$.

###### Theorem 4

There are no $\pi_{2,2}$-invariant distributions on $[0,1]$ other than those
of the form $\delta_{s}$ or $\delta_{s_{1},s_{2}}$.

By considering the endpoints of the support of an invariant distribution, and
scaling, this reduces to proving

> equivalent assertion: The only $\pi_{2,2}$-invariant distribution on $[0,1]$
> whose support contains both $0$ and $1$ is the distribution $\delta_{0,1}$.

We will prove this in two steps.

###### Lemma 5

There is no $\pi_{2,2}$-invariant distribution whose support contains $0$ and
which assigns zero weight to the point $0$.

Proof.  For a proof by contradiction, suppose such an invariant distribution
$\theta$ exists. Take $Y,Y_{1},Y_{2}$ independent with distribution $\theta$.
Invariance says that the random variable $X$ defined as

$\displaystyle X$ $\displaystyle=$ $\displaystyle Y_{2}\mbox{ if
}|Y-Y_{1}|<|Y-Y_{2}|$ $\displaystyle=$ $\displaystyle Y_{1}\mbox{ if
}|Y-Y_{2}|<|Y-Y_{1}|$

(with our usual convention about ties) will also have distribution $\theta$.
Fix $0<x<1$. From the definition we have the inclusion of events

$\\{X\leq x\\}\subseteq A_{1}\cup A_{2}\cup A_{3}$

where

$\displaystyle A_{1}$ $\displaystyle:=$ $\displaystyle\\{Y_{1}\leq x,\
Y_{2}\leq x\\}$ $\displaystyle A_{2}$ $\displaystyle:=$
$\displaystyle\\{Y_{1}\leq x,\ Y_{2}>x,\
Y\geq{\textstyle\frac{1}{2}}(Y_{1}+Y_{2})\\}$ $\displaystyle A_{3}$
$\displaystyle:=$ $\displaystyle\\{Y_{2}\leq x,\ Y_{1}>x,\
Y\geq{\textstyle\frac{1}{2}}(Y_{1}+Y_{2})\\}$

and the $(A_{i})$ are disjoint. Now note that

$A_{2}\subseteq\\{Y_{1}\leq x,\ Y\leq{\textstyle\frac{1}{2}}(x+Y_{2})\\}$

and similarly for $A_{3}$. So by independence, the distribution function $F$
of $\theta$ satisfies

$F(x)\leq F^{2}(x)+2F(x){\mathbb{P}}(Y\leq{\textstyle\frac{1}{2}}(x+Y_{2})).$
(16)

By hypothesis, $F(x)>0$ for $x>0$ and $F(x)\downarrow 0$ as $x\downarrow 0$.
So we can divide both sides of (16) by $F(x)$ and take limits as $x\downarrow
0$ and deduce

${\mathbb{P}}(Y\leq{\textstyle\frac{1}{2}}Y_{2})\geq{\textstyle\frac{1}{2}}.$

By symmetry we also have
${\mathbb{P}}(Y_{2}\leq{\textstyle\frac{1}{2}}Y)\geq{\textstyle\frac{1}{2}}$,
and so

${\mathbb{P}}({\textstyle\frac{1}{2}}Y<Y_{2}<2Y)=0.$

But this is impossible for i.i.d. samples from a distribution $\theta$ on
$(0,1]$, because it would remain true for $\theta$ conditioned on an interval
of the form $[y,3y/2]$.

Using Lemma 5 and reflection-symmetry of $[0,1]$, to prove the equivalent
assertion and hence Theorem 4 it will be sufficient to prove

###### Lemma 6

If $\theta$ is a $\pi_{2,2}$-invariant distribution and
$\theta_{0}>0,\theta_{1}>0$ then $\theta=\delta_{0,1}$.

Here we write $\theta_{s}$ for $\theta(\\{s\\})$.
Proof.  First note an elementary fact:

if $0<x<1$ and $\beta\geq 0$ and $x\geq x^{2}+2x(1-x)\beta$ then $\beta\leq
1/2$. (17)

From the construction with $(Y,Y_{1},Y_{2},X)$ we have

$\displaystyle\theta_{0}$ $\displaystyle=$
$\displaystyle\theta_{0}^{2}+2{\mathbb{P}}(Y_{1}=0,Y_{2}>0,Y>Y_{2}/2)+{\mathbb{P}}(Y_{1}=0,Y_{2}>0,Y=Y_{2}/2)$
(19) $\displaystyle=$
$\displaystyle\theta_{0}^{2}+2\theta_{0}\left({\mathbb{P}}(2Y>Y_{2},Y_{2}>0)+{\textstyle\frac{1}{2}}{\mathbb{P}}(2Y=Y_{2},Y_{2}>0)\right)$
$\displaystyle=$
$\displaystyle\theta_{0}^{2}+2\theta_{0}(1-\theta_{0})\left({\mathbb{P}}(2Y>Y_{2}|Y_{2}>0)+{\textstyle\frac{1}{2}}{\mathbb{P}}(2Y=Y_{2}|Y_{2}>0)\right)$
$\displaystyle\geq$
$\displaystyle\theta_{0}^{2}+2\theta_{0}(1-\theta_{0})\left({\mathbb{P}}(Y>1/2)+{\mathbb{P}}(Y=1/2,Y_{2}<1|Y_{2}>0)+{\textstyle\frac{1}{2}}{\mathbb{P}}(Y=1/2,Y_{2}=1|Y_{2}>0)\right)$
$\displaystyle=$
$\displaystyle\theta_{0}^{2}+2\theta_{0}(1-\theta_{0})\left({\mathbb{P}}(Y>1/2)+\theta_{1/2}\
[{\mathbb{P}}(Y_{2}<1|Y_{2}>0)+{\textstyle\frac{1}{2}}{\mathbb{P}}(Y_{2}=1|Y_{2}>0)]\right)$
$\displaystyle=$
$\displaystyle\theta_{0}^{2}+2\theta_{0}(1-\theta_{0})\left({\mathbb{P}}(Y>1/2)+\theta_{1/2}\
[{\textstyle\frac{1}{2}}+{\textstyle\frac{1}{2}}{\mathbb{P}}(Y_{2}<1|Y_{2}>0)]\right)$
$\displaystyle\geq$
$\displaystyle\theta_{0}^{2}+2\theta_{0}(1-\theta_{0})\left({\mathbb{P}}(Y>1/2)+{\textstyle\frac{1}{2}}{\mathbb{P}}(Y=1/2)\right).$

By hypothesis $\theta_{0}>0$, so by (17) we have
${\mathbb{P}}(Y>1/2)+{\textstyle\frac{1}{2}}{\mathbb{P}}(Y=1/2)\leq 1/2$.
However we have the analogous sequence of equalities and inequalities for
$\theta_{1}$, which imply
${\mathbb{P}}(Y<1/2)+{\textstyle\frac{1}{2}}{\mathbb{P}}(Y=1/2)\leq 1/2$, and
so we must have

${\mathbb{P}}(Y>1/2)+{\textstyle\frac{1}{2}}{\mathbb{P}}(Y=1/2)=1/2={\mathbb{P}}(Y<1/2)+{\textstyle\frac{1}{2}}{\mathbb{P}}(Y=1/2).$

The quantity (19) now equals $\theta_{0}$, so the inequalities at (LABEL:ttt1)
and (19) must in fact be equalities. In order for the inequality leading to
(19) to be an equality, we must have either $\theta_{1/2}=0$ or
${\mathbb{P}}(Y_{2}<1|Y_{2}>0)=0$. In the latter case, $\theta$ is supported
on $\\{0,1\\}$ and so $\theta=\delta_{0,1}$, as desired. So the remaining case
is $\theta_{1/2}=0$. In this case, for the inequality leading to (LABEL:ttt1)
to be an equality, we must have ${\mathbb{P}}(Y_{2}<2Y,Y<1/2|Y_{2}>0)=0$. But,
as at the end of the proof of Lemma 5, this can only happen if
${\mathbb{P}}(0<Y<1/2)=0$. By the analogous argument for $\theta_{1}$ we have
${\mathbb{P}}(1/2<Y<1)=0$, and so the distribution is supported on $\\{0,1\\}$
and must be $\delta_{0,1}$, as desired.

This line of argument can be extended to some other values of $(j,k)$ – see
[1].

##  6 A class of tree spaces

Figure 1: A BTL space $S$ with $|S|=7$.

In this section we consider binary222Essentially the same argument works
without the binary assumption. tree leaves (BTL), illustrated in Figure 1, as
a class of finite spaces. Here $S$ is the finite set of leaves; the edges have
lengths which serve to determine the distance between two leaves as the length
of the unique path between them; the edges also define $|S|-2$ branchpoints.
To “break symmetry” we assume

all distances $(d(s_{i},s_{j}),j\neq i)$ are distinct. (20)

We claim that, as suggested by the general picture from numerics, for $k=2$
there are no invariant measures other than the omnipresent ones. An invariant
measure supported on a subset of leaves is an invariant measure on the induced
spanning tree of that subset, so to prove that claim it suffices to prove

###### Theorem 7

On a BTL space $S$ with $|S|\geq 3$ and satisfying (20), and for $k=2$, there
are no invariant measures with full support.

Proof.  As in previous proofs, consider a $\pi_{1,2}$ \- invariant
distribution $\theta$ with full support on $S$, where $|S|\geq 3$. Take
$Y_{0},Y_{1},Y_{2}$ independent with distribution $\theta$. Invariance says
that the random variable $X$ defined as

$\displaystyle X$ $\displaystyle=$ $\displaystyle Y_{1}\mbox{ if
}d(Y_{0},Y_{1})<d(Y_{0},Y_{2})$ (21) $\displaystyle=$ $\displaystyle
Y_{2}\mbox{ if }d(Y_{0},Y_{2})<d(Y_{0},Y_{1})$

will also have distribution $\theta$. We proceed to a proof by contradiction.

We quote an elementary fact.

###### Lemma 8

For any probability distribution $\theta$ on a BTL space $S$, either
(i) $\theta(s_{0})>\frac{1}{2}$ for some $s_{0}\in S$
or (ii) there exists a centroid, that is a branchpoint such that the
associated partition $S=\cup_{i=1}^{3}A_{i}$ of leaves satisfies
$0<\theta(A_{i})\leq\frac{1}{2}$ for all $i$.

Consider case (i). That is, suppose $\theta$ is invariant and
$\theta(s_{0})\in(\frac{1}{2},1)$. From the invariance relation (21), in order
that $X=s_{0}$ it is sufficient that

$(Y_{0}\neq s_{0},Y_{1}=s_{0},Y_{2}=s_{0})\mbox{ or }(Y_{0}=s_{0},Y_{1}\mbox{
or }Y_{2}=s_{0}).$

So, setting $\theta(s_{0})=x\in(\frac{1}{2},1)$,

$x\geq(1-x)x^{2}+x(2x-x^{2}).$

Cancelling $x$, this reduces to $2x^{2}-3x+1\geq 0$, but this inequality is
false for $x\in(\frac{1}{2},1)$.

Now consider case (ii). There is a centroid branchpoint defining a partition
$S=\cup_{i=1}^{3}A_{i}$. Consider the leaf $s_{1}$ which is closest to the
centroid. We may assume $s_{1}\in A_{1}$. From the invariance relation (21),
in order that $X=s_{1}$ it is sufficient that the following condition (*)
holds:

> exactly one of $(Y_{1},Y_{2})$ equals $s_{1}$
> and
> $Y_{0}$ and the other333The leaf from $(Y_{1},Y_{2})$ that is not $s_{1}$.
> $Y$ are in different components of $\cup_{i=1}^{3}A_{i}$.

For instance, if $Y_{0}\in A_{2}$ and $Y_{1}=s_{1}$ and $Y_{2}\in A_{3}$, then
$Y_{2}$ is some leaf in $A_{3}$ which is farther from the centroid than is
$s_{1}$, so $d(Y_{0},s_{1})<d(Y_{0},Y_{2})$. The other possibilities are
similar.

By considering the three possibilities for “ different components of
$\cup_{i=1}^{3}A_{i}$” we see that the probability of (*) equals
$\theta(s_{1})$ times

$2\theta(A_{1})\theta(A_{2})+2(\theta(A_{1})-\theta(s_{1}))\theta(A_{2})$
$+2\theta(A_{1})\theta(A_{3})+2(\theta(A_{1})-\theta(s_{1}))\theta(A_{3})$
$+4\theta(A_{2})\theta(A_{3})$

which rearranges to

$\theta(s_{1})[4(\theta(A_{1})\theta(A_{2})+\theta(A_{1})\theta(A_{3})+\theta(A_{2})\theta(A_{3}))-2\theta(s_{1})(\theta(A_{2})+\theta(A_{3}))]$
$=\theta(s_{1})\cdot B,\mbox{ say.}$ (22)

A disjoint sufficient condition for $X=s_{1}$ is that $Y_{1}=Y_{2}=s_{1}$,
which has probability $\theta^{2}(s_{1})$. So

$\theta(s_{1})={\mathbb{P}}(X=s_{1})\geq\theta(s_{1})(B+\theta(s_{1})).$

Cancelling the $\theta(s_{1})$ term,

$1\geq
4(\theta(A_{1})\theta(A_{2})+\theta(A_{1})\theta(A_{3})+\theta(A_{2})\theta(A_{3}))-2\theta(s_{1})(\theta(A_{2})+\theta(A_{3})-{\textstyle\frac{1}{2}}).$

Because $\sum_{i}\theta(A_{i})=1$ we have
$\theta(A_{2})+\theta(A_{3})-{\textstyle\frac{1}{2}}={\textstyle\frac{1}{2}}-\theta(A_{1})$
and

$2(\theta(A_{1})\theta(A_{2})+\theta(A_{1})\theta(A_{3})+\theta(A_{2})\theta(A_{3}))=1-\sum_{i}\theta^{2}(A_{i})$

and the inequality above reduces to

$1\geq
2-2\sum_{i}\theta^{2}(A_{i})-2\theta(s_{1})({\textstyle\frac{1}{2}}-\theta(A_{1})).$

Because $\theta(s_{1})\leq\theta(A_{1})$, this implies

$C:=\sum_{i}\theta^{2}(A_{i})+\theta(A_{1})({\textstyle\frac{1}{2}}-\theta(A_{1}))\geq{\textstyle\frac{1}{2}}.$
(23)

We need to show that $C\geq{\textstyle\frac{1}{2}}$ cannot in fact occur under
the constraints $0<P(A_{i})\leq{\textstyle\frac{1}{2}}$ and
$\sum_{i}\theta(A_{i})=1$. Given $\theta(A_{1})=x$, the quantity $C$ is
maximized when
$(\theta(A_{2}),\theta(A_{3}))=({\textstyle\frac{1}{2}},{\textstyle\frac{1}{2}}-x)$
and so

$C\leq
x^{2}+{\textstyle\frac{1}{4}}+({\textstyle\frac{1}{2}}-x)^{2}+x({\textstyle\frac{1}{2}}-x)=x^{2}-{\textstyle\frac{1}{2}}x+{\textstyle\frac{1}{2}}.$

This implies that $C<{\textstyle\frac{1}{2}}$ on the open interval
$x\in(0,{\textstyle\frac{1}{2}})$, and we cannot have $x=0$ or
${\textstyle\frac{1}{2}}$ by the $\theta(A_{i})>0$ constraint.

## References

  * [1] David J. Aldous, Madelyn Cruz, and Shi Feng.  Markov chains and mappings of distributions on compact spaces II: Numerics and conjectures, 2024.  arXiv:2403.18153.
  * [2] David A. Levin and Yuval Peres.  Markov chains and mixing times.  American Mathematical Society, Providence, RI, 2017.  Second edition of [ MR2466937], With contributions by Elizabeth L. Wilmer, With a chapter on “Coupling from the past” by James G. Propp and David B. Wilson.


# Task-Independent Spiking Central Pattern Generator: A Learning-Based
Approach

Elie Aljalbout  Technische Universität München  Florian Walter  Technische
Universität München  Florian Röhrbein  Technische Universität München  Alfred
Kärcher SE Co. & KG  Alois Knoll  Technische Universität München

###### Abstract

Legged locomotion is a challenging task in the field of robotics but a rather
simple one in nature. This motivates the use of biological methodologies as
solutions to this problem. Central pattern generators are neural networks that
are thought to be responsible for locomotion in humans and some animal
species. As for robotics, many attempts were made to reproduce such systems
and use them for a similar goal. One interesting design model is based on
spiking neural networks. This model is the main focus of this work, as its
contribution is not limited to engineering but also applicable to
neuroscience. This paper introduces a new general framework for building
central pattern generators that are task-independent, biologically plausible,
and rely on learning methods. The abilities and properties of the presented
approach are not only evaluated in simulation but also in a robotic
experiment. The results are very promising as the used robot was able to
perform stable walking at different speeds and to change speed within the same
gait cycle.

##  1 Introduction

Robotic locomotion is a functionality that enables robots to navigate and
transport themselves. One specific form of locomotion is legged locomotion.
Despite the challenges faced to achieve this form, it presents multiple
advantages such as adaptability to multiple types of terrains (e.g. unpaved
ground, rugged terrain, stairs etc.), interaction with the physical
environment, agility in terms of generated paths complexity, and superior
failure tolerance (e.g. when a leg fails the robot can still move) [14].

In biology, humans and animals can move and transport their bodies from one
place to another in a quite efficient manner, which motivates and highlights
the necessity and importance of approaching this problem with more inspiration
from biology. In addition, imitating biological locomotion could also lead to
more human- and animal-like gaits. However, the functionality leading to
biological locomotion is not yet fully understood. Central Pattern Generators
(CPGs) are neural networks that are thought to be responsible for multiple
rhythmic behaviors, including locomotion. The exact architecture and
configuration leading to the ability of these networks to produce rhythmic
patterns are still not entirely known. Consequently, this field of research is
considered multidisciplinary, since it contributes to both robotics and
neuroscience. Another advantage of engineering CPGs is the possibility to
implement them on neuromorphic hardware [23], leading to more power
efficiency.

In robotics, many attempts were made to design CPGs for locomotion control.
These attempts can be separated into two main groups:

  1. 1.

Methods based on systems of coupled oscillators (SCO)

  2. 2.

Spiking neural networks (SNN) based methods

The second type is the one adopted in this work. As for SCOs, they offer the
ability to be mathematically analyzed due to their mathematical nature, which
made them advance faster than their counterparts that present complex dynamics
and harder analytical obstacles. Unfortunately, SCOs lack biological
plausibility, which means that their contributions to neuroscience and
interventions in neuro-prosthetics are limited. For example, when interfacing
prosthetic robotic devices to amputated humans and spinal cord injury
patients, SNN-based approaches process and receive the same kind of signals as
the CPGs of their biological counterparts [7]. The same doesn’t apply for
SCOs.

In this paper, we introduce a novel SNN-based CGG architecture that combines
biological plausibility with easy tunable output behavior. The architecture
builds on previously published models but adds three notable features: First,
the gait speed can be easily modulated through a tonic input signal even
within an ongoing cycle. Second, the generated patterns are not encoded
manually or through evolutionary algorithms but directly learned based on a
target signal. Third, the number of phases can be defined freely, which,
together with the second feature, allows the approach to be general and task-
independent. As we will show during the evaluation, this approach is very
stable and converges quickly, making it very well-suited for implementing it
in robots. Importantly, while also other CPG architectures proposed earlier
support speed modulation, pattern learning, or flexible number of phases, none
of them supports all of these properties at the same time. Our contribution
therefore can both provide valuable insight to both neuroscience-focused CPG
research and robot locomotion research.

The rest of this paper is structured as follows. In section 2, related methods
will be presented, followed by the proposed approach in section 3. Results and
conclusion will be then discussed in sections 4 and 5 respectively.

##  2 Related Work

Several approaches have been demonstrated to be capable of building CPGs using
spiking neural networks. Some of these methods are simply based on hand-tuning
of the synaptic efficacies [3, 12, 13, 15, 5, 8]. Such methods manage to
reproduce basic CPG functionalities such as stable walking and in some cases
speed control [13, 15]. The main disadvantage of these approaches is the small
size of their networks111CPGs are known to have a small number of neurons
relative to the whole nervous system. The discussed methods use networks that
are far smaller than biological CPGs. which makes the weight tuning process
easier but forbids the network from scaling to additional features such as
modulation by sensory feedback and autonomous gait transition. Another set of
approaches relies on genetic algorithms and similar methods to learn the
network weights. For instance, the work in [20] applied genetic evolution on
the network from [13] and succeeded to obtain a similar behavior, omitting the
complicated tuning of the locomotion behavior but keeping the same
disadvantages as previously discussed. An alternative approach was recently
presented in [7, 6], where the CPG network is formulated as a Christiansen
grammar and an evolutionary process is performed on it. Interestingly, this
method not only learns the network’s weights but also its architecture, which
is a great advantage offering more automation for the CPG building process.
Further approaches based on evolutionary algorithms were also proposed in [17]
and [4]. Besides the obvious disadvantage of slow convergence, evolutionary
methods don’t scale to larger networks and their training method is not
biologically plausible, which also limits their neuroscientific contribution.
Other methods also fall in the category of non-plausible to biology such as
[22], where it is proposed to perform brute-force search on the CPG parameters
space, and [11] where a gradient descent based method is presented. These
methods have multiple disadvantages such as slow convergence (especially in
the case of brute-force) and being incapable to scale to advanced CPG
features. Finally, the method in [19] presents a general framework to build
task-independent CPGs. It relies on the remote supervision method (ReSuMe)
[18] to learn the weights connecting two populations of randomly connected
neurons to the motor neurons of the CPG. This method uses biologically
plausible learning but a non-plausible network architecture. Additionally, it
only achieves stable walking at a fixed speed and omits additional CPG
features.

Figure 1: The proposed network architecture for a CPG with two phases. The
neural phase generator (NPG) generates a rhythm corresponding to the number of
desired phases. The H neuron is the representative of the activity of each
phase, Q takes care of inhibiting other phase modules and T is responsible for
the transition to another phase. This rhythm is then propagated to the pattern
forming networks (PFN) which map the NPG outputs to a higher dimensional
space. PFNs output is then fed to motor neurons which activate the robot’s
muscles.

##  3 Proposed Approach

As previous methods either relied on hand-tuning, neglected advanced CPG
features, or couldn’t scale to larger CPG networks, we present the following
approach, which is the first to successfully combine all of these
characteristics using a biologically plausible network architecture. From a
general perspective, the architecture consists of three main components. The
first one is the neural phase generator (NPG), which is adapted from [15] with
minor modifications. In analogy to the biological model in [21], this part
plays the role of the rhythm generator. The second part is the pattern forming
network (PFN), which consists of networks with random parameters. At this
level spiking patterns appear in a rhythmic fashion with respect to the NPG’s
output. Finally, the last part is the pool of motor neurons, which activates
the actuated muscles through spike activity. Each of these components will be
discussed next in more details.

Figure 1 shows the proposed architecture in the case of a CPG which encodes
only two phases. For this reason, the architecture is almost symmetric, when
omitting the randomness in the PFNs. In the case when more phases are
required, the architecture is extended in all three parts. For instance, at
the NPG level, more modules should be added corresponding to the new
requirement (Figure 2). Similarly, at the second level, random networks are
also created with respect to the number of desired phases. Finally, in the
pool of motor neurons, the number of neurons can also be different but not
necessarily. In the following, each of these components will be described
separately, in terms of their inner connections, functionality, and role.

###  3.1 Neural Phase Generator

Figure 2: Illustration of an NPG network for the case of a CPG having four
phases.

As already mentioned, this part is adapted from [15]. Briefly, each NPG module
consists of three neurons H, Q and T (as seen in Figure 2). The H neuron is
the representative of the module’s activity. In other words, whenever a
module’s H neuron is spiking, the phase this module represents is active. As
for the Q neuron, its role is to guarantee that other NPG modules are inactive
when required. It does so by inhibiting the H neurons of the corresponding
modules. As for T neurons, they are responsible for the transition from module
to module. In contrast to the NPG in [15], our modified version receives an
external tonic input, which is responsible for the start and end of activity
by exciting and inhibiting the neurons of the NPG. This input represents the
modulating inputs issued from the brain to control the CPG’s activity. Once
the NPG is active, cutting the external input won’t stop its activity.
Instead, the network will keep oscillating because of its inner mutual
connections. This behavior is compatible with the nature of CPGs being able to
oscillate even in the absence of input from the brain. To make sure this
property is available, a constraint is required: the auto-synapse of the H
neurons should have weights that are high enough to ensure a continuous
spiking behavior until the respective T neurons fire, which then results in a
phase transition. Moreover, when the tonic input has a higher frequency, the H
neuron will fire more frequently resulting in a faster transition to the next
phase. This modification gives control to the NPG over the speed and phase
properties of the produced gaits. This behavior is plausible with the proposed
biological model where this control is available at the level of the rhythm
generator. In [15], this ability is given to the second layer of the
architecture which in their case is the motor output shaping stage. Compared
to our architecture this would correspond to giving control to the PFNs.

Another modification, to the NPG of [15], is in the number of T neurons in
each module. In the previous work, two T neurons were present in every module,
as for ours, this number is variable. It is dependent on the duration of the
phase each module is responsible for, the longer the duration the bigger the
number. All T neurons except the last one, receive the same external tonic
input (as the H neuron), which is important since they each excite a part of
the PFN. However, these intermediate T neurons, in addition to their tonic
inputs are not visible in figure 2 for the sake of simplicity.

Finally, the number of modules in our version of the NPG is not fixed to two,
but instead is flexible depending on the number of desired phases, as
illustrated in figure 2. This flexibility allows the architecture to be
general and adaptable for different types of locomotion.

###  3.2 Pattern Forming Networks (PFN)

This part of the architecture consists of several separate networks. The
number is equivalent to the number of phases of the desired CPG. Each of these
networks corresponds to a module of the NPG and is excited by its H and
intermediate T neurons. The neurons within each of these networks have
randomly initialized properties, leading to richer dynamics and higher
learning abilities.

In order to simplify the learning procedure, a minor restriction is enforced:
neurons of the PFNs are required to spike only once during each CPG cycle.
This is achieved by creating another network with a similar number of neurons
for every PFN. This second network is named inhibiting network (IN).

Figure 3: Illustration of inhibiting networks. Neurons in the left column are
the inhibiting neurons, and the ones on the right are those of the PFN.

Both networks are then connected in a one-to-one fashion with excitation
issued from the PFN to the IN and with inhibition in the opposite direction.
INs have no synapses connecting their neurons, but instead, exciting auto-
synapses for all of their neurons to ensure that each one of them keeps
spiking once it’s activated. To illustrate how this works: each time a neuron
in a PFN spikes, it’ll excite a corresponding neuron in an IN leading it to
spike in a continuous fashion and inhibit reciprocally the same PFN neuron
that made it fire. This concept is illustrated in figure 3. Now to make sure
that the IN is reset in the next cycle, it’ll be inhibited by a Q neuron of an
NPG module from a different phase.

Additionally, neurons within a certain PFN are connected with random mutual
inhibition at a $10\%$ rate of connectivity (i.e. the number of synapses is
$10\%$ the one of an all-to-all connectivity scheme), which ensures that
spikes produced by different neurons within the network are well distributed
over the whole duration of the corresponding phase, and not concentrated in a
small duration right after the NPG spikes. Distribution of the PFN spikes in
time is very important for the learning, in a way that the more spikes are
distributed, the more the number of desired behaviors that can be learned
increases. The intuition behind that is simple. To illustrate it, consider the
case where no spike occurs in the time interval $[t,t+\epsilon]$ at the PFN.
This will make it less probable that any spike occurs in this same interval in
the pool of motor neurons, which are only excited by PFN neurons. In
conclusion, PFN spikes should be the most possibly scattered over the phase
duration. To enforce that, the size of the network can be increased in order
to obtain a higher spiking rate. Conceptually, the role of the PFN can be seen
as mapping the NPG spikes into a higher dimensional space in which spikes are
well distributed in time. Finally, the external control inputs of the NPG are
also propagated to this layer. When the frequency of these signals is
increased, the frequency of spiking in the H and intermediate T neurons in the
NPG will also increase. Thus, the neurons of the PFN will fire sooner, which
means that also the speed of gait will increase. The relation between the
frequency of spiking at the NPG level and the time of spikes in the PFN is not
a perfect mapping, but that is left for the learning to figure out. Concerning
the learning, it is important to note that it’s dependent on the PFN’s size.
For instance, the larger the PFN is, the more probable it is that the learning
method will converge. However, the convergence will then require longer
training time. Throughout our experiments, we used PFNs with sizes ranging
from $150$ to $500$ neurons. The exact selection of the number is also
dependent on each phase’s duration.

###  3.3 Pool of Motor Neurons

Figure 4: Pool of motor neurons, which receive excitatory and inhibitory
inputs from the PFN. The weights of the synapses connecting PFNs and these
neurons are learned using the remote supervision method ReSuMe.

This last part of the architecture is where rhythmic patterns really matter,
since the spikes at this level are responsible for muscle activation and
producing gaits for legged locomotion. In biology, the precise way these
spikes are decoded into muscle activation is not fully known, but research has
indicated that population rate coding is used [9, 16]. In other terms, the
activation of a muscle is encoded by a population of neurons. Their spike rate
is related to the level of activation of the muscle. Therefore, the number of
neurons in the pool is dependent on the desired locomotion and the number of
actuated muscles. A similar scheme was adopted in the robotic experiment
performed within this work. Moreover, neurons within the pool are not inter-
connected, but receive as input the spikes produced by the PFNs, through
synapses of both excitatory and inhibitory types. The weights of those
synapses are the only ones that are learned. They are adapted to obtain the
desired locomotion behavior. The learning procedure will be discussed later.
Finally, the control signals (external tonic input) applied at the NPG are
also propagated to the neurons in this pool. For instance, when the tonic
input has an increased frequency, the speed of the spike pattern produced by
the PFN will increase. Ideally, this will lead to the spiking behavior of the
motor neurons to be accelerated.

###  3.4 Learning

When designed faithfully to the constraints already mentioned, the previously
described architecture is capable by itself of producing rhythmic spiking
patterns in the pool of motor neurons. However, this rhythmic spiking behavior
is random and does not correspond to any locomotive behavior. Therefore, the
network should be trained somehow to adapt to the desired spiking
behavior222The desired spiking behavior is specified by the CPG designer and
can also be given by demonstration. in order for it to be used in real robot
locomotion tasks. In order to do so, the Remote Supervision Method (ReSuMe)
[18], was employed. However, the only synapses of the network which learning
is applied on, are the ones connecting the PFNs to the pool of motor neurons.

The way ReSuMe works is as follows [18]: For each of the learning neurons (in
our case the motor neurons), a teacher signal is associated with predetermined
timing, representing the desired spiking behavior. This signal is not
delivered to the learning neuron, but still plays an important role in the
learning updates of synapses terminating at it. Moreover, the modification is
based on two rules, the first one depends on the correlation between the
presynaptic and the desired spike times, as for the second, it depends on the
correlation between the presynaptic and the postsynaptic spike times [18]. The
following is the corresponding modification function for the synaptic weight
between a presynaptic neuron $k$ (from the PFN) and a postsynaptic neuron $i$
(from the motor neurons):

$\begin{split}\frac{d}{dt}w_{ki}(t)=&S^{d}(t)\Big{[}a^{d}+\int_{0}^{\infty}W^{d}(s^{d})S^{in}(t-s^{d})ds^{d}\Big{]}\\\
+&S^{l}(t)\Big{[}a^{l}+\int_{0}^{\infty}W^{l}(s^{l})S^{in}(t-s^{l})ds^{l}\Big{]}\end{split}$
(1)

Where $S^{d}$, $S^{l}$ and $S^{in}$ are respectively the target, the post and
the presynaptic spike trains, $a^{d}$ and $a^{l}$ determine the so-called non-
Hebbian processes of weight modifications, and $s^{l}$ and $s^{d}$ represent
respectively the difference between the time of spike of the postsynaptic
neuron and the presynaptic one and the difference between the time of the
teacher signal and the time of spike of the presynaptic neuron. In the case of
excitatory synapses, the terms $a^{d}$ and $a^{l}$ are positive, and negative
otherwise. As for $W^{d}$ and $W^{l}$, they represent the learning windows,
and are formulated as follows:

$W^{d}(s^{d})=\begin{cases}\ +A^{d}\cdot
exp(\frac{-s^{d}}{\tau{}^{d}}),&\text{if }s^{d}>\,0\\\ 0,&\text{if }s^{d}\leq
0\end{cases}$ (2) $W^{l}(s^{l})=\begin{cases}\ -A^{l}\cdot
exp(\frac{-s^{l}}{\tau{}^{l}}),&\text{if }s^{l}>\,0\\\ 0,&\text{if }s^{l}\leq
0\end{cases}$ (3)

Where $W^{d}$ and $W^{l}$ correspond respectively to the learning window of
the target and postsynaptic neurons. As for $A^{d}$, $A^{l}$, $\tau^{d}$ and
$\tau^{l}$, they are all constants, such that $A^{d}$, $A^{l}$ are positive
for excitatory synapses and negative otherwise and $\tau^{d}$, $\tau^{l}$ are
always positive. When setting $a^{d}=-a^{l}=a$, $\tau^{d}=\tau^{l}$ and
$A^{d}=A^{l}$, then equation (1) takes the following form:

$\frac{d}{dt}w_{ki}(t)=[S^{d}(t)-S^{l}(t)]\Big{[}a^{d}+\int_{0}^{\infty}W^{d}(s)S^{in}(t-s^{d})ds^{d}\Big{]}$
(4)

Intuitively, what the method does to the network is to strengthen (weaken)
synapses which are excitatory (inhibitory) and are incoming to a certain motor
neuron, when a spike is transmitted through these synapses within the learning
window as defined by $\tau$, whenever the subject motor neuron is desired to
spike. This rule pushes the motor neurons to spike at the desired times.
Additionally, it weakens (strengthens) synapses which are excitatory
(inhibitory) and are incoming to a certain motor neuron when a spike is
transmitted through these synapses within the learning window defined by
$\tau$, whenever the motor neuron actually spikes. As for this rule, it
forbids motor neurons from spiking at undesired times.

##  4 Experiments

We designed experiments to answer the following questions:

  * •

To which extent is the presented approach capable of learning desired spiking
patterns?

  * •

Is the targeted speed modulation property working after the learning procedure
is done?

  * •

Can the system control the locomotion of legged robots?

For the first two questions, we designed experiments in a simulation-only
environment as described in section 4.2. Additionally, we built a real-world
experiment with a legged robot (section 4.3) to answer the last question and
to validate the second one again. The setup for both experiments is described
in section 4.1.

Properties |  | External
---
input
current
(pA)
| Resting
---
membrane
potential
(mV)
| Capacity
---
of the
membrane
(pF)
| Membrane
---
time
constant
(ms)
| Spike
---
threshold
(mV)
| Reset
---
potential
(mV)
| min | max | min | max | min | max | min | max | min | max | min | max
CPG | 0.0 | -70.0 | 250.0 | 10.0 | -55.0 | -70.0
PFN | 0.0 | 150.0 | -90.0 | -70.0 | 100.0 | 300.0 | 9.0 | 30.0 | -50.0 | -30.0 | -90.0 | -60.0
Motor | 0.0 | -70.0 | 250.0 | 10.0 | -55.0 | -70.0

Table 1: Parameters of the spiking neurons used in the experiments.

###  4.1 Setup

In order to allow our results to be reproducible, we provide the
implementation details we used along with our experiments. Namely, we used the
NEST simulator [10]. In all of our experiments, the leaky integrate-and-fire
neuron model [1] was used with alpha-function shaped synaptic currents. As
shown in Table 1, all neurons used the NEST default parameters except the ones
of the PFN, which are initialized with random uniform values. Parameters that
are not mentioned in Table 1 are set to default. As for synapses, they had
different initialization ranges at different levels. At the PFN level, neurons
were connected to each other with synapses initialized with uniformly random
weights between $-3$ and $-1$. Neurons of the PFN and motor neurons are
connected together with an all-to-all connectivity pattern. We experimented
with several numbers of PFN neurons between $300$ and $3000$ and all values
worked out. It is important to note that a larger number of neurons allows a
smaller spike shift error, however leads to a longer training time. $20\%$ of
the total number of neurons at this level are inhibitory and have initial
weights initialized with uniform random values between $-25$ and $-1$. The
rest of the neurons are excitatory and have their weights initialized with
uniform random values between $1$ and $5$.

As for the learning, we tried several combinations of hyper-parameters for
ReSuMe and report the following ones as being the best in terms of speed and
accuracy of learning. For learning windows $W^{d}$ and $W^{l}$ we used $3ms$.
Additionally we found out that $a=3$, $A^{d}=A^{l}=6$ and
$\tau_{d}=\tau_{l}=2$ perform best.

###  4.2 Simulated experiments

|
---|---
$s=100$ | $s=150$
|
$s=200$ | $s=250$
Figure 5: Spiking patterns obtained for CPGs with $s$ motor neurons for a
single CPG cycle. Red dots represent desired spiking patterns randomly
generated during the testing phase. Blue dots represent the patterns learned
using the presented method.

In order to initially validate the proposed method, we tested it in
simulation. Namely, we applied the method to more than $250$ randomly
generated target spiking sequences with different numbers of motor neurons and
let the network learn to imitate these sequences. We experimented with pools
of motor neurons of sizes ranging from $5$ to $250$. As a result, the method
managed to successfully learn how to imitate all presented spiking sequences
without any need for parameter hand-tuning. The average spike shift error of
all experiments was very low ($\tt{\sim}2.1ms$). Figure 5 shows four different
spiking patterns that have been successfully learned at the level of the motor
neurons. These patterns are generated by our CPG within a single cycle. These
results show that the proposed architecture together with ReSuMe enables
learning any desired locomotive behavior with minimal design effort.
Consequently, it is possible to embed any desired locomotive gait (whether
measured from an already existing walking system or self-defined) in a spiking
neural network.

|
---|---
$s=25$, $F=250$ |  $s=25$, $F=500$
|
$s=150$, $F=250$ |  $s=150$, $F=500$
Figure 6: Speed modulation of two CPGs with different number of motor neurons
$s$ in each row. When increasing the spiking frequency $F$ of the tonic input
from 250 spikes/s to 500 spikes/s we obtain an increased speed in the motor
output. Vertical lines indicate the end of a spiking cycle.

Additionally, the obtained CPGs are capable of acceleration and deceleration
of the same spiking sequences when presented with a higher and lower frequency
tonic inputs respectively. Figure 6 illustrates this property. In the first
row, a CPG with a small number of motor neurons is shown, as it’s easier to
watch the modulation property with a sparse plot. However, our method is also
capable of modulating larger numbers of motor neurons as shown in the second
row. In most cases, such results were obtained by training the CPG network
using only one tonic input frequency. However, in some rare cases, speed
modulation was only obtained after training the network under different tonic
input frequencies.

(a) $F=250$ $spikes/s$

(b) $F=500$ $spikes/s$

Figure 7: Outputs of the motor neurons used in the robot experiment under
different tonic input frequencies $F$.

###  4.3 Robotics experiment

We tested our approach on a robot (Allbot VR408 [2] shown in figure 8). The
Allbot is a simple quadruped robot with four legs and eight joints that are
actuated by standard servo motors. It is equipped with an Arduino Mega board
that is used in the experiment to relay motion commands from the network
simulation running on the host computer to the servos. The outputs of the
motor neurons were used to encode the joint angles. Namely, we used population
rate coding, such that the conversion from spikes to angles is done based on:

$\theta_{i}=\frac{n_{i}}{N_{i}}*\pi$ (5)

where $\theta_{i}$ is the angle of the joint $i$, $n_{i}$ is the number of
spikes of the motor neurons of the joint $i$, and $N_{i}$ is the total number
of neurons representing the $i^{th}$ joint. Each joint is represented by a set
of 25 neurons. Activations within a small window of time333The time window of
spike decoding is dependent on the control frequency of the robot correspond
to the angle command. To generate the desired spiking pattern, we used
equation (5) to obtain the number of spikes for each joint $n_{i}$, and
generated $n_{i}$ spikes at different randomly chosen neurons encoding joint
$i$. The motor activations obtained in this experiment are illustrated in
Figure 7. Similarly, other control types (such as joint velocity control) can
be applied. After learning a predefined motor spiking behavior, the robot is
capable of stable walking at different speeds corresponding to different tonic
inputs. We provide a video of this
experiment444https://sites.google.com/view/task-independent-cpg/home. The code
for all experiments will be made available upon publication of this paper.

(a)

(b)

Figure 8: (a) ALLBOT VR408. (b) Joints of a single robot leg.

##  5 Conclusion

Spiking CPGs are promising techniques for robot locomotion. Previous methods
have either missed the generality feature, scalability, biological
plausibility or required parameter hand-tuning. In this work, we were able to
present the first general framework for building spiking CPGs that includes
all of the desired features. Namely, the framework allows to build CPGs for
different types of robots with a minimal design effort, and that is done using
a learning algorithm. The used learning method is the remote supervision
method which was chosen for its fit to this work’s needs and its proved
convergence and biological plausibility. Additionally, the resulting CPGs have
architectures, behaviors, and features similar to biological counterparts. For
instance, they can be modulated by external signals to change the speed of
walking even during the process of locomotion. We also showed how our CPG can
successfully learn a target spiking behavior allowing it to achieve desired
gaits on target robots. A possible extension of this work would be to use
reinforcement learning to learn the weights of the network, eliminating the
need for behavior demonstration.

## Acknowledgements

This project/research has received funding from the European Union’s Horizon
2020 Framework Programme for Research and Innovation under the Specific Grant
Agreements No. 720270 (Human Brain Project SGA1) and No. 785907 (Human Brain
Project SGA2).

## References

  * [1] Abbott, L. F. (1999). Lapicque’s introduction of the integrate-and-fire model neuron (1907). Brain research bulletin, 50(5-6), 303-304.
  * [2] A robotic toolkit for building quadruped and biped robots, consists of 3D printing and servo motors. http://www.allbot.eu
  * [3] Billard, A., Ijspeert, A.J.: Biologically inspired neural controllers for motor control in a quadruped robot. In: Proceedings of the IEEE-INNS-ENNS International Joint Conference on Neural Networks, vol. 6, pp. 637–641. IEEE (2000)
  * [4] Caamano, P., Becerra, J.A., Bellas, F., Duro, R.J.: Using spiking neural networks for the generation of coordinated action sequences in robots. In: International Conference on Neural Information Processing, pp. 1013–1020. Springer (2008)
  * [5] Donati, E., Corradi, F., Stefanini, C., Indiveri, G.: A spiking implementation of the lamprey’s central pattern generator in neuromorphic vlsi. In: Biomedical Circuits and Systems Conference, pp. 512–515. IEEE (2014)
  * [6] Espinal, A., Rostro-Gonzalez, H., Carpio, M., Guerra-Hernandez, E.I., Ornelas-Rodriguez,M., Puga-Soberanes, H.J., Sotelo-Figueroa, M., Melin, P.: Quadrupedal robot locomotion:A biologically inspired approach and its hardware implementation. Computational Intelligence and Neuroscience (2016)
  * [7] Espinal, A., Rostro-Gonzalez, H., Carpio, M., Guerra-Hernandez, E.I., Ornelas-Rodriguez,M., Sotelo-Figueroa, M.: Design of spiking central pattern generators for multiple locomotion gaits in hexapod robots by christiansen grammar evolution. Frontiers in Neurorobotics 10 (2016)
  * [8] Garcia-Saura, C.: Central pattern generators for the control of robotic systems. ArXiv preprint arXiv:1509.02417 (2015)
  * [9] Gerstner, W., Kreiter, A.K., Markram, H., Herz, A.V.: Neural codes: firing rates and beyond. Proceedings of the National Academy of Sciences 94(24), 12740–12741 (1997)
  * [10] Gewaltig, M.O., Diesmann, M.: Nest (neural simulation tool). Scholarpedia 2(4), 1430 (2007)
  * [11] Kuroe, Y., Iima, H.: A learning method for synthesizing spiking neural oscillators. In:International Joint Conference on Neural Networks, pp. 3882–3886. IEEE (2006)
  * [12] Lewis, M.A., Hartmann, M.J., Etienne-Cummings, R., Cohen, A.H.: Control of a robot leg with an adaptive avlsi CPG chip. Neurocomputing 38, 1409–1421 (2001)
  * [13] Lewis, M.A., Tenore, F., Etienne-Cummings, R.: CPG design using inhibitory networks. In: Proceedings of the 2005 IEEE International Conference on Robotics and Automation, pp. 3682–3687. IEEE (2005)
  * [14] Machado, T., Silva, M.F.: An overview of legged robots. In: International symposium on mathematical methods in engineering. MME Press Ankara, Turkey (2006)
  * [15] Maufroy, C., Kimura, H., Takase, K.: Towards a general neural controller for quadrupedal locomotion. Neural Networks 21(4), 667–681 (2008)
  * [16] Maunsell, J.H., Van Essen, D.C.: Functional properties of neurons in middle temporal visual area of the macaque monkey. i. selectivity for stimulus direction, speed, and orientation. Journal of neurophysiology 49(5), 1127–1147 (1983)
  * [17] Orchard, G., Russell, A., Mazurek, K., Tenore, F., Etienne-Cummings, R.: Configuring silicon neural networks using genetic algorithms. In: International Symposium on Circuits and Systems, pp. 1048–1051. IEEE (2008)
  * [18] Ponulak, F.: Resume - new supervised learning method for spiking neural networks. Institute of Control and Information Engineering, Poznan University of Technology (2005)
  * [19] Ponulak, F., Belter, D., Kasinski, A.: Adaptive central pattern generator based on spiking neural networks. In: Proceedings of EPFL LATSIS Symposium 2006, Dynamical Principles for Neuroscience and Intelligent Biomimetic Devices, pp. 121–122 (2006)
  * [20] Russell, A., Orchard, G., Etienne-Cummings, R.: Configuring of spiking central pattern generator networks for bipedal walking using genetic algorithms. In: International Synopsium on Circuits and Systems, pp. 1525–1528. IEEE (2007)
  * [21] Rybak, I.A., Shevtsova, N.A., Lafreniere-Roula, M., McCrea, D.A.: Modelling spinal circuitry involved in locomotor pattern generation: insights from deletions during fictive locomotion. The Journal of physiology 577(2), 617–639 (2006)
  * [22] Smolinski, T.G.: Automated Parameter Search in Small Network Central Pattern Generators, pp. 292–294. Springer New York, New York, NY (2015). DOI 10.1007/978-1-4614-6675-823.
  * [23] Walter, F., Röhrbein, F., Knoll, A.: Neuromorphic implementations of neurobiological learning algorithms for spiking neural networks. Neural Networks, 72, 152-167 (2015)