# Deep limits and cut-off phenomena for neural networks

Benny Avelin<EMAIL_ADDRESS>
Department of Mathematics
University of Uppsala
Box 256, 751 05 Uppsala, Sweden Anders Karlsson<EMAIL_ADDRESS>
Department of Mathematics
University of Uppsala
Box 256, 751 05 Uppsala, Sweden
Section de mathématiques
Université de Genève
Case Postale 64, 1211 Genève 4, Switzerland

###### Abstract

We consider dynamical and geometrical aspects of deep learning. For many
standard choices of layer maps we display semi-invariant metrics which
quantify differences between data or decision functions. This allows us, when
considering random layer maps and using non-commutative ergodic theorems, to
deduce that certain limits exist when letting the number of layers tend to
infinity. We also examine the random initialization of standard networks where
we observe a surprising cut-off phenomenon in terms of the number of layers,
the depth of the network. This could be a relevant parameter when choosing an
appropriate number of layers for a given learning task, or for selecting a
good initialization procedure. More generally, we hope that the notions and
results in this paper can provide a framework, in particular a geometric one,
for a part of the theoretical understanding of deep neural networks.

Keywords: Deep limits, neural network, deep learning, ergodic theory, metric
geometry

##  1 Introduction

In this paper, we develop a geometric toolkit which we propose as a means to
study neural networks, in particular as the depth tends to infinity. Viewed in
this sense, we can deduce properties of the limit of neural networks as the
number of layers go to infinity provided that the layers preserve certain
distance functions (metrics).

As a starting point we will consider neural networks with random layers, which
can occur from dropout Srivastava et al. (2014), Bayesian neural networks Neal
(2012), neural networks with noise (Neural SDE) Liu et al. (2019), or simply
random initialization. The assumptions on the dependence between subsequent
layers is weak and we only assume stationarity.

Our analysis shows that under the assumption of stationarity, if one can find
a metric space for which the “layer transformations” of the neural networks is
non-expansive, then the limit and its growth rate can be described using
powerful tools from ergodic theory.

###  1.1 Background

As is by now well known, certain deep networks perform better than their
shallow counterpart, see for instance He et al. (2016). Fairly recently Belkin
et al. (2019) observed a phenomenon later dubbed “Deep Double Descent” in
Nakkiran et al. (2020). The deep double descent means that for deep networks,
after a certain width threshold the generalization properties becomes better
and better even though the class of networks becomes increasingly complex,
Hornik et al. (1989). However, as mentioned above, also deeper networks seems
better in terms of generalization, which suggests there is a regularizing
effect of depth under certain conditions, similar to deep double descent.

The wide limit of neural networks is fairly well studied, see for instance
Neal (1996); Mei et al. (2018); Jacot et al. (2018); Rotskoff and Vanden-
Eijnden (2018). However the deep limit is not a particularly well defined
concept and there are many different ways to view it. One of the more
practically successful ones are the Neural ODEs, introduced by Chen et al.
(2018), which can be seen as a deep limit of residual networks (Thorpe and van
Gennip (2018); Avelin and Nyström (2021)). The discrete model for these
continuous neural networks can be formulated as

$\displaystyle x_{t_{i+1}}=x_{t_{i}}+\frac{1}{n}T_{i}(x_{t_{i}}),\quad
i=1,\ldots,n$ (1)

where $T_{i}$ represents a layer in the neural network. In the case of Neural
SDEs (Liu et al. (2019); Tzen and Raginsky (2019)), or in the Bayesian
framework (see for instance Neal (2012)) we can view each layer as being
random. Furthermore, the special case of i.i.d. random layers is present in
the random initialization of the network, and is in fact a very important
aspect to understand when it comes to training neural networks, Sutskever et
al. (2013).

A key observation in the above formulation is that the discrete form 1
represents an approximation of the ODE

$\displaystyle\frac{\partial x(t)}{\partial t}=T_{t}(x_{t}),\quad t\in[0,1],$
(2)

i.e. the discrete system is an approximation of a fixed time horizon ODE, with
a time-step of size $1/n$. Another point of view is to consider a fixed time-
step and consider the behavior of the system as $n\to\infty$, i.e. of

$\displaystyle x_{t+1}=x_{t}+T_{t}(x_{t}),\quad t=1,\ldots,n.$ (3)

In a sense the networks in 3 represent a duality of thought compared to 1, in
that either we consider a deep neural network as consisting of many distinct
layers or we consider a single (time dependent) network being repeated in a
recurrent fashion, “layer” is represented by time. This connects deep neural
networks to the concept of a recurrent neural network and in a sense they are
the same, specifically neural ODEs and neural SDEs which are even trained
using recurrent back-propagation (real-time recurrent learning), see Chen et
al. (2018); Liu et al. (2019); Robinson and Fallside (1988); Rohwer (1990);
Tzen and Raginsky (2019); Williams and Zipser (1989).

In the context of Bayesian neural networks, there has recently been some
progress in establishing the deep limits of these as certain Gaussian
processes, see Agrawal et al. (2020); Dunlop et al. (2018); Duvenaud et al.
(2014).

According to E et al. (2020), at the continuous level many machine learning
models are the gradient flow of a reasonably nice functional, and they argue
that this is a reason for models such as 1 (ResNets) are numerically stable,
Hanin (2018); Schaefer et al. (2008). They also suggest that for 3 one should
expect trouble since there is no continuum limit. This is true to some degree,
but one should keep in mind that even a standard ResNet does not have the
scaling factor $1/n$ in front of $T_{n}$, thus one could argue that it is more
reasonable to consider the limit of fixed time-step dynamics. Of course this
may not have a limit but perhaps a certain rescaling does, for instance, one
could consider

$\displaystyle\lim_{n\to\infty}\frac{1}{n}\log(x_{n}),$

or

$\displaystyle\lim_{n\to\infty}\frac{1}{n}x_{n}.$

###  1.2 Our contribution

In this paper we take the viewpoint of 3 and we rephrase the update equation
as $x_{n+1}=T_{n}(x_{n})$. The problem is now one of discrete (possibly
chaotic) dynamical systems. The main contribution of our paper is that we
develop a framework to study deep neural networks from a geometric
perspective. Specifically it allows us to read out certain stability
properties whenever there exists a metric which is preserved by the network
layers. Depending on the metrics involved, it can tell us if the networks tend
to satisfy some regularity as we go deeper, even though in principle the
networks can become arbitrarily complex mappings (Lu et al. (2017)). This
serves as an indication as to when one would observe the “Deep Double Descent”
phenomenon with respect to depth.

Note that we study the network directly without worrying about how we obtained
said network. We consider stationary sequences of layers, which can cover for
instance a sequence correlated layers for which the correlation tends to zero
as we go deeper.

In the context of independent, identically distributed (i.i.d.) random layers,
which corresponds to the random initialization of the weights in the network,
we perform a few experiments where we observe a cut-off phenomenon. This means
that for a given type of neural network there is a certain number of layers
where the network behaves very differently. We call this its _cut-off depth_.
The cut-off phenomenon was first discussed in card shuffling by Aldous and
Diaconis (1986) who showed that seven shuffles are enough. It remains to
understand the full significance of the cut-off depth for deep learning.

Figure 1: The cut-off phenomenon for neural network mixing with TanH
activation. The figure describes the total variation (TV) distance to the
equilibrium measure as a function of the number of layers. The leftmost figure
is width 1 and the rightmost figure is width 2.

In summary, part of the rationale for our study is:

  * •

understanding how the function class evolves as the number of layers go to
infinity.

  * •

understanding if the discrete time-step model leads to a cut-off phenomenon?
That is, is there a threshold amount of layers for which the space of output
functions of the network increase significantly?

  * •

the randomness of our mappings $T_{i}$ is a way to model a host neural
networks, including standard regularization techniques, Bayesian networks,
other noise injected models or random initialization.

##  2 The dynamics of deep neural networks

Let $X\subseteq\mathbb{R}^{d}$ denote the space which contains the input
information as well as the intermediate data traveling through the hidden
layers, in the notation below, $x_{n}\in X$ for all $n$. It could be the full
vector space, or a subset such as the positive cone or a unit cube. Each layer
defines a transformation $T:X\rightarrow X$ typically of the form:

$T:x\mapsto\sigma(Wx+b),$

where $W$ is a matrix, called _weights_ , and $b$ a vector, called _bias
vector_. Note that this is just an example and we can in fact have $T$ being a
small network. The _activation function_ $\sigma$ is a non-linear function
$\sigma:\mathbb{R}\rightarrow\mathbb{R}$ that is fixed for the whole network
and is applied to each coordinate. Two standard choices for $\sigma(t)$ are
$\max\\{0,t\\},$ called the rectified linear unit (ReLU), and $\tanh(t)$
(TanH). The former has been observed to often work very well in practice and
the latter has the advantage of being a diffeomorphism
$\mathbb{R}\rightarrow(-1,1)$. Other common choices are
$\min\\{1,\max\\{0,t\\}\\}$ and the sigmoid or logistic function
$1/(1+e^{-t})$.

As mentioned in the introduction, in deep learning one uses several layers,
sometimes even up to a thousand. We denote these $n$ layer transformations,
$T_{1},...,T_{n}$. In order to gain some theoretical understanding for the
role of the number of layers, the depth, in neural network, we are interested
in what happens to the neural network when $n$ is large, or
$n\rightarrow\infty$. There are now in fact two possible dynamics to look at,
first, new layers are added at the end just before the final output, or
second, new layers are added just after the initial input. For a given initial
input $x_{0}\in X$, these two dynamics correspond to, respectively,

$\displaystyle x_{n}=T_{n}T_{n-1}...T_{1}x_{0}$ (4)

and

$\displaystyle x_{n}=T_{1}T_{2}...T_{n}x_{0}.$ (5)

We can view these dynamics as a representation of transfer learning. Where 4
corresponds to adding a new layer at the end, which is the standard way
transfer learning is used. On the other hand 5 corresponds to keeping the last
layers and inserting a new layer at the beginning, this can be seen as
transfer learning in the context of domain adaptation. Note that if we take
the maps randomly, or more precisely independent and identically distributed
(i.i.d.), then for each fixed $n$ and given $x_{0}$ the distributions of
$x_{n}$ are the same. But if we study the dynamics, the evolution $x_{n}$ of
an individual $x_{0}$ then these two dynamics behave differently because of
the non-commutativity of the layer transformations.

In addition to these two “input dynamics” we will also consider “output
dynamics”. At the last hidden layer one often has a _decision function_
$f:X\rightarrow Y$, where $Y$ is also a subset of some vector space. For
example, $f$ is in many cases an indicator function of some set or a smooth
approximation thereof.

Like the transpose of matrices, or more precisely the adjoint of operators
which transforms linear functionals instead of vectors, we can look at the
effect of the layers on the output function. We can thus let the layers and
dynamics transform the decision function to a function that is then directly
applied to the initial input. In other words, we are pulling back the decision
function to the input as it were. This is the dual dynamics. In formulas,

$(T^{*}f)(x):=f(Tx).$

Note that orders get reversed just like for the transpose of matrices:

$(T_{1}^{*}T_{2}^{*}...T_{n}^{*}f)(x_{0})=f(T_{n}T_{n-1}...T_{1}x_{0})$

and corresponding to the second dynamics above:

$(T_{n}^{*}T_{n-1}^{*}...T_{1}^{*}f)(x_{0})=f(T_{1}T_{2}...T_{n}x_{0}).$

##  3 Framework and strategy

We introduce a geometric viewpoint on neural networks. In several of the most
poplar network models in deep learning, we exhibit an associated metric space
on which the layer maps act as non-expansive maps.

Once we have this one could potentially use the contraction mapping principle,
in the version of a sequence of maps, which composed has a summable
contraction constant. We refer to the review by Diaconis and Freedman (1999)
for more information.

But often in our context such strong contraction property is not available,
but then we instead have the non-commutative ergodic theorem in Gouëzel and
Karlsson (2020), recalled below as Theorem 1, as a main tool.

In order to probe the neural network to understand a bit better the role of
the number of layers, we apply these ergodic theorems to stationary sequences
of layer maps. In experiments we observe a cut-off phenomenon, see Section 7.

We will in this paper analyze the following two dynamical systems:

$x_{n}=T_{1}T_{2}...T_{n}x_{0}$

or

$f_{n}:=T_{1}^{*}T_{2}^{*}...T_{n}^{*}f.$

These cover the two dynamics above: adding layers at the beginning or at the
end, respectively. It is this order that corresponds to random walks, where
each step $y_{n}$ is not far from $y_{n-1}$. Expressed differently, these
orders will make it possible to extract some coherent asymptotic behavior of
$x_{n}$ in the first case and $f_{n}$ in the second. (For certain quantities
such as the probability distribution or the basic growth the order does not
matter. This will be seen later).

##  4 Ergodic theorems

To better understand where the tools of ergodic theory come from let us recall
the classical law of large numbers (LLN). It asserts that for i.i.d. random
variables $X_{i}$ with $E[|X_{i}|]<\infty$,

$\displaystyle\frac{1}{n}(X_{1}+X_{2}+\ldots+X_{n})\to
E[X_{1}]\quad\text{almost surely.}$

One can wonder if there is a similar law when the random variables are not
commuting, for example the product of randomly selected matrices. Notice that
in the non-commutative case it is not obvious how to form an average, and one
complication is clear _a priori_ : the limits typically will depend on the
order of the maps, in contrast to the classical LLN. The results concerning
such non-commutative ergodic theorems are of two types, subadditive ergodic
theorems and multiplicative ergodic theorems. Let us begin by describing
Kingman’s subadditive ergodic theorem Kingman (1973) as it is a generalization
of the LLN to operations that are subadditive. As a simple special case (the
Furstenberg-Kesten theorem) let us consider a sequence of i.i.d. random
matrices $A_{i}\in\mathbb{R}^{N\times N}$ and the functions
$a(i,j)=\log||A_{i}\ldots A_{i+j}||$. This is subadditive in $i$, i.e.
$a(1,i+j)\leq a(1,i)+a(i,j)$. The subadditivity follows from the basic norm
inequality. In this case, Kingman’s subadditive ergodic theorem asserts that

$\displaystyle\lim_{n\to\infty}\frac{a(1,n)}{n}\to\ell.$

The limiting value is deterministic like in the LLN but on the other hand
there is no good formula for its value.

Multiplicative ergodic theorems are stronger in the sense that they say more
about the limit, loosely speaking they give an asymptotic direction of the
limit. The prototypical multiplicative ergodic theorem is Oseledets’ theorem
Oseledets (1968), which relates to random products of matrices.

The setting for these ergodic theorems and for us here is that of integrable
ergodic cocycles. This corresponds to stationary sequences in probability
theory language, a special case of which is the i.i.d with finite first moment
setting. Another special case is the mixing case where one has asymptotic
independence.

We will call such integrable ergodic cocycles simply _stationary sequences_
(with the integrability condition implicitly understood). Technically it means
that we have an underlying probability measure space $(M,\mu)$, $\mu(M)=1$,
and a measurable transformation $L:M\rightarrow M$ that preserves the measure
$\mu(L^{-1}A)=\mu(A)$ and is ergodic, which means that every $L$-invariant set
has measure either $0$ or $1$. Finally we have a measurable map from $M$ to a
set of layer maps, $m\mapsto T_{m}$, so that for the metric under
consideration, all distances involved are measurable and

$\displaystyle\int_{M}d(T_{m}x,x)d\mu(m)<\infty.$

Note that this condition is independent of $x$ since each map $T_{m}$ is non-
expansive in the metric.

Consider now a stationary sequence of random matrices
$T_{i}\in\mathbb{R}^{N\times N}$, then Oseledets’ theorem states that

$\displaystyle\lim_{n\to\infty}\frac{1}{n}\log\|T_{1}\ldots T_{n}x\|=\lambda$

the value of $\lambda$ is not random but depends on $x$, and there is at most
$N$ different values, which are called Lyapunov exponents. A good way of
understanding the Lyapunov exponents is to consider the special case where all
the $T_{i}$’s are the same matrix, in this case we just taking powers of this
matrix and the Lyapunov exponents are just the logarithm of the absolute value
of the eigenvalues.

Oseledets’ theorem only applies to linear maps (or the derivative cocycle of a
diffeomorphsism), but we need to analyze non-linear maps since we have an
activation function. To describe and understand the dynamics of such more
general settings we should define something quantitative, like a norm or a
metric. For example, to measure how close one decision function is to another,
or the distance between two information vectors. These norms and distances
should be preserved or at least not increase when a transformation is applied
to any two points. Specifically, Let $d$ denote a metric on either $X$ or some
space $L$ of functions $X\rightarrow Y$. We are interested in metrics that are
_semi-invariant_ , i.e. metrics that for a given map $U$ satisfies

$\displaystyle d(U(z),U(w))\leq d(z,w).$

for all $z$ and $w$ in $X$ or $L$, respectively. Correspondingly we will call
a mapping that has a semi-invariant metric, as a _non-expansive_ map with
respect to said metric.

We choose layer transformations $T_{k}$ at random (a stationary sequence). We
let $z_{n}$ denote either of the two random processes,

$T_{1}T_{2}...T_{n}x_{0}\,\,\textrm{or}\,\,T_{1}^{*}T_{2}^{*}...T_{n}^{*}f.$

Fix _ $d$_ a semi-invariant metric wiht respect to all the layer maps $T_{k}$.
Like in the matrix example above, it is easy to see that $d(z_{0},z_{n})$ is
then a subadditive process, see for example Karlsson and Ledrappier (2011);
Gouëzel and Karlsson (2020) (this is most clearly written with ergodic
theoretic formalism). Then by Kingman’s subadditive ergodic theorem Kingman
(1973) we have that

$\lim_{n\rightarrow\infty}\frac{1}{n}d(z_{0},z_{n})$

exists a.s. This can be viewed as the existence of a basic regularity or
growth. This holds for all the dynamics considered, including the reverse
orders.

Multiplicative ergodic theorems refine this convergence, in the way that it
predicts a directional behavior of $z_{n}$ (compare again with the matrix case
above). The precise statement, which in fact generalizes Oseledets’ theorem
above, is as follows:

###### Theorem 1

(Gouëzel and Karlsson (2020)) For any stationary sequence of maps as above,
with $z_{n}$ denoting the orbit, there exists a.s. a metric functional $h$
(that is _a priori_ random) such that

$\lim_{n\rightarrow\infty}-\frac{1}{n}h(z_{n})=\lim_{n\rightarrow\infty}\frac{1}{n}d(z_{n},z_{0}).$

In order to apply this general theorem we need to understand what the metric
functionals are for a given metric, see Section A.

In summary, the strategy we suggest is as follows:

  1. 1.

Given a selected type of layer maps, find a metric space on which the maps act
non-expansively.

  2. 2.

Determine the metric functionals of this space

  3. 3.

Apply the non-commutative ergodic theorem and interpret the result in the
given situation.

Moreover, we believe that already the first step, the metric setting, will
have other interests for deep learning, different from the application of
ergodic theorems.

##  5 Metric spaces

A metric space is a set with a distance function. In recent decades it has
been realized that significant geometrical arguments work in such a general
setting even without any differentiability. The subject is now often called
metric geometry and has begun to infiltrate other areas, such as computer
science.

Here we now turn to describing some metrics that are relevant in our context.
The most basic metric space is the euclidean space of some finite dimension.
Here the metric is given by $d(x,y)=||x-y||$ where the norm comes from a
scalar product.

Convex cones in vector spaces admit several useful choices of metrics with
important maps being non-expansive, we refer to the excellent book Lemmens and
Nussbaum (2012) for full information. We give some special cases of this here.
Given a finite dimensional real vector space, consider the set $X$ of vectors
having all its entries positive. Thus $X$ is a generalized first quadrant, a
convex cone. Consider the following expression

$f(x,w)=\log\left(\max_{i}\frac{w(i)}{x(i)}\right)$

in terms of the coordinates $x(i)$ and $w(i)$ of the vectors $x,w\in X$. Note
that $f$ is asymmetric in its arguments, so in order to build a metric, we
need to symmetrize it. It is also clearly not necessarily positive. There are
two options, the Thompson metric and the Hilbert projective metric
(“projective” refers to that it is a distance function between lines, while
the distance between two proportional vectors are easily seen to be 0. The
triangle inequality is not obvious, but we refer again to Lemmens and Nussbaum
(2012) for that.) The Thompson’s metric is defined as:

$d(x,w)=\log\left(\max\left\\{\max_{i}\frac{x(i)}{w(i)},\max_{i}\frac{w(i)}{x(i)}\right\\}\right)$

in terms of the coordinates $x(i)$ and $w(i)$ of the vectors $x,w\in X$. This
makes the cone $X$ into a metric space and its main feature is that any order-
preserving, subhomogeneous map of the cone into itself is a non-expansive map
in this metric, for more information see Section 6.1.

The Hilbert metric is instead

$d(x,w)=\log\left(\max_{i}\frac{x(i)}{w(i)}\cdot\max_{i}\frac{w(i)}{x(i)}\right).$

This is also obviously a symmetric expression. On the other hand it clearly is
$0$ on rays, $d(x,\lambda x)=0$, $\lambda>0$. More generally if we define the
equivalence relation $x\underset{X}{\equiv}y$ if there is a $\lambda>0$ such
that $\lambda x=y$, then $((X,\underset{X}{\equiv}),d)$ is a metric space.
Furthermore, order-preserving, homogeneous maps are non-expansive in this
metric.

In one possible approach to Oseledets’ theorem, one looks at the associated
action of the matrices on the space of positive scalar products on the
underlying vector space, see Karlsson and Ledrappier (2011) and references
therein. Since our maps are not linear we cannot do the same. We will instead
suggest to look at a much larger space, namely distance functions either on
$\mathbb{R}^{N}$ or the cube. To illustrate these ideas we fix the cube
$X=[-1,1]^{N}$ and use TanH which map the whole vector space diffeomorphically
onto the open cube. Let $M$ be the set of distance functions on $X$ bi-
Lipschitz equivalent to the original distance defined by a norm
$\left\|\cdot\right\|$. Here is a metric on $M$:

$D(d_{1},d_{2})=\log\left(\max\left\\{\sup_{x\neq
y}\frac{d_{2}(x,y)}{d_{1}(x,y)},\sup_{x\neq
y}\frac{d_{1}(x,y)}{d_{2}(x,y)}\right\\}\right).$

Notice that for two distance functions that are $K$-bi-Lipschitz to each
other, $0<D(d_{1},d_{2})<\infty$. For the optimal $K>1$ their distance would
be $\log K$. This function is clearly symmetric and $D(d_{1},d_{2})=0$ if and
only if $d_{1}=d_{2}.$ The triangle inequality is also satisfied because of
the obvious properties of sup and log, like for the Thompson metric above. So
$(M,D)$ is a metric space.

If T are diffeomorphisms then obviously this metric is invariant under $T^{*}$
which maps $d$ to $T^{*}d(x,y):=d(Tx,Ty)$, since $T$ just permutes the
underlying set, leaving the supremum invariant. Suppose T is merely injective
(otherwise there would occur a division by zero), then $T$ is non-expansive in
view of the inequality:

$\log\left(\sup_{x\neq
y}\frac{d_{2}(Tx,Ty)}{d_{1}(Tx,Ty)}\right)\leq\log\left(\sup_{x\neq
y}\frac{d_{2}(x,y)}{d_{1}(x,y)}\right)$

since $TX\subset X$.

From this perspective, it seems that injectivity is important, TanH and
invertible matrices gives (non-surjective) isometries of $(M,D)$. Other
activation functions would also be possible, but not ReLU.

In one dimensional dynamics, say in the study of diffeomorphisms $f$ of a
finite interval $I$, one finds the following quantitative measures of
distortion and distance. First out is

$\text{Var}(\log(\text{D}\,f))=\int_{I}\,\left|\frac{\text{D}^{2}\,f}{\text{D}\,f}\right|\,dx$

This has subadditive properties under composition of maps, see Navas (2018).
In particular if one takes a random composition and divides by n, this
converges a.e. to a deterministic value, called the _asymptotic variation_.

Another measure of distortion is

$D(f)=\sup_{x,y\in
I}\left|\log\left|\frac{f^{\prime}(x)}{f^{\prime}(y)}\right|\right|,$

see Deroin et al. (2007). $D$ is subadditive with respect to compositions
$f^{n}$ and symmetric with respect to $f$ and the inverse $f^{-1}$. We can
make this into a metric as follows, _the distortion metric_ :

$d(f,g)=\sup_{x,y\in
I}\left|\log\left|\frac{g^{\prime}(x)}{f^{\prime}(x)}\frac{f^{\prime}(y)}{g^{\prime}(y)}\right|\right|$

but one can also consider a Thompson version. Furthermore, Theorem 1 also
holds for asymmetric distances, see Gouëzel and Karlsson (2020), this allows
us to even consider just half of it, i.e taking away $y$.

We can extend the distortion metric to higher dimension by considering the
Jacobians instead of the derivatives. Let us consider diffeomorphisms
$f,g:\Omega\to\Omega$ in $\mathbb{R}^{N}$, then we can consider the Jacobian
matrix $J_{f}=\left\\{\frac{\partial f_{i}}{\partial x_{j}}\right\\}$ and the
Jacobian determinant $|J_{f}|$. Let us define the pseudo-metric

$\displaystyle
D(f,g)=\sup_{x,y\in\Omega}\left|\log\frac{|J_{f}(x)|}{|J_{g}(x)|}\frac{|J_{g}(y)|}{|J_{f}(y)|}\right|$

Diffeomorphisms $h:\Omega\to\Omega$ are isometries, i.e.

$\displaystyle d(f\circ h,g\circ h)$
$\displaystyle=\sup_{x,y\in\Omega}\left|\log\frac{|(J_{f}\circ
h)(x)J_{h}(x)|}{|(J_{g}\circ h)(x)J_{h}(x)|}\frac{|(J_{g}\circ
h)(y)J_{h}(y)|}{|(J_{f}\circ h)(y)J_{h}(y)|}\right|$
$\displaystyle=\sup_{x,y\in\Omega}\left|\log\frac{|(J_{f}\circ
h)(x)||J_{h}(x)|}{|(J_{g}\circ h)(x)||J_{h}(x)|}\frac{|(J_{g}\circ
h)(y)||J_{h}(y)|}{|(J_{f}\circ h)(y)||J_{h}(y)|}\right|$
$\displaystyle=\sup_{x,y\in\Omega}\left|\log\frac{|J_{f}(x)|}{|J_{g}(x)|}\frac{|J_{g}(y)|}{|J_{f}(y)|}\right|=d(f,g)$

The second to last follows from $\det(AB)=\det(A)\det(B)$ and the last step
follows from the fact that $h$ is a diffeomorphism.

In the case where $h:\Omega\to\Omega$ is not a diffeomorphism with non-
singular Jacobian, they are non-expansive in the view of

$\displaystyle d(f\circ h,g\circ h)$
$\displaystyle=\sup_{x,y\in\Omega}\left|\log\frac{|(J_{f}\circ
h)(x)||J_{h}(x)|}{|(J_{g}\circ h)(x)||J_{h}(x)|}\frac{|(J_{g}\circ
h)(y)||J_{h}(y)|}{|(J_{f}\circ h)(y)||J_{h}(y)|}\right|$
$\displaystyle\leq\sup_{x,y\in\Omega}\left|\log\frac{|J_{f}(x)|}{|J_{g}(x)|}\frac{|J_{g}(y)|}{|J_{f}(y)|}\right|=d(f,g).$

##  6 Main results

In this section we will use our previously outlined strategy to derive
conclusions about the deep limit of neural networks. We will be employing both
subadditive and multiplicative ergodic theorems to do so.

###  6.1 Positive models

We begin our exposition into explicit examples, by considering what we call
positive models. That is, layers that can only produce positive output. To be
specific, let us take for $X$ the cone of vectors in $\mathbb{R}^{N}$ with all
coordinates $\geq 0$. The layer maps $T(x)=\sigma(Wx+b)$ are such that $W$ is
a matrix with every entry $\geq 0$ and same for $b$. Finally $\sigma$ is an
activation function which is increasing and satisfies $\sigma(\lambda
x)\leq\lambda\sigma(x)$ for every $\lambda>0$ and $x\geq 0$, (note that this
implies that $\sigma(0)=0$). For example, ReLu, TanH, and the sigmoid.

Note the following properties:

  * •

$W$ preserves the cones, i.e. maps $X$ into $X.$ So does $b$ and finally also
$\sigma.$ Therefore $T:X\rightarrow X$. (With ReLU this is true without
assumptions on $W$ and $b$.)

  * •

In fact, more is true, if $x\leq y$ in the partial order defined by the cone
(i.e. all components of $x$ are smaller or equal to those of $y)$ then by the
positivity of $W$ and $b$ as well as the increasing property of $\sigma$, it
holds that $Tx\leq Ty$.

Such maps are called _order-preserving_.

###### Definition 2

A map $f:X\to X$ is called _subhomogeneous_ if $\lambda f(x)\leq f(\lambda x)$
for all $x\in X$ and $0<\lambda<1$.

###### Example 1

Let us consider the $1$-dimensional case, in this case $a,b,x\in\mathbb{R}$.
Let us consider the sigmoid $\sigma(x)=\frac{1}{1+e^{-x}}$ and prove that
$\sigma(ax+b)$ is subhomogeneous for any $b>-2$ if $a>0$.

First assume that $\lambda\in(0,1]$. Call $f_{\lambda}(x)=\lambda\sigma(ax+b)$
and $g_{\lambda}(x)=\sigma(a\lambda x+b)$. We will prove the strict inequality
$f_{\lambda}(x)<g_{\lambda}(x)$, assume now that for an $x_{0}$ there is a
$\lambda_{x_{0}}\in(0,1]$ such that
$f_{\lambda_{x_{0}}}(x_{0})=g_{\lambda_{x_{0}}}(x_{0})$ then let us
differentiate w.r.t $\lambda$ and see

$\displaystyle\frac{d}{d\lambda}f_{\lambda}(x_{0})\bigg{\lvert}_{\lambda=\lambda_{x_{0}}}$
$\displaystyle=\sigma(ax_{0}+b)=g_{1}(x_{0})$
$\displaystyle\frac{d}{d\lambda}g_{\lambda}(x_{0})\bigg{\lvert}_{\lambda=\lambda_{x_{0}}}$
$\displaystyle=ax_{0}\sigma^{\prime}(a\lambda_{x_{0}}x_{0}+b)=ax_{0}g_{\lambda_{x_{0}}}(x_{0})(1-\sigma(a\lambda_{x_{0}}x_{0}+b))$
$\displaystyle=ax_{0}f_{\lambda_{x_{0}}}(x_{0})(1-\sigma(a\lambda_{x_{0}}x_{0}+b))$
$\displaystyle=ax_{0}\lambda_{x_{0}}(1-\sigma(a\lambda_{x_{0}}x_{0}+b))g_{1}(x_{0})<g_{1}(x_{0})$

where in the above we have used that

$\displaystyle x(1-\sigma(x+b))<1$ (6)

which only holds for $b>-2$ and can be proven by straightforward computation.
Note that this implies that for $\lambda<\lambda_{x_{0}}$ we have
$f_{\lambda}(x_{0})<g_{\lambda}(x_{0})$. Now note that when $\lambda=1$ we
have $f(x)=g(x)$ for all $x$, by the above argument the strict inequality
carries to all $\lambda\in(0,1)$ for all $x$.

Notice that this example generalizes to positive matrices $a$ and vectors $x$
as long as the activation function is applied component-wise (as usual) and
the vector $b$ is greater than $-2$ in each component.

###### Example 2

Let us again consider the case $\sigma(x)=\frac{1}{1+e^{-x}}$ and prove that
given a positive matrix $W$ then for _any_ $b$ the mapping $Tx=\sigma(Wx+b)$
is order preserving. By the definition of order preserving we need to prove
that given $x\leq y$ we have $Tx\leq Ty$. Even though $Wx+b$ might not be
mapped into $X$ if $b$ is not a positive vector, we still have that
$(Wy+b)-(Wx+b)\in X$. This together with the monotonicity of $\sigma$ gives
that $T$ is order preserving on the positive cone.

Collecting the above, we may thus formulate:

###### Proposition 3

Given a matrix $W$ with positive entries, a vector $b$ with entries larger
than $-2$, and $\sigma(x)=\frac{1}{1+e^{-x}}$ the sigmoid activation function.
Then the associated layer map $T(x)=\sigma(Wx+b)$ is non-expansive with
respect to the Thompson metric on the standard positive cone.

###### Example 3

Consider now the ReLU activation function $\sigma(x)=\max\\{x,0\\}$. We affirm
that the mapping $Tx=\sigma(Wx+b)$ is subhomogeneous if $b\geq 0$ and $W$ is
arbitrary. To see this: for any $0<\lambda<1$, we have

$\displaystyle\lambda T(x)=\sigma(W\lambda x+\lambda B)\leq\sigma(W\lambda
x+B)=T(\lambda x).$

And as already remarked, if all entries of $W$ is positive, then $T$ is order-
preserving as well.

We thus have some interesting examples of order-preserving and subhomogeneous
layer transformations and since these properties are preserved under
composition we get a rich bank of examples. As was recalled above, the
positive cone admits a metric $d$, the Thompson metric, which is semi-
invariant under such maps, i.e.

$d(Tx,Tx^{\prime})\leq d(x,x^{\prime})$

for all $x,x^{\prime}\in X$.

In the following theorem we consider the mappings $T:X\rightarrow X$ to be
random in a stationary way and quite general, but we will keep the above
examples in mind. The point is here that in order for our dynamics to have a
well defined limit we consider the “layers” added to the beginning of the
network instead of at the end, i.e. we are interested in

$\displaystyle x_{n}=T_{1}T_{2}\ldots,T_{n}x_{0}.$

In conclusion, applying Theorem 1 we get:

###### Theorem 4

Let $X$ be the positive cone in $\mathbb{R}^{N}$ and let $T_{i}:X\to X$ be a
stationary sequence of maps that is order preserving and subhomogeneous. Let
$x_{n}=T_{1}T_{2}\ldots T_{n}x_{0}$, for a fixed $x_{0}\in X$. Then

$\lim_{n\rightarrow\infty}\sup_{i}|x_{n}(i)|^{1/n}=e^{\lambda}$

and there is a (random) coordinate $1\leq i_{0}\leq d$ such that

$\lim_{n\rightarrow\infty}|x_{n}(i_{0})|^{1/n}=e^{\lambda}.$

The latter statement excludes a certain spiraling inside the cone.

###  6.2 Unitary case

In the case of positive models that are order preserving and subhomogeneous,
we get two things, first of all there is essentially exponential growth of the
components of the trajectories, secondly we where able to determine the
direction of such a trajectory. The subhomogeneity and order preserving
transforms allowed us to find a metric which made the maps semi-invariant, if
we on the other hand loosen the restrictions and consider general transforms
we have the problem of finding good metrics. However if we instead restrict
the matrices $W$ to be unitary (spectral norm $1$) then even if we consider
fairly general norms we still get a non-expansive mapping. In this case we can
also deduce the specific form of the metric functional which gives us a very
concrete result.

To describe our situation let us begin by taking $X=\mathbb{R}^{N}$ together
with a norm $\|\cdot\|_{N}$ and layer maps $T(x)=W^{T}\sigma(Wx+b)$ such that
the corresponding operator norm $\left\|W\right\|_{N}\leq 1$, $b$ general and
$\sigma$ which satisfy
$\left\|\sigma(x)-\sigma(y)\right\|_{N}\leq\left\|x-y\right\|_{N}$ (i.e.
Lipschitz with constant $1$, or non-expansive) when applied to vectors
(component-wise as always). The point of having the layer transformations in
the form given above, is that it is a popular layer type that is used in for
instance ResNets (He et al. (2016)) and provides us with a layer that can span
the entire vector space. Let us now remark on the $1$-Lipschitz condition of
the activation function in our norm. Begin by noting that most used activation
functions are $1$-Lipschitz with respect to the standard absolute value, then
if we assume that the norm is monotone, i.e. $\|x\|_{N}\leq\|y\|_{N}$ if
$|x_{i}|\leq|y_{i}|$ for all $i$, the activation function becomes
$1$-Lipschitz in the norm $\|\cdot\|_{N}$. For the theorem below we also
assume the unit ball in the norm is a strictly convex set. We get by applying
Theorem 1 (in the form of Corollary 1.5 in Gouëzel and Karlsson (2020) ):

###### Theorem 5

Let $(X=\mathbb{R}^{N},\|\cdot\|_{N})$ be a normed vector space which has the
above monotonicity property and strictly convex unit ball. Consider a
stationary sequence of layer maps $T_{n}$ of the form
$T(x)=W^{T}\sigma(Wx+b)$, $\left\|W\right\|_{N}\leq 1$, $b\in X$, and $\sigma$
is $1$–Lipschitz when applied componentwise in
$(X=\mathbb{R}^{N},\|\cdot\|_{N})$. Then as $n\rightarrow\infty$ it holds that
a.s. there exists a vector $v$ such that

$\frac{1}{n}T_{1}T_{2}...T_{n}x_{0}\rightarrow v.$

The vector v is a priori random but independent of the initial data $x_{0}$.
The norm of $v$ is deterministic.

The above theorem is a consequence of the same theorem that gave rise to
Theorem 4, but in this case, since we have a norm we get that the metric
functional reduces to a dot-product (Proposition 8) and from this fact we can
read off the explicit convergence in the above theorem, see Karlsson (2019)
for more details.

###  6.3 An example of the reverse order

Let $\Omega=[-1,1]^{N}$ and use TanH which map the whole vector space
diffeomorphically onto the open cube. Let $X$ be the set of distance functions
on $\Omega$ bi-Lipschitz equivalent to the original distance $d_{0}$ defined
by a norm $\left\|\cdot\right\|$. The metric $D$ on $X$ is:

$D(d_{1},d_{2})=\log\left(\max\left\\{\sup_{x\neq
y}\frac{d_{2}(x,y)}{d_{1}(x,y)},\sup_{x\neq
y}\frac{d_{1}(x,y)}{d_{2}(x,y)}\right\\}\right).$

which was discussed in Section 5, see also Section A.3 for more details about
the corresponding metric functionals.

We consider maps $T$ of the usual type except for insisting on that the
matrices $W$ are invertible and also having fixed the activation function TanH
so that $T:\Omega\rightarrow\Omega$ and which will then be a diffeomorphism,
and as remarked in Section 5, leaving the distance of the metric space $(X,D)$
invariant for the induced action $T^{*}$.

We consider the reverse dynamics $T_{n}T_{n-1}...T_{0}x_{0}$, and shift the
point of view by instead studying

$d_{n}:=T_{1}^{*}T_{2}^{*}...T_{n}^{*}d_{0}.$

In the above we see that the maps are “inserted” just before $d_{0}$, which is
again the order which corresponds to random walks.

Thanks to the subadditive ergodic theorem and the other theorems in Gouëzel
and Karlsson (2020) $T_{1}^{*}T_{2}^{*}...T_{n}^{*}d_{0}$ will have some
regular behavior when $n\rightarrow\infty$ especially in terms of metric
functionals on the space $X$. We have the following result proven in Section
A.3.

###### Theorem 6

Under the above assumptions there is a number $\lambda$ so that

$\lim_{n\rightarrow\infty}\left(\sup_{x\neq
y}\frac{\|T_{n}T_{n-1}...T_{1}x-T_{n}T_{n-1}...T_{1}y\|}{\|x-y\|}\right)^{1/n}=e^{\lambda}$

Moreover, in case $\lambda>0$ there exists a point $x\in\Omega$ and a sequence
$z_{i}=(x_{i},y_{i})\in\\{(x,y):x,y\in\Omega,x\neq y\\}$ such that
$z_{i}\to(x,x)$ and for any $\varepsilon>0$ there is a number $N$ so that for
$n>N$

$\displaystyle\frac{\|T_{n}\ldots T_{1}x_{i}-T_{n}\ldots
T_{1}y_{i}\|}{\|x_{i}-y_{i}\|}\geq e^{(\lambda-\varepsilon)n},$

for all sufficiently large $i$.

The second assertion means in words that there is a random point $x\in\Omega$
and points approaching $x$ which the maps $T_{n}T_{n-1}...T_{1}$ separate with
maximal exponential rate.

###  6.4 Dynamics of decision functions

Let $\Omega\subset\mathbb{R}^{N}$ be a compact set. We consider the dynamics
$T_{n}T_{n-1}...T_{1}x$ for $x\in\Omega$ but shift the point of view to study
instead

$f_{n}:=T_{1}^{*}T_{2}^{*}...T_{n}^{*}f,$

where $f$ is the original decision function defined on $\Omega$. There are
several possibilities here, especially with TanH and $\Omega=[-1,1]^{N}$, but
we keep it simple and general to illustrate what can be shown. We assume that
$f$ and all layer maps are diffeomorphisms $\Omega\rightarrow\Omega$.

The maps are chosen in a stationary way as before. They preserve our Jacobi
distortion metric $D$, from Section 5. It is then a standard fact that
$d(f,f_{n})$ is a subadditive process. The subadditive ergodic theorem applies
and gives since $J_{f}$ is bounded and bounded away from $0$:

###### Theorem 7

In this situation there is a well-defined exponential growth rate $\lambda$ of
the distortion of the decision functions $f_{n}(x)=f(T_{n}T_{n-1}...T_{1}x)$,
more precisely,

$\lim_{n\rightarrow\infty}\sup_{x,y}|J_{f_{n}}(x)J_{f_{n}}(y)^{-1}|^{1/n}=e^{\lambda}.$

Via the theory of random dynamical systems as exposed in Arnold (1995), there
is an approach to the last theorem via Oseledets’ theorem, see also Li (2018)
for comparison.

##  7 Short time (layer) behavior for random weight initialization

When training deep neural networks an important aspect is how to do the random
weight initialization such that we get a network that can actually be trained.
In this section we will explore this concept a bit as it relates to our random
layer transformations in this paper. However, in contrast to the previous
theory we will study the “short time” behavior.

We will consider neural networks of the following simple type

$\displaystyle X_{i+1}=\sigma(W_{i}\cdot X_{i}),\quad i=0,\ldots$

where $W_{i}$ is i.i.d. from some distribution and $X_{0}$ is some starting
point of the network. If we consider this as a Markov chain, i.e. for a fixed
point we would take random weights at each step and consider the distribution
of the output for a fixed $X_{0}$.

We will view this from the context of mixing and as such we need to frame our
discussion with some terminology.

Consider a Markov chain on a domain of size $n$ (could be a graph or group for
instance) and let $P_{n}^{t}(x,\cdot)$ be the distribution of the Markov chain
started at $x$ at time $t$. Suppose that the Markov chain has a stationary
distribution $\pi_{n}$. It is well known that if the Markov chain is
irreducible and aperiodic we get exponential convergence, i.e. we get

$\displaystyle\max_{x}\|P_{n}^{t}(x,\cdot)-\pi_{n}\|_{TV}\leq e^{-ct}$

for some constant $c>0$. Here $\|\cdot\|_{TV}$ is the total variation (TV)
distance, which is defined for measures $\mu,\nu$ on $\Omega$ (finite state-
space $\Omega$)

$\displaystyle\|\mu-\nu\|_{TV}=\sum_{x\in\Omega}|\mu(x)-\nu(y)|.$

If we define

$\displaystyle d_{n}(t):=\max_{x}\|P_{n}^{t}(x,\cdot)-\pi_{n}\|_{TV},$

then we say that the Markov chain exhibits a cut-off at $t_{n}$ with window
$w_{n}$ if $w_{n}=o(t_{n})$, and

$\displaystyle\lim_{\alpha\to\infty}\liminf_{n\to\infty}d_{n}(t_{n}-\alpha
w_{n})$ $\displaystyle=1,$
$\displaystyle\lim_{\alpha\to\infty}\limsup_{n\to\infty}d_{n}(t_{n}+\alpha
w_{n})$ $\displaystyle=0.$

Intuitively speaking, this means that $d_{n}$ is close to $1$ for times just
below $t_{n}$ and $d_{n}$ is close to 0 for times slightly larger than
$t_{n}$, at least for large values of $n$. Note: the mixing time for the
Markov chain will be close to $t_{n}$ for large $n$. The mixing time for a
Markov chain is defined as

$\displaystyle t_{mix}(\epsilon)$
$\displaystyle:=\min\\{t:d(t)\leq\epsilon\\}$ $\displaystyle t_{mix}$
$\displaystyle:=t_{mix}(1/4).$

The cut-off phenomenon is stronger than the concept of fast mixing as it is a
double sided property.

###  7.1 Neural network induced Markov chains

In the examples that we will simulate below, the limiting distribution is
actually the point-mass at 0, and to make the total variation distance easier
to define we work with finite precision, which makes the state-space finite.
Note that if we have two measures $\mu$ and $\nu$ on a finite set $\Omega$
then,

$\displaystyle\|\mu-\nu\|_{TV}=\sum_{x\in\Omega}|\mu(x)-\nu(y)|$

i.e. the total variation distance is equivalent to the $L^{1}$ distance of the
densities, which is easier to compute. We now come to our first example,
namely a fully connected neural network with TanH activation and revisit the
study done in Glorot and Bengio (2010).

Consider the following simple Markov chain of neural network type (with
heuristic initialization, see Glorot and Bengio (2010))

$\displaystyle X_{i+1}=\tanh(W_{i}\cdot X_{i}),\quad i=0,\ldots$ (7)

where $X_{0}\neq 0,X_{0}\in\mathbb{R}^{N}$, where $W_{i}$ are i.i.d.
$W_{i}\sim\text{unif}\left(\left[-\frac{1}{\sqrt{N}},\frac{1}{\sqrt{N}}\right]^{N}\right)$,
and the TanH is applied componentwise (as is customary).

Above we can think of $N$ as the “size” of the Markov chain in the sense
above. The result of the simulation can be found in Fig. 2, where we worked
with a finite precision of $0.001$ and measure the total variation distance to
the point-mass at $0$. With the heuristic scaling factor introduced i.e. that
$W_{i}\sim\text{unif}\left(\left[-\frac{1}{\sqrt{N}},\frac{1}{\sqrt{N}}\right]^{N}\right)$
we see that they exhibit a cut-off at pretty much the same level. The cut-off
implies that the behavior of this random initialization is markedly different
for a layer count of around $3$ to layer counts above $10$.

It seems that this phenomenon occurs even for asymmetric activation functions,
like the Sigmoid-weighted Linear Unit (SiLU or Swish) (Ramachandran et al.
(2017); Elfwing et al. (2018))

$\displaystyle\sigma(x)=\frac{x}{1+e^{-x}}$ (8)

with the same setup as in 7 with the above activation, see Fig. 3. It does not
however seem to occur for non-smooth activations, like ReLU.

Figure 2: The cut-off phenomenon for Neural network mixing. The leftmost
figure is width 1 and the rightmost figure is width 2. Mixing times are 9 and
11 respectively. Figure 3: The cut-off phenomenon for Neural network mixing.
The figure corresponds to activation function 8. Mixing time is 7.

##  8 Conclusion and outlook

We have shown that some aspect of the understanding of deep learning networks
have a very natural dynamical interpretation, where recent ergodic theorems
can be applied. Indeed the tools are quite general and therefore there is a
wealth of possibilities. In other mathematical contexts this versatility has
already been demonstrated.

In the context of deep learning, one should translate the meaning of the
metrics and their functionals (say in terms of notions of complexity). This we
have achieved for several choices of metrics.

A question that arises from the cut-off phenomenon that we demonstrate
experimentally is its relevance for training of the network, that is, what
difference it makes in practice from choosing fewer layers than the cut-off
depth vis-à-vis choosing more layers. In fact, the fast convergence towards
the point mass at $0$ hints that the last layers in a deep network will have
an activation close to $0$, which is the linear regime of the activation
function, implying that the deeper layers are basically linear mappings.
Furthermore, in Glorot and Bengio (2010) they considered instead the
normalization which in our case becomes $\frac{\sqrt{3}}{\sqrt{N}}$ (as the
input is the same as the output dimension) which actually only delays the cut-
off to higher layer counts, it is however still there. One could speculate
that the variance of the initialization can as such be used to control how
nonlinear the initialized neural network is.

We believe that the metric setting will in forthcoming work have other
interests in deep learning, not only from the application of the non-
commutative ergodic theorem. This could ultimately inform the choice of best
design of the neural network for a given practical task.

##  A Metric functionals

The mathematical abstraction of distance is that of a metric space $X$ which
is a set with a distance function $d(x,y)$ that is symmetric, positive and
zero if and only if $x=y$. Moreover there is the fundamental triangle
inequality: the distance between $x$ and $y$ cannot be larger than the sum of
the distances from $x$ to $z$ and from $z$ to $y$, for any point $z$.
Sometimes it is useful and natural to relax the condition of symmetry.

We will now define the metric space analogs of linear functionals, affine
hyperplanes and half spaces. Let $(X,d)$ denote a metric space, fix $x_{0}\in
X$. Let $F(X,\mathbb{R})$ be the space of continuous functions
$X\rightarrow\mathbb{R}$ equipped with the topology of pointwise convergence.
We define the continuous injection

$\displaystyle\Phi:X\hookrightarrow F(X,\mathbb{R})$ $\displaystyle x\mapsto
h_{x}(\cdot):=d(\cdot,x)-d(x_{0},x)$

The functions $h_{x}$ are all non-expansive with respect to $d$ and vanish at
$x_{0}$. The image $\Phi(X)$ can be identified with a subset of a product of
compact intervals, which is compact by Tychonoff’s theorem. The closure of the
image $\overline{\Phi(X)}$ will therefore be compact (similar to the
compactness in the weak topology of functional analysis). See for example
Gaubert and Vigeral (2012) or Karlsson (2019) for details. We will call the
elements in this compact space _metric functionals_. In particular, to each
point $x$ there is the corresponding metric functional $h_{x}$. For metric
functionals that are genuine limits, their level-sets are called horospheres,
and sublevel sets are called horoballs. These two concepts are the metric
analogs of affine hyperplanes and half-spaces in the linear vector space
setting.

The first metric space to look at is the the finite dimensional euclidean
space, where the metric is given by $d(x,y)=||x-y||$ where the norm comes from
a scalar product. In this case the metric functionals are up to a constant:
the distance to a point in $X$ (in which case sublevel sets are balls) or
linear functionals of norm 1 (in which case sublevel sets are half-spaces).
Horospheres and horoballs are affine hyperplanes and half-spaces respectively.

###  A.1 Metric functionals in the case of smooth norms

Below we provide a characterization of the metric functionals in case of
norms. For another proof for $p$-norms see Gutiérrez (2019).

###### Proposition 8

Let $\|\cdot\|_{a}$ be a norm on $\mathbb{R}^{N}$ that is $C^{2}$ as a
function on the corresponding unit sphere. Consider the function

$\displaystyle h_{y^{n}}(x)=\|y^{n}-x\|_{a}-\|y^{n}\|_{a}$

where $\|y^{n}\|_{a}\to\infty$ as $n\rightarrow\infty$. Then there is a
subsequence and a vector $w$ with $\|w\|_{a}=1$ such that

$\displaystyle\lim_{n\to\infty}h_{y^{n}}(x)=-x\cdot\left(\nabla\|\cdot\|_{a}\bigg{\rvert}_{w}\right).$

Proof  Begin by first noting that $w^{n}=y^{n}/\|y^{n}\|_{a}$ is a vector on
the unit sphere and as such there is a subsequence of $\\{w^{n}\\}$ converging
to a vector $w$ s.t. $\|w\|_{a}=1$. We will dispense with notation for
subsequences for simplicity. Let us rewrite

$\displaystyle
h_{y^{n}}(x)=\|y^{n}-x\|_{a}-\|y^{n}\|_{a}=\frac{\left\|w^{n}-\frac{x}{\|y^{n}\|_{a}}\right\|_{a}-\|w^{n}\|_{a}}{1/\|y^{n}\|_{a}}$

Taylor expanding the norm $\|\cdot\|_{a}$ around $w^{n}$

$\displaystyle\|z\|_{a}=\|w^{n}\|_{a}+(z-w^{n})\cdot\nabla\|\cdot\|_{a}\rvert_{w}+f_{n}(z)\|z-w^{n}\|_{2}$

where $f_{n}(z)\to 0$ as $z\to w^{n}$, and $\|\cdot\|_{2}$ is the standard
Euclidean norm. For simplicity let us relabel our sequence such that
$\|y^{n}\|_{a}=n$, then

$\displaystyle\frac{\left\|w^{n}-\frac{x}{n}\right\|_{a}-\|w^{n}\|_{a}}{1/n}$
$\displaystyle=-x\cdot\left(\nabla\|\cdot\|_{a}\bigg{\rvert}_{w}\right)+\frac{f_{n}(x/n)\|x/n\|_{2}}{1/n}$
$\displaystyle=-x\cdot\left(\nabla\|\cdot\|_{a}\bigg{\rvert}_{w}\right)+f_{n}(x/n)\|x\|_{2}$

Together with a diagonal argument it is now clear that the $C^{2}$ regularity
of $\|\cdot\|_{a}$ gives us the limit as $n\to\infty$ (for a subsequence)

$\displaystyle\lim_{n\to\infty}h_{y^{n}}(x)=-x\cdot\left(\nabla\|\cdot\|_{a}\bigg{\rvert}_{w}\right).$

###  A.2 Metric functionals for the Thompson metric

In the following we make the identification of a vector
$x\in\mathbb{R}_{+}^{m}$ with a positive function
$x:I=\\{1,\ldots,m\\}\to\mathbb{R}_{+}$. For example, $\mathbbm{1}(i)=1$ for
all $i$. Consider the following

$\displaystyle d_{1}(x,y)=\log\sup_{I}\left(\frac{x}{y}\right)$

we wish to derive the limit of

$\displaystyle d_{1}(x,y_{n})-d_{1}(\mathbbm{1},y_{n})$

where $d(\mathbbm{1},y_{n})\to\infty$ as $n\to\infty$.

Consider the following

$\displaystyle d_{1}(x,y_{n})-d_{1}(\mathbbm{1},y_{n})$
$\displaystyle=\sup_{I}\left(\log(x)+\log(\mathbbm{1}/y_{n})-\sup_{I}\log(\mathbbm{1}/y_{n})\right)$
$\displaystyle=\sup_{I}\left(\log(x)+\log\left(\frac{\mathbbm{1}/y_{n}}{e^{\sup_{I}\log(\mathbbm{1}/y_{n})}}\right)\right)$

the function

$\displaystyle
u_{n}=\frac{\mathbbm{1}/y_{n}}{e^{\sup_{I}\log(\mathbbm{1}/y_{n})}}$

satisfies $0<u_{n}\leq 1$ and $\sup_{I}u_{n}=1$. Therefore there is a
subsequence of $u_{n}$ that converges (since $I$ is finite) to $u$ and we have
that the same subsequence satisfies

$\displaystyle
d_{1}(x,y_{n})-d_{1}(\mathbbm{1},y_{n})\to\sup_{I}\left(\log(x)+\log(u)\right)$

The other part of the Thompson metric satisfies a similar relation, i.e. let

$\displaystyle d_{2}(x,y)=\log\sup_{I}\left(\frac{y}{x}\right)$

we wish to derive the limit of

$\displaystyle d_{2}(x,y_{n})-d_{2}(\mathbbm{1},y_{n})$

where $d(\mathbbm{1},y_{n})\to\infty$ as $n\to\infty$.

Consider the following

$\displaystyle
d_{2}(x,y_{n})-d_{2}(\mathbbm{1},y_{n})=\sup_{I}\left(\log(1/x)+\log\left(\frac{y_{n}}{e^{\sup_{I}\log(y_{n})}}\right)\right)$

the function

$\displaystyle v_{n}=\frac{y_{n}}{e^{\sup_{I}\log(y_{n})}}$

satisfies $0<v_{n}\leq 1$ and $\sup_{I}v_{n}=1$. Therefore there is a
subsequence of $v_{n}$ that converges to $v$ (since $I$ is finite) and we have
that the same subsequence satisfies

$\displaystyle
d_{2}(x,y_{n})-d_{2}(\mathbbm{1},y_{n})\to\sup_{I}\left(\log(1/x)+\log(v)\right)$

This discussion can be compared with a similar one in Gutiérrez (2019). We
summarize everything in the following proposition:

###### Proposition 9

The metric functionals of the either half of the (or the full metric) Thompson
metric that arise as limits (also called horofunctions) are given as follows:
For $d_{1}$ we get that there exists a non-zero function $u:I\to[0,1]$ such
that

$\displaystyle h_{u}(x)=\sup_{I}\log(xu)$

and $\sup_{I}u=1$. For $d_{2}$ we get the existence of a non-zero function
$v:I\to[0,1]$ such that

$\displaystyle h_{v}(x)=\sup_{I}\log(v/x)$

and $\sup_{I}v=1$. In the case of $d=\max\\{d_{1},d_{2}\\}$ then we have

$\displaystyle
h_{u,v}(x)=\max\\{\sup_{I}\left(\log(x)+\log(u)\right),\sup_{I}\left(\log(1/x)+\log(v)\right\\},$

where $uv=0$ and $\sup_{I}\max\\{u,v\\}=1$.

####  A.2.1 Extension of the ideas to infinite dimensional spaces

Identifying the metric functionals in cases of infinite dimensional spaces is
more subtle, see for example Gutiérrez (2020). Especially for the $L^{\infty}$
space. We will take a look at a special case when we can actually determine
the boundary for the Thompson version of the distortion metric.

$\displaystyle d(f,g)=\max\left\\{\sup_{x\in
I}\log\left|\frac{g^{\prime}(x)}{f^{\prime}(x)}\right|,\sup_{x\in
I}\log\left|\frac{f^{\prime}(x)}{g^{\prime}(x)}\right|\right\\}$

We choose again as basepoint in our function space the identity function
$\mathbbm{I}(x)=x$. Note that

$\displaystyle d(\mathbbm{I},g)=\|\log(|g^{\prime}|)\|_{\infty}.$

###### Lemma 10

Consider the Thompson version of the distortion metrics used above,

$\displaystyle d(f,g)=\max\left\\{\sup_{x\in
I}\log\left|\frac{g^{\prime}(x)}{f^{\prime}(x)}\right|,\sup_{x\in
I}\log\left|\frac{f^{\prime}(x)}{g^{\prime}(x)}\right|\right\\}$

then if $g_{n}$ is s sequence of $C^{2}$ functions such that
$d(\mathbbm{I},g_{n})\to\infty$, satisfying the following differential
relation

$\displaystyle\sup_{I}|g_{n}^{\prime\prime}|\leq
C_{0}\sup_{I}|g_{n}^{\prime}|,$

for some $C_{0}>0$ independent of $n$, then there exists functions $u,v$ such
that there is a subsequence that converges as

$\displaystyle(d(f,g_{n})-d(\mathbbm{I},g_{n}))\to\sup_{I}\max\left\\{\log(|f^{\prime}|u),\log\left(\frac{1}{|f^{\prime}|}v\right)\right\\}.$

Furthermore the functions $u,v$ are Lipschitz with constant $C_{0}$ such that
$\sup_{I}\max\\{u,v\\}=1$.

Proof  First note that

$\displaystyle d(f,g_{n})-d(\mathbbm{I},g_{n})=\max\left\\{\sup_{x\in
I}\log\left|\frac{g_{n}^{\prime}(x)}{e^{d(\mathbbm{I},g_{n})}}\frac{1}{f^{\prime}(x)}\right|,\sup_{x\in
I}\log\left|f^{\prime}(x)\frac{1/g_{n}^{\prime}(x)}{e^{d(\mathbbm{I},g_{n})}}\right|\right\\}.$

Consider

$\displaystyle u_{n}$
$\displaystyle=\frac{g_{n}^{\prime}}{\exp(d(\mathbbm{I},g_{n}))},$
$\displaystyle v_{n}$
$\displaystyle=\frac{(1/g_{n}^{\prime})}{\exp(d(\mathbbm{I},g_{n}))},$

then

$\displaystyle|u_{n}|,|v_{n}|\leq\frac{e^{\pm\log(|g_{n}^{\prime}|)}}{\exp(d(\mathbbm{I},g_{n}))}\leq
1.$

This implies that

$\displaystyle|u_{n}^{\prime}|$
$\displaystyle=\frac{|g_{n}^{\prime\prime}|}{\exp(d(\mathbbm{I},g_{n}))}\leq\frac{C_{0}\sup_{I}|g_{n}^{\prime}|}{\exp(d(\mathbbm{I},g_{n}))}=C_{0}\sup_{I}|u_{n}|\leq
C_{0}$ $\displaystyle|v_{n}^{\prime}|$
$\displaystyle=\frac{\left|\frac{g_{n}^{\prime\prime}}{(g_{n}^{\prime})^{2}}\right|}{\exp(d(\mathbbm{I},g_{n}))}\leq\frac{C_{0}\frac{1}{\sup_{I}|g_{n}^{\prime}|}}{\exp(d(\mathbbm{I},g_{n}))}=C_{0}\sup_{I}|v_{n}|\leq
C_{0}.$

From this we get that there is a subsequence of $u_{n},v_{n}$ that converges
uniformly to $u_{\infty},v_{\infty}$ which are Lipschitz with constant
$C_{0}$. Furthermore, $\sup_{I}\max\\{u_{\infty},v_{\infty}\\}=1$.

###### Remark 11

Actually if we see the above proof, then the only thing we need is to make
sure that $u_{n}^{\prime},v_{n}^{\prime}$ is uniformly continuous. This
follows for instance if the modulus of continuity $\omega$ of $g^{\prime}$ is
bounded by $|g^{\prime}|$.

###  A.3 A Thompson metric for distance functions: Proof of Theorem 6

Let $\Omega$ be a compact subset of a finite dimensional vector space with a
norm $\|\cdot\|$. Let us consider the space $X$ of metrics on $\Omega$ which
are bi-Lipschitz equivalent to $d_{0}(x,y):=\|x-y\|$. Specifically this means
that $d\in X$ iff there exists a constant $C>1$ such that

$\displaystyle\frac{1}{C}\|x-y\|\leq d(x,y)\leq C\|x-y\|,\quad\forall
x,y\in\Omega.$

On the space $X$ we can consider a Thompson type metric, defined as

$\displaystyle D(d_{1},d_{2})=\log\left(\max\left\\{\sup_{x\neq
y}\frac{d_{1}(x,y)}{d_{2}(x,y)},\sup_{x\neq
y}\frac{d_{2}(x,y)}{d_{1}(x,y)}\right\\}\right),\quad d_{1},d_{2}\in X.$

To understand this metric, note that

$\displaystyle e^{-D(d_{1},d_{2})}d_{1}(x,y)\leq d_{2}(x,y)\leq
e^{D(d_{1},d_{2})}d_{1}(x,y).$

It is now easy to see that $X$ is complete under the metric $D$. Consider the
mapping $\mathcal{F}:X\to L^{\infty}(\Omega\times\Omega)\cap
C(\Omega\times\Omega\setminus\\{x=y\\})$, defined as

$\displaystyle\mathcal{F}(d)(x,y)=\log(d(x,y))-\log(\|x-y\|),$

denote $Y=\mathcal{F}(X)$. From the bi-Lipschitz condition we see that
$\sup_{x,y}|\mathcal{F}(d)(x,y)|\leq\infty$, furthermore

$\displaystyle\|\mathcal{F}d_{1}-\mathcal{F}d_{2}\|_{L^{\infty}(\Omega\times\Omega)}=D(d_{1},d_{2}).$

As such, we see that the mapping $\mathcal{F}$ is an isometric mapping of $X$
into $L^{\infty}(\Omega\times\Omega)$ with respect to the canonical metric on
the Banach space $L^{\infty}(\Omega\times\Omega)$.

Note also, again since the $L^{\infty}$-norm does not change on a null set,
that we may write

$\displaystyle\|f\|=\max\\{p(f),p(-f)\\},$

where $p(f)$ is the hemi-norm (cf. Gaubert and Vigeral (2012))

$\displaystyle
p(f)=\operatorname*{ess\,sup}_{z\in\Omega\times\Omega}e_{z}(f),$

where $e_{z}$ is the evaluation functional $e_{z}(f)=f(z)$. In the case where
$f$ is continuous on $\Omega\times\Omega\setminus\\{x=y\\}$ we can write

$\displaystyle p(f)=\sup_{z\in\Omega\times\Omega\setminus\\{x=y\\}}e_{z}(f).$

Now we introduce layer maps. More precisely, we consider maps
$T:\Omega\rightarrow\Omega$ which are injective. As explained above these
induce non-expansive maps in the metric $D$:

$\displaystyle(T^{*}d)(x,y):=d(Tx,Ty).$

To any mapping $U:X\to X$ there is a corresponding map
$U\mathcal{F}(d)=\mathcal{F}(Ud)$, and, due to the isometry property of
$\mathcal{F}$, a non-expansive mapping on $X$ becomes a non-expansive mapping
on $Y$. We take as usual a stationary sequence of such layer maps $T_{n}$ and
denote by $f_{n}=\mathcal{F}(T_{n}...T_{1})^{*}(d_{0})$, where $d_{0}$ denotes
the metric coming from the initial norm. Note that $f_{0}=\mathcal{F}d_{0}=0$.

We note that $a(0,n)=p(f_{n}-f_{0})=p(f_{n})$ is a subadditive process, or
subadditive cocycle in the terminology of Gouëzel and Karlsson (2020) (hemi-
metrics work the same since only the triangle inequality and the non-
expansiveness is used). By the subadditive ergodic theorem there is a number
$\lambda$ such that

$\displaystyle\lim_{n\rightarrow\infty}\frac{1}{n}p(f_{n})=\lambda.$

Moreover, (Gouëzel and Karlsson, 2020, Theorem 1.1) asserts that for any
decreasing positive sequence $\delta_{n}\rightarrow 0$ there are times
$n_{i}\rightarrow\infty$ such that for every $i$ and $n<n_{i}$

$\displaystyle p(f_{n_{i}}-f_{n})-p(f_{n_{i}})\leq-n(\lambda-\delta_{n})$

By compactness we may moreover assume that
$h_{n_{i}}(g)=p(f_{n_{i}}-g)-p(f_{n_{i}})$ converges to a metric functional
$h$.

We now follow a similar reasoning to (Gaubert and Vigeral, 2012, p. 349). From
the continuity off the diagonal for elements in $Y$, we see that, given
$f,g\in Y$ and $\delta>0$ there is an off-diagonal point $\hat{z}=(x,y)$
independent of $g$ such that

$\displaystyle p(f-g)-p(f)\geq p(f-g)-e_{\hat{z}}(f)-\delta\geq
e_{\hat{z}}(f-g)-e_{\hat{z}}(f)-\delta=e_{\hat{z}}(-g)-\delta.$ (9)

Let $z_{i}$ be the off-diagonal points from 9 corresponding to $h_{n_{i}}$ in
the inequality above now with $\delta_{n_{i}}$. By again passing to a
subsequence we can ensure that the corresponding sequence of points $z_{i}$
converges to a point $z$ thanks to compactness of $\Omega\times\Omega$ and
$\delta_{i}\rightarrow 0$.

We now end up with two cases, either the limit point $z$ is on the diagonal or
it is off diagonal. Let us begin with the off-diagonal case: If $z=(x,y)$ is
off-diagonal, then 9 gives that

$\displaystyle h(g)\geq-e_{z}(g).$

This implies in view of the multiplicative ergodic theorem in Gouëzel and
Karlsson (2020) that

$\displaystyle\lim_{n\rightarrow\infty}\frac{1}{n}f_{n}(z)=\lambda,$

and more concretely,

$\displaystyle\lim_{n\rightarrow\infty}d_{n}(x,y)^{1/n}=e^{\lambda}.$

Since
$d_{n}(x,y)=(T_{n}...T_{1})^{*}(d_{0})(x,y)=\|T_{n}...T_{1}x-T_{n}...T_{1}y\|$
and $\Omega$ is bounded in $d_{0}$, we can only have this conclusion in case
$\lambda=0$.

In the second case, when $z=(x,x)$ is on the diagonal, then in the notation
$z_{i}=(x_{i},y_{i})$ we get from above on the one hand, for fixed $f_{n}$ and
all $n_{i}>n$

$\displaystyle p(f_{n_{i}}-f_{n})-p(f_{n_{i}})\geq-
e_{z_{i}}(f_{n})-\delta_{n}$

and on the other hand

$\displaystyle p(f_{n_{i}}-f_{n})-p(f_{n_{i}})\leq-n(\lambda-\delta_{n}).$

This implies that

$\displaystyle-n(\lambda-\delta_{n})\geq-\log
d_{n}(x_{i},y_{i})+\log\|x_{i}-y_{i}\|-\delta_{n}$ $\displaystyle
d_{n}(x_{i},y_{i})\geq\|x_{i}-y_{i}\|e^{n(\lambda-\delta_{n}(1+1/n))}$

Choose $\epsilon$, then choose $N$ s.t. $\delta_{N}(1+1/N)<\epsilon$, then for
$n_{i}>n>N$ we get

$\displaystyle
d_{n}(x_{i},y_{i})\geq\|x_{i}-y_{i}\|e^{n(\lambda-\delta_{n}(1+1/n))}\geq\|x_{i}-y_{i}\|e^{n(\lambda-\epsilon)}.$

In words this means that there are sequences of points $x_{i}\rightarrow x$
and $y_{i}\rightarrow x$ which realize the growth rate of the Lipschitz
constant, or put even more strikingly, there is a point $x\in\Omega$ such that
nearby points are separated by the maximum amount ($\lambda$) by the maps
$T_{n}T_{n-1}...T_{1}$.

## Acknowledgments

The first author was supported by the Swedish Research Council grant dnr:
2019-04098. The second author was partly supported by Swiss NSF grant
200020_159581.

## References

  * Agrawal et al. (2020) D. Agrawal, T. Papamarkou, and J. Hinkle.  Wide neural networks with bottlenecks are deep gaussian processes.  _Journal of Machine Learning Research_ , 21(175):1–66, 2020.
  * Aldous and Diaconis (1986) D. Aldous and P. Diaconis.  Shuffling cards and stopping times.  _The American Mathematical Monthly_ , 93(5):333–348, 1986.
  * Arnold (1995) L. Arnold.  Random dynamical systems.  In _Dynamical systems_ , pages 1–43. Springer, 1995.
  * Avelin and Nyström (2021) B. Avelin and K. Nyström.  Neural ODEs as the deep limit of ResNets with constant weights.  _Analysis and Applications_ , 19(03):397–437, 2021.
  * Belkin et al. (2019) M. Belkin, D. Hsu, S. Ma, and S. Mandal.  Reconciling modern machine-learning practice and the classical bias–variance trade-off.  _Proceedings of the National Academy of Sciences_ , 116(32):15849–15854, 2019.
  * Chen et al. (2018) R. T. Chen, Y. Rubanova, J. Bettencourt, and D. Duvenaud.  Neural ordinary differential equations.  In _Proceedings of the 32nd International Conference on Neural Information Processing Systems_ , pages 6572–6583, 2018.
  * Deroin et al. (2007) B. Deroin, V. Kleptsyn, A. Navas, et al.  Sur la dynamique unidimensionnelle en régularité intermédiaire.  _Acta mathematica_ , 199(2):199–262, 2007.
  * Diaconis and Freedman (1999) P. Diaconis and D. Freedman.  Iterated random functions.  _SIAM review_ , 41(1):45–76, 1999.
  * Dunlop et al. (2018) M. M. Dunlop, M. A. Girolami, A. M. Stuart, and A. L. Teckentrup.  How deep are deep gaussian processes?  _Journal of Machine Learning Research_ , 19(54):1–46, 2018.
  * Duvenaud et al. (2014) D. Duvenaud, O. Rippel, R. Adams, and Z. Ghahramani.  Avoiding pathologies in very deep networks.  In _Artificial Intelligence and Statistics_ , pages 202–210. PMLR, 2014.
  * E et al. (2020) W. E, C. Ma, and L. Wu.  Machine learning from a continuous viewpoint, I.  _Science China Mathematics_ , 63(11):2233–2266, 2020.
  * Elfwing et al. (2018) S. Elfwing, E. Uchibe, and K. Doya.  Sigmoid-weighted linear units for neural network function approximation in reinforcement learning.  _Neural Networks_ , 107:3–11, 2018.
  * Gaubert and Vigeral (2012) S. Gaubert and G. Vigeral.  A maximin characterisation of the escape rate of non-expansive mappings in metrically convex spaces.  _Math. Proc. Cambridge Philos. Soc._ , 152(2):341–363, 2012.
  * Glorot and Bengio (2010) X. Glorot and Y. Bengio.  Understanding the difficulty of training deep feedforward neural networks.  In _Proceedings of the thirteenth international conference on artificial intelligence and statistics_ , pages 249–256. JMLR Workshop and Conference Proceedings, 2010.
  * Gouëzel and Karlsson (2020) S. Gouëzel and A. Karlsson.  Subadditive and multiplicative ergodic theorems.  _Journal of the European Mathematical Society_ , 22(6):1893–1915, 2020.
  * Gutiérrez (2019) A. W. Gutiérrez.  The horofunction boundary of finite-dimensional $l_{p}$ spaces.  _Colloquium Mathematicum_ , 155(1):51–65, 2019\.
  * Gutiérrez (2020) A. W. Gutiérrez.  Characterizing the metric compactification of $l_{p}$ spaces by random measures.  _Annals of Functional Analysis_ , 11(2):227–243, 2020.
  * Hanin (2018) B. Hanin.  Which neural net architectures give rise to exploding and vanishing gradients?  In _Proceedings of the 32nd International Conference on Neural Information Processing Systems_ , pages 580–589, 2018.
  * He et al. (2016) K. He, X. Zhang, S. Ren, and J. Sun.  Deep residual learning for image recognition.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 770–778, 2016.
  * Hornik et al. (1989) K. Hornik, M. Stinchcombe, and H. White.  Multilayer feedforward networks are universal approximators.  _Neural networks_ , 2(5):359–366, 1989.
  * Jacot et al. (2018) A. Jacot, F. Gabriel, and C. Hongler.  Neural tangent kernel: convergence and generalization in neural networks.  In _Proceedings of the 32nd International Conference on Neural Information Processing Systems_ , pages 8580–8589, 2018.
  * Karlsson (2019) A. Karlsson.  Elements of a metric spectral theory.  _arXiv preprint arXiv:1904.01398_ , 2019.
  * Karlsson and Ledrappier (2011) A. Karlsson and F. Ledrappier.  Noncommutative ergodic theorems.  In _Geometry, rigidity, and group actions_ , Chicago Lectures in Math., pages 396–418. Univ. Chicago Press, Chicago, IL, 2011.
  * Kingman (1973) J. F. C. Kingman.  Subadditive ergodic theory.  _Annals of Probability_ , 1(6):883–909, 1973\.
  * Lemmens and Nussbaum (2012) B. Lemmens and R. Nussbaum.  _Nonlinear Perron-Frobenius Theory_ , volume 189.  Cambridge University Press, 2012.
  * Li (2018) H. Li.  Analysis on the nonlinear dynamics of deep neural networks: Topological entropy and chaos.  _arXiv preprint arXiv:1804.03987_ , 2018.
  * Liu et al. (2019) X. Liu, S. Si, Q. Cao, S. Kumar, and C.-J. Hsieh.  Neural SDE: Stabilizing neural ode networks with stochastic noise.  _arXiv preprint arXiv:1906.02355_ , 2019.
  * Lu et al. (2017) Z. Lu, H. Pu, F. Wang, Z. Hu, and L. Wang.  The expressive power of neural networks: A view from the width.  In _Advances in Neural Information Processing Systems_ , volume 30. Curran Associates, Inc., 2017.
  * Mei et al. (2018) S. Mei, A. Montanari, and P.-M. Nguyen.  A mean field view of the landscape of two-layer neural networks.  _Proceedings of the National Academy of Sciences_ , 115(33):E7665–E7671, 2018.
  * Nakkiran et al. (2020) P. Nakkiran, G. Kaplun, Y. Bansal, T. Yang, B. Barak, and I. Sutskever.  Deep double descent: Where bigger models and more data hurt.  In _8th International Conference on Learning Representations, $\\{$ICLR$\\}$ 2020,_, 2020.
  * Navas (2018) A. Navas.  On conjugates and the asymptotic distortion of 1-dimensional $c^{1+bv}$ diffeomorphisms.  _arXiv preprint arXiv:1811.06077_ , 2018.
  * Neal (1996) R. M. Neal.  Priors for infinite networks.  In _Bayesian Learning for Neural Networks_ , pages 29–53. Springer, 1996.
  * Neal (2012) R. M. Neal.  _Bayesian learning for neural networks_ , volume 118.  Springer Science & Business Media, 2012.
  * Oseledets (1968) V. I. Oseledets.  A multiplicative ergodic theorem. Characteristic Ljapunov, exponents of dynamical systems.  _Trudy Moskovskogo Matematicheskogo Obshchestva_ , 19:179–210, 1968.
  * Ramachandran et al. (2017) P. Ramachandran, B. Zoph, and Q. V. Le.  Searching for activation functions.  _arXiv preprint arXiv:1710.05941_ , 2017.
  * Robinson and Fallside (1988) A. J. Robinson and F. Fallside.  Static and dynamic error propagation networks with application to speech coding.  In _Neural information processing systems_ , pages 632–641. Citeseer, 1988.
  * Rohwer (1990) R. Rohwer.  The “moving targets” training algorithm.  In _European Association for Signal Processing Workshop_ , pages 100–109. Springer, 1990.
  * Rotskoff and Vanden-Eijnden (2018) G. M. Rotskoff and E. Vanden-Eijnden.  Trainability and accuracy of neural networks: An interacting particle system approach.  _arXiv preprint arXiv:1805.00915_ , 2018.
  * Schaefer et al. (2008) A. M. Schaefer, S. Udluft, and H.-G. Zimmermann.  Learning long-term dependencies with recurrent neural networks.  _Neurocomputing_ , 71(13-15):2481–2488, 2008\.
  * Srivastava et al. (2014) N. Srivastava, G. Hinton, A. Krizhevsky, I. Sutskever, and R. Salakhutdinov.  Dropout: a simple way to prevent neural networks from overfitting.  _The journal of machine learning research_ , 15(1):1929–1958, 2014.
  * Sutskever et al. (2013) I. Sutskever, J. Martens, G. Dahl, and G. Hinton.  On the importance of initialization and momentum in deep learning.  In _International conference on machine learning_ , pages 1139–1147. PMLR, 2013.
  * Thorpe and van Gennip (2018) M. Thorpe and Y. van Gennip.  Deep limits of residual neural networks.  _arXiv preprint arXiv:1810.11741_ , 2018.
  * Tzen and Raginsky (2019) B. Tzen and M. Raginsky.  Neural stochastic differential equations: Deep latent Gaussian models in the diffusion limit.  _arXiv preprint arXiv:1905.09883_ , 2019.
  * Williams and Zipser (1989) R. J. Williams and D. Zipser.  A learning algorithm for continually running fully recurrent neural networks.  _Neural computation_ , 1(2):270–280, 1989.


Further author information: (Send correspondence to J. F.)
J. F.: E-mail<EMAIL_ADDRESS>
# The future looks dark: improving high contrast imaging with hyper-parameter
optimization for data-driven predictive wavefront control

J. Fowler  Department of Astronomy & Astrophysics, University of California,
Santa Cruz, CA, USA  Rebecca Jensen-Clem  Department of Astronomy &
Astrophysics, University of California, Santa Cruz, CA, USA  Maaike A. M. van
Kooten  National Research Council Canada, Herzberg Astronomy and Astrophysics
Research Center, Victoria, Canada  Vincent Chambouleryon  Department of
Astronomy & Astrophysics, University of California, Santa Cruz, CA, USA
Sylvain Cetre  W. M. Keck Observatory, Waimea, Hawaii, USA

###### Abstract

The direct imaging and characterization of exoplanets requires extreme
adaptive optics (XAO), achieving exquisite wavefront correction (upwards of
90$\%$ Strehl) over a narrow field of view (a few arcseconds). For these XAO
systems the temporal error is often a leading term in the error budget,
wherein the wavefront evolves faster than the lag between wavefront sensing
and control. For atmospheres with high-velocity wind layers, this can result
in a wind-driven halo in the coronagraphic dark-zone, limiting sensitivity to
faint, close-in companions. The AO system’s lag-time is often limited by the
wavefront sensor exposure time, especially in the case of fainter guidestars.
Predictive control mitigates the temporal error by predicting the shape of the
wavefront by time the system correction is applied. One such method of
prediction is empirical orthogonal functions (EOF), wherein previous states in
the wavefront sensor history are used to learn linear correlations with a
minimization problem. This method has been demonstrated on-sky at
Subaru/SCExAO and Keck/NIRC2, but has yet to be optimized. With this work as a
starting point, we explore the optimal filter hyper-parameter space for
implementing EOF on-sky, study its stability under varying atmospheric
parameters, and discuss future paths for facilitization of predictive control.
This work not only offers a pathway to optimize Keck and Subaru observing, but
also acts as a pathfinder for predictive control methods with extremely large
telescopes.

###### keywords:

predictive control, temporal error, empirical orthogonal functions (EOF),
hyper-parameter optimization, simulated annealing, Keck, Subaru

##  1 INTRODUCTION

Predictive wavefront control, first suggested by Dessenne in 1997 [1], has
been explored through many implementations over the years, ranging from
classic and predictive Linear Quadratic Gaussian controllers [2, 3], neural
network approaches [4, 5], through reinforcement learning [6, 7, 8, 9], and
through classic minimization problems [10, 11, 12, 13]; for a more complete
review of predictive methods, see (Fowler, 2023)[14]. Predictive controllers
consistently show promising performance with on-sky testing [15, 16, 17], and
in particular, empirical orthogonal functions (EOF) at Keck has shown up to 3x
improvement in contrast at separations of 3-7 $\lambda$/D, and 2x for 3
$\lambda$/D[18].

However, no facilitized implementations of predictive control are running on-
sky and little has been published in the way of optimizing its performance for
regular usage under varying atmospheric conditions. The implementation at
Keck[18] explored the optimal mixing between an integrator and predictor and
optimal gain for the high order correction, but built a predictive filter with
a consistent set of hyper-parameters. Similarly, one work[4] explored the
optimal amount of data for an EOF filter to train on with Subaru/SCExAO
telemetry, but predictive control at SCExAO is run on-sky with a static set of
filter hyper-parameters (private communication, SCExAO team).

In the effort to move predictive control from an on-sky demonstration to a
facilitized mode, we must address how to consistently set up the predictive
controller for varying nightly conditions and make recommendations for a
hyper-parameter optimization strategy during observing nights. Here, we
address the tuning of individual hyper-parameters in the predictive filter as
the first step towards developing such an on-sky strategy.

##  2 TUNEABLE HYPER-PARAMETERS IN THE EOF FILTER

Empirical orthogonal functions (EOF)[10] builds a linear filter F by
assembling l previous states of wavefront sensor information (or history
vectors h) into a training matrix D, and comparing those to known future
states P. The history vector contains n frames of wavefront sensor data (each
with m points for each mode) flattened into a single vector; this is shown
visually in Figure 1. The filter matrix encodes the wind velocities and other
linear system errors.

The logic to build a filter matrix is shown below; the matrix inversion is
done with a least-squares inverse, and uses $\alpha$ as a regularization
parameter.

$\displaystyle min||\mathbf{D}^{T}\mathbf{F}^{T}-\mathbf{P}^{T}||^{2}$ (1)
$\displaystyle\mathbf{F}=((\mathbf{D}^{T})^{\dagger}\mathbf{P}^{T})^{T}$ (2)
$\displaystyle\mathbf{F}=\mathbf{P}\mathbf{D}^{T}(\mathbf{D}\mathbf{D}^{T}+\alpha\mathbf{I})^{-1}$
(3) $\displaystyle\textrm{prediction}=\mathbf{F}\mathbf{h}$ (4)

At every iteration of the real-time-controller, we multiply the filter matrix
F, by the history vector h to predict an iteration one lag-time in the future.
This full logic is visualized in Figure 2.

Figure 1: The history vector contains n frames of wavefront sensor data
flattened into a single vector.

Figure 2: To predict a future iteration with m modes, we multiply a filter
matrix that is m by m$\cdot$n (with m modes and n points in time) by the
current history vector h.

The filter matrix provides a good correction while the atmosphere is driven by
the same linear forces on which it was trained, therefore we must plan to
update the filter matrix every few minutes or have a rolling update of the
filter matrix for on-sky operations. To generate a predictive filter, we
consider the number of frames in our training data l, the number of frames in
our history vector n, and the regularization constant $\alpha$ as tuneable
hyper-parameters of the EOF filter. In simulation, our matrix is well
conditioned, and we can run with an $\alpha=1$, so the simulations presented
in the rest of this work aim to explore optimal values for l and n.

##  3 OPTIMIZING FILTER HYPER-PARAMETERS IN SIMULATION

We simulate the performance of EOF with varying hyper-parameters l and n as
applied to atmospheres with varying windspeeds and Fried parameters. We
simulate a classic single-layer Kolmogorov[19] phase screen with the Taylor
frozen flow approximation[20] with a given Fried parameter r0 and windspeed
for a 10 meter telescope at a resolution of 48x48 pixels. We subtract piston,
tip, and tilt modes, and apply a circular aperture. We sample and correct the
atmosphere at 1kHz, with a lag of two time steps (2 ms), with perfect
wavefront sensing and perfect correction – so the error is indicative only of
error from the control itself.

For each hyper-parameter grid search, we build a predictive EOF filter with
varying numbers of frames of training data (l) and varying history vector
lengths (n), and build a heat map of the the median performance in root mean
square (RMS) error. We compute a residual by subtracting the prediction of the
full phase-screen at each iteration from the known injected phase at each
iteration, and take an RMS of that difference at each point in time, leaving
us with the RMS for each of 10000 iterations. We then take the median of that
time-series distribution for the final value plotted and used for
optimization. Figures 3 and 4 show examples of differing r0 with the same
windspeed and differing windspeed with the same r0 respectively. We note that
the optimal history vector length changes based on different turbulence
conditions.

---
Figure 3:  Heatmap of EOF predictor performance (in median RMS error in nm) as
training data length l and history vector length n vary. Left: r0 of 10 cm.
Middle: r0 of 15 cm. Right: r0 of 20 cm. All of the above simulations have a
windspeed of 11 m/s. Note that the optimal history vector length changes based
on r0.
---
Figure 4:  Heatmap of EOF predictor performance (in median RMS error in nm) as
training data length l and history vector length n vary. Left: windspeed of 7
m/s. Middle: windspeed of 11 m/s. Right: windspeed of 14 m/s. All of the above
simulations have an r0 of 15 cm. Note that the optimal history vector length
is different for the slowest atmosphere. Given the two right-most plots find
an optimal in the corner of our grid, they may also find different optimal
values given an expanded search space.

While these simulations show that the optimal history vector length changes,
simulations consistently show that the optimal amount of training data is the
maximum amount available. We expect that predictor performance due to
different lengths of training data in the filter matrix would be more
impactful on real data, where including more frames of data introduces more
noise into the training (our simulations are lacking measurement noise), and
when the velocity layer will vary over some timescale.

##  4 ON-THE-FLY OPTIMIZATION WITH SIMULATED ANNEALING

If the desired strategy for EOF implementation is to do on-sky hyper-parameter
tuning, building a full grid of filter hyper-parameter runs is likely not
feasible. For these initial grids, we used 10 seconds of data and built a
filter based on up to 24 seconds of training data; even for a system optimized
to run in real time, to do this 600 times would take $\sim$ 6 hours. We
present an alternative way to explore this hyper-parameter space more
optimally with simulated annealing, a probabilistic method for finding the
global optimum of a non-differentiable function.

---
Figure 5:  Simulated annealing samples many values, but it does not accept
every step. This is the run of a simulated annealing algorithm applied to a
grid with a windspeed of 14 m/s and an r0 of 15 cm, chosen because the grid
had the most variation. Left: Entire walk of the simulated annealing to find
the optimal solution. Right: Every point sampled (whether or not it was
accepted) over the 400 steps.

Simulated annealing [21] is a probabilistic way to explore a discrete
parameter space. The algorithm picks a random point of a given step size away
from the current point and accepts it if the point provides a better solution,
or if a Metropolis criterion is met with a given temperature. This process is
named after metal annealing, where metal is heated and cooled strategically
for optimal performance. Pseudo-code of this algorithm and more discussion of
the Metropolis criterion are presented in Appendix A. Figure 5 shows how the
algorithm walks through one of the hyper-parameter search grids shown above.

---
Figure 6:  Simulated annealing is applied to the previously shown dataset,
with a windspeed of 14 m/s and an r0 of 15 cm. The input temperature is 119
and the step size is 5. The dashed line indicates when the algorithm reaches
within $0.1$ nm of the optimal solution. Left: full run. Right: closer detail
on the first few iterations.

We applied simulated annealing to the previously shown filter hyper-parameter
grid for a windspeed of 14 m/s and an r0 of 15 cm, and found that the
algorithm could repeatably find within $0.1$ nm of the optimal value within 10
steps and repeatably recover the exact solution within 400 steps. For an on-
sky strategy, we expect we could use closer to 10 steps, to find a very good
if not exactly optimal solution. However, how quickly the simulated annealing
algorithm converges to near its minimum will be dependent on the performance
of each filter over a set of conditions, so testing with lab-data and on-sky
telemetry will be needed before specific strategies can be recommended.

Figure 6 shows 1000 runs of simulated annealing on the same data set. This
example was run with a step size of 5 (which maps to a single history vector
frame and 50000 training data frames), and a temperature of 119. Initial
testing found that higher temperatures (100-200) worked better for this data
set but no quantitative optimization was done for temperature and step size in
this case, due to how well the algorithm performed with these parameters. A
strength of simulated annealing is its effective application to problems
without significant optimization or previous knowledge of a parameter space.

##  5 PRELIMINARY LAB RESULTS

To further this work we plan to examine $\alpha$ as third hyper-parameter to
optimize, as well as explore how more realistic conditions (e.g., imperfect
noisy data) impact the optimal training data length. To do this, we plan to
run similar grids of optimization on the SEAL (Santa Cruz Extreme AO
Laboratory) testbed [22]. Along with more realistic conditions, SEAL has a
coronagraphic branch with a vector vortex coronagraph [23]; moving to the
bench will allow us to optimize hyper-parameters based on coronagraphic
contrast as well as wavefront error.

Figures 7 and 8 show a preliminary implementation of EOF predictive control on
the SEAL testbed, comparing the performance of EOF (with a filter built with a
5 frame history vector and 60000 frames of training data) to an integrator
with a gain of 0.2. We simulate a single layer atmosphere with a D/r0 of 100
(e.g., a 10 meter telescope with an r0 of 10 cm) and a pupil crossing time of
1s (e.g., 10 m/s wind layer) and sample its phase every millisecond. We apply
turbulence with the same D/r0 in the lab, but with the amplitude scaled down
to the linear regime of the Boston MEMS DM (within 2$\pi$ of phase), and then
applied and corrected with the MEMS DM for 10000 iterations.

We use a high speed Thorlabs Shack-Hartmann Wavefront Sensor with 640 slopes
and a Boston Micromachines MEMS deformable mirror with 24 actuators across the
pupil. Initial results show improvement over a classic integrator both in
reduction of wavefront error (Figure 8), and in the the focal plane PSF images
(Figure 7).

Figure 7: PSFs of predictive control and a classic integrator side-by-side on
the same color-scale. Note: these are preliminary results without robust image
calibration and serve only as an initial visual comparison.

Figure 8: Bench results of an EOF controller as compared to a classic
integrator. RMS wavefront is estimated from the SHWFS residuals as converted
to DM space. Input turbulence is a single wind layer with a pupil crossing
time of 1s and an atmosphere with a D/r0 of 100, with amplitude scaled down
the linear $2\pi$ dynamic range of the Boston MEMS.

Plans for future lab testing include recreating these hyper-parameter grid
searches in an environment more indicative of a true AO system; we will apply
turbulence with a Meadowlark spatial light modulator of 1116 pixels across,
and correct with a woofer/tweeter system of a 97 actuator ALPAO deformable
mirror and the aforementioned Boston MEMS. This will give us the ability to
mimic more realistic effects of measurement error (with turbulence applied as
a much finer spatial scale than we can correct.)

##  6 CONCLUSIONS

In conclusion, we explore future facilitization of predictive wavefront
control by looking at optimal hyper-parameters for the training data and
history vector length of the predictive filter. We find that in simulation,
optimal filter hyper-parameters vary based on the atmospheric turbulence, and
expect that we will need to optimize on-sky. We explore simulated annealing as
a way to probe filter hyper-parameters more efficiently, and find it can
repeatably converge to within 0.1 nm RMS of the optimal solution within 10
iterations.

We present preliminary lab results, that compare EOF to a classic integrator,
and see promising performance from EOF. This lab implementation is a precursor
to future exploration of $\alpha$ in the EOF matrix inversion, as well as an
opportunity to study the impact of these corrections on the coronagraphic dark
hole. We also hope to test this optimization on-sky with Keck/NIRC2 or
Subaru/SCExAO.

Facilitizing predictive control for regular on-sky use will not only improve
current high contrast imaging for large telescopes, but will show technology
maturation for its use in upcoming extremely large telescopes. When predictive
control is in use, the future of our coronagraphic dark hole is looking dark.

##  Appendix A SIMULATED ANNEALING ALGORITHM

Here we further describe the Metropolis criterion, and present Python pseudo-
code of the simulated annealing algorithm to demonstrate the application of
this algorithm to our data. For each iteration, the algorithm picks a new
point, evaluates if that point is a better solution, and if not decides to
keep the new point only if it meets the Metropolis criterion, which is met by
comparing a randomly generated number to a function built from a cooling input
temperature and how far the evaluation of the new point is from the previous.
We note that simulated annealing and the Metropolis criterion make this
algorithm vulnerable to getting trapped in a deep local minimum, and plan to
study this impact with more realistic lab data. For a given iteration i:

$\displaystyle\textrm{cooled temperature}=\frac{\textrm{temperature}}{i+1}$
(5)
$\displaystyle\textrm{difference}=\textrm{evaluation}(i-1)-\textrm{evaluation}(i)$
(6)
$\displaystyle\textrm{metropolis}=\exp\left(\frac{\textrm{difference}}{\textrm{cooled
temperature}}\right)$ (7)

This logic runs for a given number of steps, but on-sky this could also be a
cutoff based on how much the algorithm is still improving the solution or some
minimum allowable wavefront error.

⬇

def simulated_annealing(objective, bounds, n_iterations, step_size, temp):

# start with an initial guess

best = preset or random guess

# evaluate it and set it as a starting point

best_eval = objective(best)

current, current_eval = best, best_eval

# for every iteration

for i in range(n_iterations):

# take a step

candidate = move from current

# evaluate the new candidate

candidate_eval = objective(candidate)

# check for new best solution

if candidate_eval < best_eval:

# store new best point

best, best_eval = candidate, candidate_eval

# otherwise, use metropolis to see if we keep the new point

# or stay with the previous guess

# this determines where the guess next iteration starts from

diff = candidate_eval - current_eval

t = temp/(i + 1)

# calculate metropolis acceptance criterion

metropolis = np.exp(-diff / t)

# check if we should keep the new point

if diff < 0 or random number < metropolis:

# store the new current point

current, current_eval = candidate, candidate_eval

# after however many iterations return the best solution

return best, best_eval

Figure 9: Pseduo-code of simulated annealing algorithm, heavily influenced by
Machine Learning Mastery. The details of making consistent integer steps over
a two dimensional grid make this algorithm more complex, so we have chosen to
present representative code as opposed to the literal implementation.

###### Acknowledgements.

This work was supported by NSF ATI Grant 2008822. Thank you to the folks at
Keck Observatory and members of the SCExAO team for discussing their
implementations of EOF.

## References

  * [1] Dessenne, C., Madec, P. Y., and Rousset, G., “Modal prediction forclosed-loop adaptive optics,” Optics Letters 22, 1535–1537 (Oct. 1997).
  * [2] Sivo, G., Kulcsár, C., Conan, J.-M., Raynaud, H.-F., Gendron, É., Basden, A., Vidal, F., Morris, T., Meimon, S., Petit, C., Gratadour, D., Martin, O., Hubert, Z., Sevin, A., Perret, D., Chemla, F., Rousset, G., Dipper, N., Talbot, G., Younger, E., Myers, R., Henry, D., Todd, S., Atkinson, D., Dickson, C., and Longmore, A., “First on-sky SCAO validation of full LQG control with vibration mitigation on the CANARY pathfinder,” Optics Express 22, 23565 (Sept. 2014).
  * [3] Poyneer, L. A., Ammons, S. M., Kim, M. K., Bauman, B., Terrel-Perez, J., Lemmer, A. J., and Nguyen, J., “Laboratory demonstration of the prediction of wind-blown turbulence by adaptive optics at 8 kHz with use of LQG control,” Applied Optics 62, 1871 (Mar. 2023).
  * [4] Wong, A. P., Norris, B. R. M., Tuthill, P. G., Scalzo, R., Lozi, J., Vievard, S., and Guyon, O., “Predictive control for adaptive optics using neural networks,” Journal of Astronomical Telescopes, Instruments, and Systems 7, 019001 (Jan. 2021).
  * [5] Swanson, R., Lamb, M., Correia, C. M., Sivanandam, S., and Kutulakos, K., “Closed loop predictive control of adaptive optics systems with convolutional neural networks,” Monthly Notices of the RAS 503, 2944–2954 (May 2021).
  * [6] Haffert, S. Y., Males, J. R., Close, L. M., Van Gorkom, K., Long, J. D., Hedglen, A. D., Guyon, O., Schatz, L., Kautz, M., Lumbres, J., Rodack, A., Knight, J. M., Sun, H., and Fogarty, K., “Data-driven subspace predictive control of adaptive optics for high-contrast imaging,” arXiv e-prints , arXiv:2103.07566 (Mar. 2021).
  * [7] Landman, R., Haffert, S. Y., Radhakrishnan, V. M., and Keller, C. U., “Self-optimizing adaptive optics control with reinforcement learning for high-contrast imaging,” Journal of Astronomical Telescopes, Instruments, and Systems 7, 039002 (July 2021).
  * [8] Nousiainen, J., Engler, B., Kasper, M., Rajani, C., Helin, T., Heritier, C. T., Quanz, S. P., and Glauser, A. M., “Laboratory experiments of model-based reinforcement learning for adaptive optics control,” Journal of Astronomical Telescopes, Instruments, and Systems 10, 019001 (Jan. 2024).
  * [9] Pou, B., Smith, J., Quinones, E., Martin, M., and Gratadour, D., “Model-free reinforcement learning with a non-linear reconstructor for closed-loop adaptive optics control with a pyramid wavefront sensor,” in [Adaptive Optics Systems VIII ], Schreiber, L., Schmidt, D., and Vernet, E., eds., Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 12185, 121852U (Aug. 2022).
  * [10] Guyon, O. and Males, J., “Adaptive Optics Predictive Control with Empirical Orthogonal Functions (EOFs),” arXiv e-prints , arXiv:1707.00570 (July 2017).
  * [11] Jensen-Clem, R., Bond, C. Z., Cetre, S., McEwen, E., Wizinowich, P., Ragland, S., Mawet, D., and Graham, J., “Demonstrating predictive wavefront control with the Keck II near-infrared pyramid wavefront sensor,” in [Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series ], Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 11117, 111170W (Sept. 2019).
  * [12] van Kooten, M. A. M., Doelman, N., and Kenworthy, M., “Robustness of prediction for extreme adaptive optics systems under various observing conditions. An analysis using VLT/SPHERE adaptive optics data,” Astronomy and Astrophysics 636, A81 (Apr. 2020).
  * [13] Fowler, J., Van Kooten, M. A. M., and Jensen-Clem, R., “Battle of the predictive wavefront controls: comparing data and model-driven predictive control for high contrast imaging,” in [Adaptive Optics Systems VIII ], Schreiber, L., Schmidt, D., and Vernet, E., eds., Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 12185, 1218582 (Aug. 2022).
  * [14] Fowler, J. and Landman, R., “Tempestas ex machina: A review of machine learning methods for wavefront control,” arXiv e-prints , arXiv:2309.00730 (Sept. 2023).
  * [15] Guyon, O., Sevin, A., Ferreira, F., Ltaief, H., Males, J., Deo, V., Gratadour, D., Cetre, S., Martinache, F., Lozi, J., Vievard, S., Fruitwala, N., Bos, S., and Skaf, N., “Adaptive optics real-time control with the compute and control for adaptive optics (Cacao) software framework,” in [Adaptive Optics Systems VII ], Schreiber, L., Schmidt, D., and Vernet, E., eds., Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 11448, 114482N (Dec. 2020).
  * [16] Dessenne, C., Madec, P.-Y., and Rousset, G., “Optimization of a Predictive Controller for Closed-Loop Adaptive Optics,” Applied Optics 37, 4623–4633 (July 1998).
  * [17] Marquis, L., Raynaud, H., Galland, N., De La Rosa, J., Garcia, I., Tubio, O., Garcia-Talavera, M., Lombardi, G., Lopez, M., Reverte, D., and Kulcsar, C., “First on-sky tests of lqg control for a 10m-class telescope: prelude on the gran telescopio canarias adaptive optics system,” Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series in same proceedings (2024).
  * [18] van Kooten, M. A. M., Jensen-Clem, R., Cetre, S., Ragland, S., Bond, C. Z., Fowler, J., and Wizinowich, P., “Predictive wavefront control on Keck II adaptive optics bench: on-sky coronagraphic results,” Journal of Astronomical Telescopes, Instruments, and Systems 8, 029006 (Apr. 2022).
  * [19] Kolmogorov, A., “The Local Structure of Turbulence in Incompressible Viscous Fluid for Very Large Reynolds’ Numbers,” Akademiia Nauk SSSR Doklady 30, 301–305 (Jan. 1941).
  * [20] Taylor, G. I., “The Spectrum of Turbulence,” Proceedings of the Royal Society of London Series A 164, 476–490 (Feb. 1938).
  * [21] Laarhoven, P. J. M., [Simulated annealing: Theory and applications ] (1987).
  * [22] Jensen-Clem, R., Dillon, D., Gerard, B., van Kooten, M. A. M., Fowler, J., Kupke, R., Cetre, S., Sanchez, D., Hinz, P., Laguna, C., Doelman, D., and Snik, F., “The Santa Cruz Extreme AO Lab (SEAL): design and first light,” in [Techniques and Instrumentation for Detection of Exoplanets X ], Shaklan, S. B. and Ruane, G. J., eds., Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 11823, 118231D (Sept. 2021).
  * [23] Moreno, A., Chambouleyron, R., Jensen-Clem, R., Dillon, D., Hinz, P., and Macintosh, B., “Implementation and characterization of the vector vortex coronagraph on the seal testbed,” Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series in same proceedings (2024).