# The Mori-Zwanzig formulation of deep learning

Daniele Venturi<EMAIL_ADDRESS>Xiantao Li<EMAIL_ADDRESS>Department of
Applied Mathematics, UC Santa Cruz, Santa Cruz, CA 95064  Department of
Mathematics, Pennsylvania State University, State College, PA 16801

###### Abstract

We develop a new formulation of deep learning based on the Mori-Zwanzig (MZ)
formalism of irreversible statistical mechanics. The new formulation is built
upon the well-known duality between deep neural networks and discrete
stochastic dynamical systems, and it allows us to directly propagate
quantities of interest (conditional expectations and probability density
functions) forward and backward through the network by means of exact linear
operator equations. Such new equations can be used as a starting point to
develop new effective parameterizations of deep neural networks, and provide a
new framework to study deep-learning via operator theoretic methods. The
proposed MZ formulation of deep learning naturally introduces a new concept,
i.e., the memory of the neural network, which plays a fundamental role in low-
dimensional modeling and parameterization. By using the theory of contraction
mappings, we develop sufficient conditions for the memory of the neural
network to decay with the number of layers. This allows us to rigorously
transform deep networks into shallow ones, e.g., by reducing the number of
neurons per layer (using projection operators), or by reducing the total
number of layers (using the decay property of the memory operator).

††journal: Research in the Mathematical Sciences

##  1 Introduction

It has been recently shown that new insights on deep learning can be obtained
by regarding the process of training a deep neural network as a discretization
of an optimal control problem involving nonlinear differential equations [11,
10, 15]. One attractive feature of this formulation is that it allows us to
use tools from dynamical system theory such as the Pontryagin maximum
principle or the Hamilton-Jacobi-Bellman equation to study deep learning from
a rigorous mathematical perspective [26, 16, 31]. For instance, it has been
recently shown that by idealizing deep residual networks as continuous-time
dynamical systems it is possible to derive sufficient conditions for universal
approximation in $L^{p}$, which can also be understood as approximation
theories using flow maps generated by dynamical systems [27].

In the spirit of modeling a deep neural network as a flow of a discrete
stochastic dynamical system, in this paper we develop a new formulation of
deep learning based on the Mori-Zwanzig (MZ) formalism. The MZ formalism was
originally developed in statistical mechanics [35, 57] to formally integrate
under-resolved phase variables in nonlinear dynamical systems by means of a
projection operator. One of the main features of such formulation is that it
allows us to systematically derive exact evolution equations for quantities of
interest, e.g., macroscopic observables, based on microscopic equations of
motion [8, 17, 19, 6, 4, 48, 9, 54, 55].

In the context of deep learning, the MZ formalism can be used to reduce the
total number of degrees of freedom of the neural network, e.g., by reducing
the number of neurons per layer (using projection operators), or by
transforming deep networks into shallows networks, e.g., using some
approximation of the MZ memory operator111As is well known, the approximation
properties of neural networks usually improve as the size of the network, and
therefore the number of weights and biases (i.e., the degrees of freedom)
increases.. Computing the solution of the MZ equation for deep learning is not
an easy task. One of the main challenges is the approximation of the memory
term and the fluctuation (noise) term, which encode the interaction between
the so-called orthogonal dynamics and the dynamics of the quantity of
interest. In the context of neural networks, the orthogonal dynamics is
essentially a discrete high-dimensional flow governed by a difference equation
that is hard to solve. Despite these difficulties, the MZ equation of deep
learning is formally exact, and can be used as a starting point to build
useful approximations and parameterizations that target the output function
directly. Moreover, it provides a new framework to study deep-learning via
operator theoretic approaches. For example, the analysis of the memory term in
the MZ formulation may shed light on the behaviour of recent neural network
architectures such as the long short-term memory (LSTM) network [43, 14].

This paper is organized as follows. In section 2 we briefly review the
formulation of deep learning as a control problem involving a discrete
stochastic dynamical system. In section 3 we introduce the composition and
transfer operators associated with the neural network. Such operators are the
discrete analogues of the stochastic Koopman and Frobenious-Perron operators
in classical (continuous time) nonlinear dynamics. In the neural network
setting the composition and transfer operators are kernel integral operators,
the kernel being the transition probability density (conditional probability)
between one layer and the next. In section 4.3 we present the MZ formulation
of deep learning and derive the operator equations at the basis of our theory.
In section 5 we study a particular class of projection operators, i.e., Mori’s
projections [55], and derive sufficient conditions under which the MZ memory
term of the neural network decays with the number of layers. This allows us to
approximate the MZ memory term with with just a few term and re-parameterize
the network accordingly. The main findings are summarized in section 6. We
also include two appendices in which we establish theoretical results
concerning the composition and transfer operators for various types of inter-
layer random noise distributions, and different training paradigms. In
particular, we discuss the “train-over-weights” versus the “train-over-noise”
paradigms and notice that the “train-over-noise” paradigm opens the
possibility to approximately encode, send, and decode secret messages in a
fully trained neural network by choosing the probability distribution of the
random noise appropriately.

##  2 Modeling neural networks as discrete stochastic dynamical systems

We model the input-output map of a neural network as a discrete stochastic
dynamical system of the form (see Figure 1)

$\bm{X}_{n+1}=\bm{F}_{n}(\bm{X}_{n},\bm{w}_{n})+\bm{M}_{n}(\bm{X}_{n},\bm{\beta}_{n})\bm{\xi}_{n},\quad{n=0,1,\cdots,L-1}.$
(1)

Here the index $n$ labels a specific layer in the network,
$\bm{X}_{n}\in\mathbb{R}^{d_{n}}$ represents the input of the $(n+1)$-th layer
(vector of dimension $d_{n}$, the dimension of each layer can vary with $n$),
$\bm{X}_{n+1}\in\mathbb{R}^{d_{n+1}}$ represents the output of the $(n+1)$-th
layer, $\bm{M}_{n}(\bm{X}_{n},\bm{\beta}_{n})\in\mathbb{R}^{d_{n}\times
m_{n}}$ is a state-dependent matrix and
$\\{\bm{\xi}_{0},\ldots,\bm{\xi}_{L-1}\\}$ is a discrete (vector-valued)
Markov process. We allow the initial state $\bm{X}_{0}$ to be random as well,
which can be directly connected to a data set in a training algorithm. The
mappings
$\bm{F}_{n}:\mathbb{R}^{d_{n}}\times\mathbb{R}^{q_{n}}\mapsto\mathbb{R}^{d_{n+1}}$
and
$\bm{M}_{n}:\mathbb{R}^{d_{n}}\times\mathbb{R}^{r_{n}}\mapsto\mathbb{R}^{d_{n+1}\times{m_{n}}}$
depend on the parameters $\bm{w}_{n}\in\mathbb{R}^{q_{n}}$ and
$\bm{\beta}_{n}\in\mathbb{R}^{r_{n}}$, respectively. Note that the dimension
of the random vectors $\bm{X}_{n}$ and $\bm{X}_{n+1}$ can vary from layer to
layer, e.g., in encoding or decoding neural networks. In a supervised learning
setting, the degrees of freedom

$\bm{w}=\\{\bm{w}_{0},\ldots,\bm{w}_{L-1}\\},\qquad\bm{\beta}=\\{\bm{\beta}_{0},\ldots,\bm{\beta}_{L-1}\\},$
(2)

are usually determined by optimizing a suitable performance metric depending
on the network output.

Figure 1: Sketch of the stochastic neural network model (1). The network is
driven by a discrete Markov random process
$\\{\bm{x}_{0},\ldots,\bm{\xi}_{L-1}\\}$. Note that the dimension of the
vectors $\bm{X}_{n}$ and $\bm{X}_{n+1}$ can vary from layer to layer, e.g., in
encoding or decoding neural networks.

A neural network is said to be recurrent if
$\bm{F}_{n}(\bm{X},\bm{w}_{n})=\bm{F}(\bm{X},\bm{w}_{n})$ for all
$n=0,\ldots,L-1$. In this setting, the domain and image of $\bm{F}$ have the
same dimension, and the discrete dynamical system (1) can be written as222A
more general form of (4) is
$\bm{X}_{n+1}=\bm{F}(\bm{X}_{n},\bm{w}_{n},\bm{\xi}_{n}).$ (3)

$\bm{X}_{n+1}=\bm{F}(\bm{X}_{n},\bm{w}_{n})+\bm{M}(\bm{X}_{n},\bm{\beta}_{n})\bm{\xi}_{n}.$
(4)

Furthermore, if $\bm{w}_{n}=\bm{w}$ (each layer of the network shares the same
weights) and $\bm{M}(\bm{X}_{n},\bm{\beta}_{n})=\bm{0}$ then the network is
termed time-invariant [28] (or autonomous). The mapping $\bm{F}$ is usually
written as

$\bm{F}(\bm{X}_{n},\bm{w}_{n})=\bm{\varphi}(\bm{W}_{n}\bm{X}_{n}+\bm{b}_{n})\qquad\bm{w}_{n}=\\{\bm{W}_{n},\bm{b}_{n}\\},$
(5)

where $\bm{\varphi}:\mathbb{R}^{N}\mapsto\mathbb{R}^{N}$ is the activation
function of the network, $\bm{W}_{n}$ is a $N\times N$ weight matrix and
$\bm{b}_{n}$ is a bias vector. In this paper we will study mostly the case
where $\bm{M}(\bm{X}_{n},\bm{\beta}_{n})$ in (4) is the identity matrix. This
yields a neural network model of the form (see Figure 2)

$\bm{X}_{n+1}=\bm{F}(\bm{X}_{n},\bm{w}_{n})+\bm{\xi}_{n}.$ (6)

We emphasize that (6) and (1) are different from recently proposed stochastic
neural network models in which weights and biases are assumed to be random
[13, 51]. By adding random noise to the output of each neural network layer we
are essentially adding an infinite number of degrees of freedom to our system.
This gives us flexibility and allows us to rethink the process of training the
neural network from a probabilistic perspective. For instance, random noise
allows us to approximately encode secret messages in fully trained neural
networks by selecting an appropriate transition probability for the noise
process (see B). The neural network output

$f(\bm{x})=\bm{\alpha}\cdot\mathbb{E}\left[\bm{X}_{L}|\bm{X}_{0}=\bm{x}\right],$
(7)

($\bm{\alpha}$ output weights,
$\mathbb{E}\left[\bm{X}_{L}|\bm{X}_{0}=\bm{x}\right]$ conditional expectation
of $\bm{X}_{L}$ given $\bm{X}_{0}=\bm{x}$) can represent a deterministic
multivariate function, e.g., an approximation of a nonlinear functional
defined on a compact subset of a Hilbert or a Banach space [49, 47].

Figure 2: Sketch of the stochastic neural network model (6). Note that the
transfer function $\bm{F}$ is the same in each layer. This implies that the
random vectors $\\{\bm{X}_{1},\ldots,\bm{X}_{L}\\}$ all have the same
dimension.

##  3 Composition and transfer operators for neural networks

Given a continuous-time Markov process $\bm{X}(t)$ and a measurable phase
space function $\bm{u}(\bm{x})$, the conditional expectation of
$\bm{u}(\bm{X}(t))$ given $\bm{X}(s)=\bm{y}$ is defined as (see, e.g., [38])

$\mathbb{E}\\{\bm{u}(\bm{X}(t))|\bm{X}(s)=\bm{y}\\}=\int\bm{u}(\bm{x})p_{t|s}(\bm{x}|\bm{y})d\bm{x},$
(8)

where $p_{t|s}(\bm{x}|\bm{y})$ is the conditional transition density333Recall
that if $\bm{X}(t)$ is the solution to a continuous-time SDEs driven by a
multiplicative Wiener process, then the conditional transition density
satisfies the well-known backward Kolmogorov equation [42], and the dynamics
of the conditional expectation (8) is generated by a hypoelliptic
deterministic operator (see [56, 44] and [46, Proposition 6]). of $\bm{X}(t)$
given $\bm{X}(s)=\bm{y}$. It follows from the theory of Markov processes that
for all $t\geq s$

$\mathbb{E}\\{\bm{u}(\bm{X}(t))|\bm{X}(s)=\bm{y}\\}=\int\mathbb{E}\\{\bm{u}(\bm{X}(t))|\bm{X}(s)=\bm{x}\\}p_{t|s}(\bm{x}|\bm{y})d\bm{x}.$
(9)

The transition density $p_{t|s}(\bm{x}|\bm{y})$ allows us to compute also the
PDF $p(\bm{x},t)$ of the random vector $\bm{X}(t)$ at time $t$ given the PDF
$p(\bm{x},s)$ of $\bm{X}(s)$ at any time $s\leq t$. In particular, we have

$p(\bm{x},t)=\int p_{t|s}(\bm{x}|\bm{y})p(\bm{y},s)d\bm{y}.$ (10)

The integral operators defined in equations (10) and (9) are known,
respectively, as Frobenious-Perron (or transfer) operator and stochastic
Koopman [46, 56] (or composition) operator associated with the Markov process
$\bm{X}(t)$. Both operators have the same conditional transition density
$p_{t|s}(\bm{x}|\bm{y})$ as kernel, and they are adjoint to one another in
$L^{2}$. However, note that in (9) the conditional transition density is
integrated with respect to the variable $\bm{x}$, while in (10) integration is
with respect to $\bm{y}$. Equations (9) and (10) embody two parallel
descriptions of statistical properties associated with Markov chains, and can
be applied repeatedly to progressively keep track of the dynamics of the
conditional expectation and the PDF of the continuous-time Markov process
$\bm{X}(t)$

Let us now derive the composition and transfer operators associated with the
neural network model (1), which map, respectively, the conditional expectation
$\mathbb{E}\left\\{\bm{u}(\bm{X}_{L})|\bm{X}_{n}=\bm{x}\right\\}$ and
$p_{n}(\bm{x})$ (the probability density of $\bm{X}_{n}$) forward and backward
across the network. To this end, we assume that
$\\{\bm{\xi}_{0},\ldots,\bm{\xi}_{L-1}\\}$ in (4) is a (vector-valued)
discrete Markov process. This implies that
$\\{\bm{X}_{0},\bm{X}_{1},\bm{X}_{2},\ldots,\bm{X}_{L}\\}$ is also a discrete
Markov process. In this assumption, the joint probability density function
(PDF) of the random variables $\\{\bm{X}_{0},\ldots,\bm{X}_{L}\\}$, i.e.,
joint PDF of the state of the entire neural network, can be factored444In
equation (11) we used the shorthand notation $p_{i|j}(\bm{x}|\bm{y})$ to
denote the conditional probability density function of the random vector
$\bm{X}_{i}$ given $\bm{X}_{j}=\bm{x}_{j}$. With this notation we have that
the conditional probability density of $\bm{X}_{i}$ given $\bm{X}_{i}=\bm{y}$
is $p_{i|i}(\bm{x}|\bm{y})=\delta(\bm{x}-\bm{y})$, where $\delta(\cdot)$ is
the Dirac delta function. as

$p(\bm{x}_{0},\ldots,\bm{x}_{L})=p_{L|L-1}(\bm{x}_{L}|\bm{x}_{L-1})p_{L-1|L-2}(\bm{x}_{L-1}|\bm{x}_{L-2})\cdots
p_{1|0}(\bm{x}_{1}|\bm{x}_{0})p_{0}(\bm{x}_{0}).$ (11)

By using the identity

$p(\bm{x}_{k+1},\bm{x}_{k})=p_{k+1|k}(\bm{x}_{k+1}|\bm{x}_{k})p_{k}(\bm{x}_{k})=p_{k|k+1}(\bm{x}_{k}|\bm{x}_{k+1})p_{k+1}(\bm{x}_{k+1})$
(12)

we see that the chain of transition probabilities (11) can be reverted,
yielding

$p(\bm{x}_{0},\ldots,\bm{x}_{L})=p_{0|1}(\bm{x}_{0}|\bm{x}_{1})p_{1|2}(\bm{x}_{1}|\bm{x}_{2})\cdots
p_{L-1|L}(\bm{x}_{L-1}|\bm{x}_{L})p_{L}(\bm{x}_{L}).$ (13)

From these expression, it follows that

$p_{n|q}(\bm{x}|\bm{y})=\int
p_{n|j}(\bm{x}|\bm{z})p_{j|q}(\bm{z}|\bm{y})d\bm{z},$ (14)

for all indices $n$, $j$ and $q$ in $\\{0,\ldots,L\\}$, excluding $n=j=q$ (see
footnote 4). The transition probability equation (14) is known as discrete
Chapman-Kolmogorov equation and it allows us to define the transfer operator
mapping the PDF $p_{n}(\bm{x}_{n})$ into $p_{n+1}(\bm{x}_{n+1})$, together
with the composition operator for the conditional expectation
$\mathbb{E}\\{\bm{u}(\bm{x}_{L})|\bm{X}_{n}=\bm{x}_{n}\\}$. As we shall will
see hereafter, the discrete composition and transfer operators are adjoint to
one another.

###  3.1 Transfer operator

Let us denote by $p_{q}(\bm{x})$ the PDF of $\bm{X}_{q}$, i.e., the output of
the $q$-th neural network layer. We first define the operator that maps
$p_{q}(\bm{x})$ into $p_{n}(\bm{x})$. By integrating the joint probability
density of $\bm{X}_{n}$ and $\bm{X}_{q}$, i.e.,
$p_{n|q}(\bm{x}|\bm{y})p_{q}(\bm{y})$ with respect to $\bm{y}$ we immediately
obtain

$p_{n}(\bm{x})=\int p_{n|q}(\bm{x}|\bm{y})p_{q}(\bm{y})d\bm{y}.$ (15)

At this point, it is convenient to define the operator

$\mathcal{N}(n,q)f(\bm{x})=\int p_{n|q}(\bm{x}|\bm{y})f(\bm{y})d\bm{y}.$ (16)

$\mathcal{N}(n,q)$ is known as transfer operator [9]. From a mathematical
viewpoint $\mathcal{N}(n,q)$ is an integral operator with kernel
$p_{n|q}(\bm{x},\bm{y})$, i.e., the transition density integrated “from the
right”. It follows from the Chapman-Kolmogorov identity (14) that the set of
integral operators $\\{\mathcal{N}(n,q)\\}$ forms a group. Namely,

$\mathcal{N}(n,q)=\mathcal{N}(n,j)\mathcal{N}(j,q),\qquad\mathcal{N}(j,j)=\mathcal{I},\qquad\forall
n,j,q\in\\{0,\ldots,L\\}.$ (17)

The operator $\mathcal{N}$ allows us to map the one-layer PDF, e.g., the PDF
of $\bm{X}_{q}$, either forward or backward across the neural network (see
Figure 3). As an example, consider a network with four layers with states
$\bm{X}_{0}$ (input), $\bm{X}_{1}$, $\bm{X}_{2}$, $\bm{X}_{3}$, and
$\bm{X}_{4}$ (output). Then Eq. (16) implies that,

$\displaystyle
p_{2}(\bm{x})=\underbrace{\mathcal{N}(2,1)\mathcal{N}(1,0)}_{\mathcal{N}(2,0)}p_{0}(\bm{x})=\underbrace{\mathcal{N}(2,3)\mathcal{N}(3,4)}_{\mathcal{N}(2,4)}p_{4}(\bm{x}).$

In summary, we have

$p_{n}(\bm{x})=\mathcal{N}(n,q)p_{q}(\bm{x})\qquad\forall
n,q\in\\{0,\ldots,L\\},$ (18)

where

$\mathcal{N}(n,q)p_{q}(\bm{x})=\int
p_{n|q}(\bm{x}|\bm{y})p_{q}(\bm{y})d\bm{y}.$ (19)

We emphasize that modeling PDF dynamics via neural networks has been studied
extensively in machine learning, e.g., in the theory of normalizing flows for
density estimation or variational inference [41, 23, 45].

###  3.2 Composition operator

For any measurable deterministic function $\bm{u}(\bm{x})$, the conditional
expectation of $\bm{u}(\bm{X}_{j})$ given $\bm{X}_{n}=\bm{x}$ is defined as

$\mathbb{E}\left\\{\bm{u}(\bm{X}_{j})|\bm{X}_{n}=\bm{x}\right\\}=\int\bm{u}(\bm{y})p_{j|n}(\bm{y}|\bm{x})d\bm{y}.$
(20)

A substitution of (14) into (20) yields

$\mathbb{E}\left\\{\bm{u}(\bm{X}_{j})|\bm{X}_{n}=\bm{x}\right\\}=\int\mathbb{E}\left\\{\bm{u}(\bm{X}_{j})|\bm{X}_{q}=\bm{y}\right\\}p_{q|n}(\bm{y}|\bm{x})d\bm{y},$
(21)

which holds for all $j,n,q\in\\{0,\ldots,L-1\\}$. At this point we define the
integral operator

$\mathcal{M}(n,q)f(\bm{x})=\int f(\bm{y})p_{q|n}(\bm{y}|\bm{x})d\bm{y},$ (22)

which is known as composition [9] or “stochastic Koopman” [46, 56] operator.
Thanks to the Chapman-Kolmogorov identity (14), the set of operators
$\\{\mathcal{M}(q,j)\\}$ forms a group, i.e.,

$\mathcal{M}(n,q)=\mathcal{M}(n,j)\mathcal{M}(j,q),\qquad\mathcal{M}(j,j)=\mathcal{I},\qquad\forall
n,j,q\in\\{0,\ldots,L\\}.$ (23)

Equation (23) allows us to map the conditional expectation (20) of any
measurable phase space function $\bm{u}(\bm{X}_{j})$ forward or backward
through the network. As an example, consider again a neural network with four
layers and states $\\{\bm{X}_{0},\dots,\bm{X}_{4}\\}$. We have

$\displaystyle\mathbb{E}\\{\bm{u}(\bm{X}_{j})|\bm{X}_{2}=\bm{x}\\}=$
$\displaystyle\mathcal{M}(2,3)\mathcal{M}(3,4)\mathbb{E}\\{\bm{u}(\bm{X}_{j})|\bm{X}_{4}=\bm{x}\\}$
$\displaystyle=$
$\displaystyle\mathcal{M}(2,1)\mathcal{M}(1,0)\mathbb{E}\\{\bm{u}(\bm{X}_{j})|\bm{X}_{0}=\bm{x}\\}.$
(24)

Equation (24) holds for every $j\in\\{0,..,4\\}$. Of particular interest in
machine-learning context is the conditional expectation of
$\bm{u}(\bm{X}_{L})$ (network output) given $\bm{X}_{0}=\bm{x}$, which can be
computed as

$\displaystyle\mathbb{E}\\{\bm{u}(\bm{X}_{L})|\bm{X}_{0}=\bm{x}\\}$
$\displaystyle=\mathcal{M}(0,L)\bm{u}(\bm{x}),$
$\displaystyle=\mathcal{M}(0,1)\mathcal{M}(1,2)\cdots\mathcal{M}(L-1,L)\bm{u}(\bm{x}),$
(25)

i.e., by propagating
$\bm{u}(\bm{x})=\mathbb{E}\\{\bm{u}(\bm{X}_{L})|\bm{X}_{L}=\bm{x}\\}$ backward
through the neural network using single layer operators $\mathcal{M}(i-1,i)$.
Similarly, we can compute, e.g.,
$\mathbb{E}\\{\bm{u}(\bm{X}_{0})|\bm{X}_{L}=\bm{x}\\}$ as

$\mathbb{E}\\{\bm{u}(\bm{X}_{0})|\bm{X}_{L}=\bm{x}\\}=\mathcal{M}(L,0)\bm{u}(\bm{x}).$
(26)

For subsequent analysis, it is convenient to define

$\bm{q}_{n}(\bm{x})=\mathbb{E}\\{\bm{u}(\bm{X}_{L})|\bm{X}_{L-n}=\bm{x}\\}.$
(27)

In this way, if $\mathbb{E}\\{\bm{u}(\bm{X}_{L})|\bm{X}_{n}=\bm{x}\\}$ is
propagated backward through the network by $\mathcal{M}(n-1,n)$, then
$\bm{q}_{n}(x)$ is propagated forward by the operator

$\displaystyle\mathcal{G}(n,q)=$ $\displaystyle\mathcal{M}(L-n,L-q).$ (28)

In fact, equations (27)-(28) allow us to write (25) in the equivalent form

$\displaystyle\bm{q}_{L}(\bm{x})=$
$\displaystyle\mathcal{G}(L,L-1)\bm{q}_{L-1}(\bm{x})$ $\displaystyle=$
$\displaystyle\mathcal{G}(L,L-1)\cdots\mathcal{G}(1,0)\bm{q}_{0}(\bm{x}),$
(29)

i.e., as a forward propagation problem (see Figure 3). Note that we can write
(29) (or (25)) explicitly in terms of iterated integrals involving single
layer transition densities as

$\displaystyle\bm{q}_{L}(\bm{x})=$
$\displaystyle\int\bm{u}(\bm{y})p_{0|L}(\bm{y}|\bm{x})d\bm{y}$
$\displaystyle=$ $\displaystyle\int\bm{u}(\bm{y})\left(\int\cdots\int
p_{L|L-1}(\bm{y}|\bm{x}_{L-1})\cdots
p_{2|1}(\bm{x}_{2}|\bm{x}_{1})p_{1|0}(\bm{x}_{1}|\bm{x})d\bm{x}_{L-1}\cdots
d\bm{x}_{1}\right)d\bm{y}.$ (30)

Remark: We emphasize that the connection between continuous-time and discrete-
time composition operators can be rigorously established based on the concept
of suspension of a diffeomorphism. Specifically, it was shown in [37,
Proposition 3.7] that every diffeomorphism between compact manifolds is the
Poincaré map associated with a global transversal section of some vector
field. In other words, assuming that the continuous time conditional
expectation $\mathbb{E}\\{\bm{u}(\bm{X}(t))|\bm{X}(s)=\bm{x}\\}$ is a
diffeomorphism in $s$ and $t$, then there exists a continuous time-dynamical
system generating $\mathbb{E}\\{\bm{u}(\bm{X}(t))|\bm{X}(s)=\bm{x}\\}$ and a
global transversal section $\Sigma$ of the phase space such that
$\mathbb{E}\\{\bm{u}(\bm{X}_{L})|\bm{X}_{n}=\bm{x}\\}$ are the intersections
of the orbits $\mathbb{E}\\{\bm{u}(\bm{X}(t_{L}))|\bm{X}(s)=\bm{x}\\}$
(parameterized by $s$) with $\Sigma$. The composition operator for the
discrete dynamical system can be obtained by simply evaluating the continuous-
time composition operator at the appropriate “returning times” $t_{n}$.

Figure 3: Sketch of the forward/backward integration process for probability
density functions (PDFs) and conditional expectations. The transfer operator
$\mathcal{N}(n+1,n)$ maps the PDF $p_{n}(\bm{x})$ of the state $\bm{X}_{n}$
into $p_{n+1}(\bm{x})$ forward through the neural network. On the other hand,
the composition operator $\mathcal{M}$ maps the conditional expectation
$\mathbb{E}\left[\bm{u}(\bm{X}_{L})|\bm{X}_{n+1}=\bm{x}\right]$ backwards to
$\mathbb{E}\left[\bm{u}(\bm{X}_{L})|\bm{X}_{n}=\bm{x}\right]$. By defining the
operator $\mathcal{G}(n,m)=\mathcal{M}(L-n,L-m)$ we can transform the backward
propagation problem for
$\mathbb{E}\left[\bm{u}(\bm{X}_{L})|\bm{X}_{n}=\bm{x}\right]$ into a forward
propagation problem for
$\bm{q}_{n}(\bm{x})=\mathbb{E}\left[\bm{u}(\bm{X}_{L})|\bm{X}_{L-n}=\bm{x}\right]$.

###  3.3 Relation between composition and transfer operators

The integral operators $\mathcal{M}$ and $\mathcal{N}$ defined in (22) and
(16) involve the same kernel function, i.e., the multi-layer transition
probability density $p_{q|n}(\bm{x},\bm{y})$. In particular, we noticed that
$\mathcal{M}(n,q)$ integrates $p_{q|n}$ “from the left”, while
$\mathcal{N}(q,n)$ integrates it “from the right”. It is easy to show that
$\mathcal{M}(n,q)$ and $\mathcal{N}(q,n)$ are adjoint to each other relative
to the standard inner product in $L^{2}$ (see [9] for the continuous-time
case). In fact,

$\displaystyle\mathbb{E}\\{\bm{u}(\bm{X}_{k})\\}=$
$\displaystyle\int\mathbb{E}\\{\bm{u}(\bm{X}_{k})|\bm{X}_{q}=\bm{x}\\}p_{q}(\bm{x})d\bm{x}$
$\displaystyle=$
$\displaystyle\int\left[\mathcal{M}(q,j)\mathbb{E}\\{\bm{u}(\bm{X}_{k})|\bm{X}_{j}=\bm{x}\\}\right]p_{q}(\bm{x})d\bm{x}$
$\displaystyle=$
$\displaystyle\int\mathbb{E}\\{\bm{u}(\bm{X}_{k})|\bm{X}_{j}=\bm{x}\\}\mathcal{N}(j,q)p_{q}(\bm{x})d\bm{x}.$
(31)

Therefore

$\mathcal{M}(q,j)^{*}=\mathcal{N}(j,q)\qquad\forall q,j\in\\{0,\ldots,L\\},$
(32)

where $\mathcal{M}(q,j)^{*}$ denotes the operator adjoint of
$\mathcal{M}(q,j)$ with respect to the $L^{2}$ inner product. By invoking the
definition (28), we can also write (32) as

$\mathcal{G}(L-q,L-j)^{*}=\mathcal{N}(j,q),\qquad\forall
j,q\in\\{0,\ldots,L\\}.$ (33)

In A we show that if the cumulative distribution function of each random
vector $\bm{\xi}_{n}$ in the noise process has partial derivatives that are
Lipschitz in $\mathscr{R}(\bm{\xi}_{n})$ (range of $\bm{\xi}_{n}$), then the
composition and transfer operators defined in Eqs. (22) and 16 are bounded in
$L^{2}$ (see Proposition 16 and Proposition 17). Moreover, is possible to
choose the probability density of $\bm{\xi}_{n}$ such that the single layer
composition and transfer operators become strict contractions. This property
will be used in section 4.3 to prove that the memory of a stochastic neural
network driven by particular types of noise decays with the number of layers.

###  3.4 Multi-layer conditional transition density

We have seen that the composition and the transfer operators $\mathcal{M}$ and
$\mathcal{N}$ defined in Eqs. (22) and (16), allow us to push forward and
backward conditional expectations and probability densities across the entire
neural network. Moreover such operators are adjoint to one another (see
equation (32)) [9, 53, 3], and also have the same kernel, i.e., the transition
density $p_{n|q}(\bm{x}_{n}|\bm{x}_{q})$. In this section, we determine an
explicit expression for such transition density. To this end, we fist derive
analytical formulas for the one-layer transition density
$p_{n+1|n}(\bm{x}_{n+1}|\bm{x}_{n})$ for various types of neural network
models. The multi-layer transition density $p_{n|q}(\bm{x}_{n}|\bm{x}_{q})$ is
then obtained by composing one-layer transition densities as follows

$p_{n|q}(\bm{x}_{n}|\bm{x}_{q})=\int\cdots\int
p_{n|n-1}(\bm{x}_{n}|\bm{x}_{n-1})\cdots
p_{q+1|q}(\bm{x}_{q+1}|\bm{x}_{q})d\bm{x}_{n-1}\cdots d\bm{x}_{q+1}.$ (34)

Let us first consider the neural network model

$\bm{X}_{n+1}=\bm{F}(\bm{X}_{n},\bm{w}_{n})+\bm{\xi}_{n},$ (35)

where $\\{\bm{\xi}_{0},\bm{\xi}_{1},\ldots\\}$ is a discrete (vector-valued)
Markov process indexed by “$n$”. By (35), $\bm{X}_{n+1}$ is the sum of two
independent random vectors555Recall that $\bm{X}_{n}$ and $\bm{\xi}_{n}$ are
statistically independent random vectors. Hence,
$\bm{F}(\bm{X}_{n},\bm{w}_{n})$ and $\bm{\xi}_{n}$ are statistically
independent random vectors., i.e., $\bm{F}(\bm{X}_{n},\bm{w}_{n})$ and
$\bm{\xi}_{n}$. Given any measurable function $h(\bm{x})$ we clearly have

$\displaystyle\mathbb{E}\left\\{h(\bm{X}_{n+1})\right\\}=$ $\displaystyle\int
h(\bm{x})p_{n+1}(\bm{x})d\bm{x}$ $\displaystyle=$ $\displaystyle\int\int
h(\bm{F}(\bm{x},\bm{w}_{n})+\bm{\xi})p_{n}(\bm{x})\rho_{n}(\bm{\xi})d\bm{x}d\bm{\xi}$
$\displaystyle=$ $\displaystyle\int\int
h(\bm{y}))\underbrace{\rho_{n}(\bm{x}-\bm{F}(\bm{y},\bm{w}_{n}))}_{p_{n+1|n}(\bm{x}|\bm{y})}p_{n}(\bm{x})d\bm{y}d\bm{x},$
(36)

where $\rho_{n}(\bm{x})$ denotes the probability density of the random vector
$\bm{\xi}_{n}$. Therefore, the one-layer transition density for the neural
network model (35) is666Equation (39) can be derived in a more general setting
by recalling the conditional probability identity
$p_{n+1|n}(\bm{x}_{n+1}|\bm{x}_{n})=\int
p_{\bm{X}_{n+1}|\bm{X}_{n},\bm{\xi}_{n}}(\bm{x}_{n+1}|\bm{x}_{n},\bm{z})\rho_{n}(\bm{z})d\bm{z}.$
(37) The conditional density of $\bm{X}_{n+1}$ given $\bm{X}_{n}=\bm{x}_{n}$
and $\bm{\xi}_{n}=\bm{z}$, i.e.,
$p_{\bm{X}_{n+1}|\bm{X}_{n},\bm{\xi}_{n}}(\bm{x}_{n+1}|\bm{x}_{n},\bm{z})$,
can be immediately computed by using (35) as
$p_{\bm{X}_{n+1}|\bm{X}_{n},\bm{\xi}_{n}}(\bm{x}_{n+1}|\bm{x}_{n},\bm{z})=\delta(\bm{x}_{n+1}-\bm{F}(\bm{x}_{n},\bm{w}_{n})-\bm{z}),$
(38) where $\delta(\bm{x})$ is the multivariate Dirac delta function.
Substituting (38) into (37), and integrating over $\bm{z}$ yields (39).

$p_{n+1|n}(\bm{x}_{n+1}|\bm{x}_{n})=\rho_{n}(\bm{x}_{n+1}-\bm{F}(\bm{x}_{\bm{n}},\bm{w}_{n})).$
(39)

Note that such transition density depends on the PDF of the random noise
$\rho_{n}$, the neural activation function $\bm{F}$, and the neural network
weights $\bm{w}_{n}$. Next, consider the more general recurrent neural network
model (4), i.e.,

$\bm{X}_{n+1}=\bm{F}(\bm{X}_{n},\bm{w}_{n})+\bm{M}(\bm{X}_{n},\bm{\beta}_{n})\bm{\xi}_{n},$
(40)

which includes a multiplicative noise term
$\bm{M}(\bm{X}_{n},\bm{\beta}_{n})\bm{\xi}_{n}$. Recall that
$\bm{M}(\bm{X}_{n},\bm{\beta}_{n})$ is a $N\times m$ matrix, where $N$ (number
of neurons) can be larger or smaller than $m$ (dimension of each vector
$\bm{\xi}_{n}$). By following the same steps as before we obtain

$\displaystyle
p_{n+1|n}(\bm{x}_{n+1}|\bm{x}_{n})=\int\delta\left(\bm{x}_{n+1}-\bm{F}(\bm{x}_{n},\bm{w}_{n})-\bm{M}(\bm{x}_{n},\bm{\beta}_{n})\bm{\xi}\right)\rho_{n}(\bm{\xi})d\bm{\xi}.$
(41)

By using well-known properties of the Dirac delta function [22] it is possible
to re-write the function composition
$\delta\left(\bm{x}_{n+1}-\bm{F}(\bm{x}_{n},\bm{w}_{n})-\bm{M}(\bm{x}_{n},\bm{\beta}_{n})\bm{\xi}\right)$
in a way that makes the integral (41) easily computable. For instance, if
$\bm{M}(\bm{x}_{n},\bm{\beta}_{n})$ is a full rank square matrix ($m=N$) then

$\delta\left(\bm{x}_{n+1}-\bm{F}(\bm{x}_{n},\bm{w}_{n})-\bm{M}(\bm{x}_{n},\bm{\beta}_{n})\bm{\xi}\right)=\frac{1}{\left|\det(\bm{M}(\bm{x}_{n},\bm{\beta}_{n}))\right|}\delta\left(\bm{\xi}-\bm{M}(\bm{x}_{n})^{-1}\left[\bm{x}_{n+1}-\bm{F}(\bm{x}_{n},\bm{w}_{n})\right]\right),$
(42)

which yields

$p_{n+1|n}(\bm{x}_{n+1}|\bm{x}_{n})=\frac{1}{\left|\det(\bm{M}(\bm{x}_{n},\bm{\beta}_{n}))\right|}\rho_{n}\left(\bm{M}(\bm{x}_{n},\bm{\beta}_{n})^{-1}\left[\bm{x}_{n+1}-\bm{F}(\bm{x}_{n},\bm{w}_{n})\right]\right).$
(43)

Other cases where $\bm{M}$ is not square can be handled similarly [38, 22]. By
following the same procedure, it is possible also to derive the transition
density for the general neural network model (1).

Remark: The transition density (39) associated with the neural network model
(35) can be computed explicitly once we choose a probability model for
$\bm{\xi}_{n}\in\mathbb{R}^{N}$. For instance, if we assume that
$\\{\bm{\xi}_{0},\bm{\xi}_{1},\ldots,\\}$ are i.i.d. Gaussian random vectors
with PDF

$\rho_{n}(\bm{\xi})=\frac{1}{(2\pi)^{N/2}}e^{-\bm{\xi}^{T}\bm{\xi}/2}\qquad\text{for
all $n=0,\ldots,L$}$ (44)

then we can explicitly write the one-layer transition density (39) as

$p_{n+1|n}(\bm{x}_{n+1}|\bm{x}_{n})=\frac{1}{(2\pi)^{N/2}}\exp\left[-\frac{\left[\bm{x}_{n+1}-\bm{F}(\bm{x}_{\bm{n}},\bm{w}_{n})\right]^{T}\left[\bm{x}_{n+1}-\bm{F}(\bm{x}_{\bm{n}},\bm{w}_{n})\right]}{2}\right].$
(45)

In A we provide an analytical example of transition density for a neural
network with two layers (one neuron per layer), $\tanh(\cdot)$ activation
function, and uniformly distributed random noise.

###  3.5 The zero noise limit

An important question is what happens to the neural network as we send the
amplitude of the noise to zero. To answer this question consider the system
(35) and the introduce the parameter $\epsilon\geq 0$, i.e.,

$\bm{X}_{n+1}=\bm{F}(\bm{X}_{n},\bm{w}_{n})+\epsilon\bm{\xi}_{n},$ (46)

We are interested in studying the orbits of this system as
$\epsilon\rightarrow 0$. To this end, we assume the $\bm{\xi}_{n}$ to be
independent random vectors each having the same density $\rho(\bm{x})$. This
implies that for all $n=0,\ldots,L-1$, the PDF of $\epsilon\bm{\xi}_{n}$ is

$\epsilon\bm{\xi}_{n}\sim\frac{1}{\epsilon}\rho_{n}\left(\frac{\bm{x}}{\epsilon}\right).$
(47)

It is shown in [24, Proposition 10.6.1] that the operator $\mathcal{N}(n+1,n)$
defined in (16)

$\displaystyle p_{n+1}(\bm{x})=$
$\displaystyle\mathcal{N}(n+1,n)p_{n}(\bm{x})$ $\displaystyle=$
$\displaystyle\int\frac{1}{\epsilon}\rho_{n}\left(\frac{\bm{x}-\bm{F}(\bm{z},\bm{w}_{n})}{\epsilon}\right)p_{n}(\bm{z})d\bm{z}$
(48)

converges in norm to the Frobenious-Perron operator corresponding to
$\bm{F}(\bm{X}_{n},\bm{w}_{n})$ as $\epsilon\rightarrow 0$. Indeed, in the
limit $\epsilon\rightarrow 0$ we have, formally

$\lim_{\epsilon\rightarrow
0}p_{n|n+1}\left(\bm{x}_{n+1}|\bm{x}_{n}\right)=\lim_{\epsilon\rightarrow
0}\int\frac{1}{\epsilon}\rho_{n}\left(\frac{\bm{x}_{n+1}-\bm{F}(\bm{x}_{n},\bm{w}_{n})}{\epsilon}\right)=\delta\left(\bm{x}_{n+1}-\bm{F}(\bm{x}_{n},\bm{w}_{n})\right).$
(49)

Substituting this expression into (16), one gets,

$p_{n+1}(\bm{x})=\mathcal{N}(n+1,n)p_{n}(\bm{x})=\int\delta\big{(}\bm{x}-\bm{F}(\bm{z},\bm{w}_{n})\big{)}p_{n}(\bm{z})d\bm{z}.$
(50)

Similarly, a substitution into equation (29) yields

$\bm{q}_{n+1}(\bm{x})=\mathcal{G}(n+1,n)\bm{q}_{n}(\bm{x})=q_{n}\left(\bm{F}(\bm{x},\bm{w}_{L-n-1})\right),$
(51)

i.e, the familiar function composition representation of neural network
mappings

$\bm{q}_{n+1}=\bm{q}_{0}\Big{(}\bm{F}\big{(}\bm{F}(\cdots\bm{F}(\bm{x},\bm{w}_{L-n})\cdots,\bm{w}_{L-1}),\bm{w}_{L}\big{)}\Big{)}.$
(52)

##  4 The Mori-Zwanzig formulation of deep learning

In section 3 we defined two linear operators, i.e., $\mathcal{N}(n,q)$ and
$\mathcal{M}(n,q)$ (Eqs. (16) and (22)), mapping the probability density of
the state $\bm{X}_{n}$ and the conditional expectation of a phase space
function $\bm{u}(\bm{X}_{n})$ forward or backward across different layers of
the neural network. In particular, we have shown that

$\displaystyle p_{n+1}(\bm{x})$
$\displaystyle=\mathcal{N}(n+1,n)p_{n}(\bm{x}),$ (53)
$\displaystyle\mathbb{E}\\{\bm{u}(\bm{X}_{L})|\bm{X}_{n}=\bm{x}\\}$
$\displaystyle=\mathcal{M}(n,n+1)\mathbb{E}\\{\bm{u}(\bm{X}_{L})|\bm{X}_{n+1}=\bm{x}\\}.$
(54)

Equation (53) maps the PDF of the state $\bm{X}_{j}$ forward through the
neural network, i.e., from the input to the output as $n$ increases, while
(54) maps the conditional expectation backward. We have also shown in section
3.2 that upon definition of

$\bm{q}_{n}(\bm{x})=\mathbb{E}\\{\bm{u}(\bm{X}_{L})|\bm{X}_{L-n}=\bm{x}\\}$
(55)

we can rewrite (54) as a forward propagation problem, i.e.,

$\bm{q}_{n+1}(\bm{x})=\mathcal{G}(n+1,n)\bm{q}_{n}(\bm{x}),$ (56)

where $\mathcal{G}(n,q)=\mathcal{M}(L-n,L-q)$ and $\mathcal{M}$ is defined in
(22). The function $\bm{q}_{n}(\bm{x})$ is defined on the domain

$\mathscr{R}(\bm{X}_{L-n})=\\{\bm{X}_{L-n}(\omega)\in\mathbb{R}^{N}\,:\,\omega\in\mathcal{S}\\},$
(57)

i.e., on the range of the random variable $\bm{X}_{L-n}(\omega)$ (see
Definition (141)). $\mathscr{R}(\bm{X}_{L-n})$ is a deterministic subset of
$\mathbb{R}^{N}$.

Eqs. (54) constitute the basis for developing the Mori-Zwanzig (MZ)
formulation of deep neural networks. The MZ formulation is a technique
originally developed in statistical mechanics [35, 57] to formally integrate
out phase variables in nonlinear dynamical systems by means of a projection
operator. One of the main features of such formulation is that it allows us to
systematically derive exact equations for quantities of interest, e.g., low-
dimensional observables, based on the equations of motion of the full system.
In the context of deep neural networks such equations of motion are Eqs.
(53)-(54), and (56).

To develop the Mori-Zwanzig formulation of deep learning, we introduce a
layer-dependent orthogonal projection operator $\mathcal{P}_{n}$ together with
the complementary projection $\mathcal{Q}_{n}=\mathcal{I}-\mathcal{P}_{n}$.
The nature and properties of $\mathcal{P}_{n}$ will be discussed in detail in
section 5. For now, it suffices to assume only that $\mathcal{P}_{n}$ is a
self-adjoint bounded linear operator, and that
$\mathcal{P}_{n}^{2}=\mathcal{P}_{n}$, i.e., $\mathcal{P}_{n}$ is idempotent.
To derive the MZ equation for neural networks, let us consider a general
recursion,

$\bm{g}_{n+1}(\bm{x})=\mathcal{R}(n+1,n)\bm{g}_{n}(\bm{x}),$ (58)

where $\\{\bm{g}_{n},\mathcal{R}(n+1,n)\\}$ can be either
$\\{p_{n},\mathcal{N}(n+1,n)\\}$ or $\\{\bm{q}_{n},\mathcal{G}(n+1,n)\\}$,
depending on the context of the application.

###  4.1 The projection-first and propagation-first approaches

We apply the projection operators $\mathcal{P}_{n}$ and $\mathcal{Q}_{n}$ to
(58) to obtain the following coupled system of equations

$\displaystyle\bm{g}_{n+1}=$
$\displaystyle\mathcal{R}(n+1,n)\mathcal{P}_{n}\bm{g}_{n}+\mathcal{R}(n+1,n)\mathcal{Q}_{n}\bm{g}_{n},$
(59) $\displaystyle\mathcal{Q}_{n+1}\bm{g}_{n+1}=$
$\displaystyle\mathcal{Q}_{n+1}\mathcal{R}(n+1,n)\mathcal{P}_{n}\bm{g}_{n}+\mathcal{Q}_{n+1}\mathcal{R}(n+1,n)\mathcal{Q}_{n}\bm{g}_{n}.$
(60)

By iterating the difference equation (60), we obtain the following
formula777Note that the difference equation (60) can be written as
$\bm{h}_{n+1}=\mathcal{A}_{n}\bm{h}_{n}+\bm{c}_{n},$ (61) where
$\bm{h}_{n}=\mathcal{Q}_{n}\bm{g}_{n}$,
$\bm{c}_{n}=\mathcal{Q}_{n+1}\mathcal{R}(n+1,n)\mathcal{P}_{n}\bm{g}_{n}$, and
$\mathcal{A}_{n}=\mathcal{Q}_{n+1}\mathcal{R}(n+1,n)$. As is well-known, the
solution to (61) is
$\bm{h}_{n}=\prod_{k=0}^{n-1}\mathcal{A}_{k}\bm{h}_{0}+\sum_{j=0}^{n-1}\Phi(n,j+1)\bm{c}_{j},\quad\text{where}\quad\Phi(n,m)=\mathcal{A}_{n}\cdots\mathcal{A}_{m}.$
(62) A substitution of $\mathcal{A}_{n}$, $\bm{h}_{n}$ and $\bm{c}_{j}$ into
(62) yields (63). for $\mathcal{Q}_{n}\bm{g}_{n}$

$\mathcal{Q}_{n}\bm{g}_{n}=\Phi_{\mathcal{R}}(n,0)\mathcal{Q}_{0}\bm{g}_{0}+\sum_{m=0}^{n-1}\Phi_{\mathcal{R}}(n,m)\mathcal{P}_{m}\bm{g}_{m},$
(63)

where $\Phi_{\mathcal{R}}(n,m)$ is the (forward) propagator of the orthogonal
dynamics, i.e.,

$\Phi_{\mathcal{R}}(n,m)=\mathcal{Q}_{n}\mathcal{R}(n,n-1)\cdots\mathcal{Q}_{m+1}\mathcal{R}(m+1,m).$
(64)

Since $\bm{g}_{0}$ is arbitrary, one can deduce the operator identity

$\mathcal{Q}_{n}\mathcal{R}(n,0)=\Phi_{\mathcal{R}}(n,0)\mathcal{Q}_{0}+\sum_{m=0}^{n-1}\Phi_{\mathcal{R}}(n,m)\mathcal{P}_{m}\mathcal{R}(m,0).$
(65)

A substitution of (63) into (59) yields the Mori-Zwanzig equation

$\bm{g}_{n+1}=\underbrace{\mathcal{R}(n+1,n)\mathcal{P}_{n}\bm{g}_{n}}_{\text{streaming}}+\underbrace{\mathcal{R}(n+1,n)\Phi_{\mathcal{R}}(n,0)\mathcal{Q}_{0}\bm{g}_{0}}_{\text{noise}}+\underbrace{\mathcal{R}(n+1,n)\sum_{m=0}^{n-1}\Phi_{\mathcal{R}}(n,m)\mathcal{P}_{m}\bm{g}_{m}}_{\text{memory}}.$
(66)

We call the first term at the right hand side of (66) streaming (or Markovian)
term. Such term reflects the change of $\bm{g}_{n}$ as we go from one layer to
the next. The third term represents the memory of the neural network, and it
encodes the interaction between the projected dynamics and its entire history.

The second term is the contribution from the orthogonal dynamics and it is
often regarded as noise. Note that if $\bm{g}_{0}$ is in the image of
$\mathcal{P}_{0}$, i.e., if $\mathcal{P}_{0}\bm{g}_{0}=\bm{g}_{0}$, then the
second term drops out, yielding a simplified MZ equation,

$\bm{g}_{n+1}=\mathcal{R}(n+1,n)\mathcal{P}_{n}\bm{g}_{n}+\mathcal{R}(n+1,n)\sum_{m=0}^{n-1}\Phi_{\mathcal{R}}(n,m)\mathcal{P}_{m}\bm{g}_{m}.$
(67)

To integrate (67) forward, i.e., from one layer to the next, we first project
$\bm{g}_{m}$ using $\mathcal{P}_{m}$ (for $m=0,\ldots,n$), then apply the
evolution operator $\mathcal{R}(n+1,n)$ to $\mathcal{P}_{n}\bm{g}_{n}$, and
the memory operator $\Phi_{\mathcal{R}}$ to the entire history of $\bm{g}_{m}$
(memory of the network). It is also possible to construct an MZ equation based
on the reversed mechanism, i.e., by projecting $\mathcal{R}(n+1,n)\bm{g}_{n}$
rather than $\bm{g}_{m}$. To this end, we first write (59) as

$\bm{g}_{n+1}=\mathcal{P}_{n+1}\mathcal{R}(n+1,n)\bm{g}_{n}+\mathcal{Q}_{n+1}\mathcal{R}(n+1,n)\bm{g}_{n},$
(68)

i.e., propagation via $\mathcal{R}(n+1,n)$ precedes projection (propagation-
first approach). By applying the variation of constant formula (65) to (68) we
arrive at a slightly different form of the MZ equation, namely

$\bm{g}_{n+1}=\mathcal{P}_{n+1}\mathcal{R}(n+1,n)\bm{g}_{m}+\Phi_{\mathcal{R}}(n+1,0)\bm{g}_{0}+\sum_{m=0}^{n-1}\Phi_{\mathcal{R}}(n+1,m+1)\mathcal{P}_{m+1}\mathcal{R}(m+1,m)\bm{g}_{m}.$
(69)

###  4.2 Discrete Dyson’s identity

Another form of the MZ equation (66) can be derived based on a discrete
version of the Dyson identity888For continuous-time autonomous dynamical
systems the Dyson’s identity can be written as [53, 54, 56, 48, 5]
$e^{t\mathcal{L}}=e^{t\mathcal{Q}\mathcal{L}}+\int_{0}^{t}e^{(t-s)\mathcal{L}}\mathcal{P}\mathcal{L}e^{s\mathcal{Q}\mathcal{L}}ds.$
(70) where $\mathcal{L}$ is the (time-independent) Liouvillian of the system.
The discrete Dyson identity and the corresponding discrete MZ formulation was
first derived by Lin and Lu in [29] for autonomous discrete dynamical
systems.. To derive such identity, consider the sequence

$\displaystyle\bm{y}_{n+1}=$
$\displaystyle\mathcal{Q}_{n+1}\mathcal{R}(n+1,n)\bm{y}_{n}$ (71)
$\displaystyle=$
$\displaystyle\mathcal{R}(n+1,n)\bm{y}_{n}-\mathcal{P}_{n+1}\mathcal{R}(n+1,n)\bm{y}_{n}.$
(72)

By using the discrete variation of constant formula, we can rewrite (72) as

$\bm{y}_{n}=\mathcal{R}(n,0)\bm{y}_{0}-\sum_{m=0}^{n-1}\mathcal{R}(n,m+1)\mathcal{P}_{m+1}\mathcal{R}(m+1,m)\bm{y}_{m}.$
(73)

Similarly, solving (71) yields

$\bm{y}_{n}=\Phi_{\mathcal{R}}(n,0)\bm{y}_{0},$ (74)

where $\Phi_{\mathcal{R}}$ is defined in (64). By substituting (74) into (73)
for both $\bm{y}_{n}$ and $\bm{y}_{m}$, and observing that $\bm{y}_{0}$ is
arbitrary, we obtain

$\mathcal{R}(n,0)=\Phi_{\mathcal{R}}(n,0)+\sum_{m=0}^{n-1}\mathcal{R}(n,m+1)\mathcal{P}_{m+1}\mathcal{R}(m+1,m)\Phi_{\mathcal{R}}(m,0).$
(75)

The operator identity (75) is the discrete version of the well-known
continuous-time Dyson’s identity. A substitution of (75) into
$\bm{g}_{n}=\mathcal{R}(n,0)\bm{g}_{0}$ yields the following alternative form
of the MZ equation (66)

$\bm{g}_{n+1}=\mathcal{P}_{n+1}\mathcal{R}(n+1,n)\Phi_{\mathcal{R}}(n,0)\bm{g}_{0}+\Phi_{\mathcal{R}}(n+1,0)\bm{g}_{0}+\sum_{m=0}^{n-1}\mathcal{R}(n+1,m+1)\mathcal{P}_{m+1}\mathcal{R}(m+1,m)\Phi_{\mathcal{R}}(m,0)\bm{g}_{0}.$
(76)

Here we have arranged the terms in the same way as in (66).

###  4.3 Mori-Zwanzig equations for probability density functions

The evolution of probability densities from one layer to the other across the
neural network is governed by the operator $\mathcal{N}(n,q)$ in (16).
Replacing $\mathcal{R}$ with $\mathcal{N}$ in (66) yields the following Mori-
Zwanzig equation for the PDF of $\bm{X}_{n}$

$p_{n+1}=\underbrace{\mathcal{N}(n+1,n)\mathcal{P}_{n}p_{n}}_{\text{streaming
term}}+\underbrace{\mathcal{N}(n+1,n)\Phi_{\mathcal{N}}(n,0)\mathcal{Q}_{0}p_{0}}_{\text{noise
term}}+\underbrace{\mathcal{N}(n+1,n)\sum_{m=0}^{n-1}\Phi_{\mathcal{N}}(n,m)\mathcal{P}_{m}p_{m}}_{\text{memory
term}},$ (77)

Alternatively, by using the MZ equation (76), we can write

$p_{n}=\Phi_{\mathcal{N}}(n,0)p_{0}+\sum_{m=0}^{n-1}\mathcal{N}(n,m+1)\mathcal{P}_{m+1}\mathcal{N}(m+1,m)\Phi_{\mathcal{N}}(m,0)p_{0},$
(78)

where

$\Phi_{\mathcal{N}}(n,m)=\mathcal{Q}_{n}\mathcal{N}(n,n-1)\cdots\mathcal{Q}_{m+1}\mathcal{N}(m+1,m).$
(79)

###  4.4 Mori-Zwanzig equation for conditional expectations

Next, we discuss MZ equations in neural nets propagating conditional
expectations

$\bm{q}_{n}(\bm{x})=\mathbb{E}\\{\bm{u}(\bm{X}_{L})|\bm{X}_{L-n}=\bm{x}\\}$
(80)

backward across the network, i.e., from $\bm{q}_{0}(\bm{x})=\bm{u}(\bm{x})$
into
$\bm{q}_{L}(\bm{x})=\mathbb{E}\\{\bm{u}(\bm{X}_{L})|\bm{X}_{0}=\bm{x}\\}$. To
simplify the notation, we denote the projection operators in the space of
conditional expectations with the same letters as in the space of PDFs, i.e.,
$\mathcal{P}_{n}$ and $\mathcal{Q}_{n}$999 The orthogonal projection for
conditional expectations is the operator adjoint of the projection
$\mathcal{P}_{m}$ that operates on probability densities, i.e.,
$\int\mathbb{E}\\{\bm{u}(\bm{X}_{k})|\bm{X}_{q}=\bm{x}\\}\mathcal{P}_{q}p_{q}(\bm{x})d\bm{x}=\int\mathcal{P}_{q}^{*}\mathbb{E}\\{\bm{u}(\bm{X}_{k})|\bm{X}_{q}=\bm{x}\\}p_{q}(\bm{x})d\bm{x}.$
(81) Such adjoint relation is the same that connects the composition and
transfer operators ($\mathcal{M}(q,n)$ and $\mathcal{N}(n,q)$ in equation
(32)). The connection between projections for probability densities and
conditional expectations was extensively discussed in [9] in the setting of
operator algebras.. Replacing $\mathcal{R}$ with $\mathcal{G}$ in (66) yields
the following MZ equation for the conditional expectations

$\bm{q}_{n+1}=\underbrace{\mathcal{G}(n+1,n)\mathcal{P}_{n}\bm{q}_{n}}_{\text{streaming
term}}+\underbrace{\mathcal{G}(n+1,n)\Phi_{\mathcal{G}}(n,0)\mathcal{Q}_{0}\bm{q}_{0}}_{\text{noise
term}}+\underbrace{\mathcal{G}(n+1,n)\sum_{m=0}^{n-1}\Phi_{\mathcal{G}}(n,m)\mathcal{P}_{m}\bm{q}_{m}}_{\text{memory
term}},$ (82)

where

$\Phi_{\mathcal{G}}(n,m)=\mathcal{Q}_{n}\mathcal{G}(n,n-1)\cdots\mathcal{Q}_{m+1}\mathcal{G}(m+1,m).$
(83)

Equation (82) can be equivalently written by incorporating the streaming term
into the summation of the memory term

$\bm{q}_{n+1}=\mathcal{G}(n+1,n)\Phi_{\mathcal{G}}(n,0)\mathcal{Q}_{0}\bm{q}_{0}+\mathcal{G}(n+1,n)\sum_{m=0}^{n}\Phi_{\mathcal{G}}(n,m)\mathcal{P}_{m}\bm{q}_{m}.$
(84)

An alternative formulation of the MZ equation which follows from Eq. (76) is

$\bm{q}_{n}=\Phi_{\mathcal{G}}(n,0)\bm{q}_{0}+\sum_{m=0}^{n-1}\mathcal{G}(n,m+1)\mathcal{P}_{m+1}\mathcal{G}(m+1,m)\Phi_{\mathcal{G}}(m,0)\bm{q}_{0}.$
(85)

Remark: The Mori-Zwanzig equations (77)-(78) and (82)-(85) allow us to perform
dimensional reduction within each layer of the network (number of neurons per
layer, via projection), or across different layers (total number of layers,
via memory approximation). The MZ formulation is also useful to perform
theoretical analysis of deep learning by using tools from operator theory. As
we shall see in the next section, the memory of the neural network can be
controlled by controlling the discrete Markov process
$\\{\bm{\xi}_{0},\bm{\xi}_{1},\ldots\\}$.

##  5 Mori-Zwanzig projection operator

Suppose that the neural network model (6) is perturbed by a discrete Markov
process $\\{\bm{\xi}_{n}\\}$ with bounded range. In this hypothesis we have
shown in A that the range of each random vector $\bm{X}_{m}$ is bounded. With
specific reference to MZ equations for scalar conditional expectations (i.e.,
conditional averages of scalar quantities of interest)

$q_{m}(\bm{x})=\mathbb{E}\left[u(\bm{X}_{L})|\bm{X}_{L-m}=\bm{x}\right]$

and recalling that101010The operator $\mathcal{M}(L-m,L-m-1)$ defined in (22)
is an operator from $L^{2}(\mathscr{R}(\bm{X}_{L-m+1}))$ onto
$L^{2}(\mathscr{R}(\bm{X}_{L-m}))$. Hence we define $\mathcal{P}_{m}$ in (87)
as an operator on $L^{2}(\mathscr{R}(\bm{X}_{L-m}))$.

$\displaystyle q_{m}(\bm{x})=$ $\displaystyle\mathcal{G}(m,m-1)q_{m-1}$
$\displaystyle=$ $\displaystyle\mathcal{M}(L-m,L-m+1)q_{m-1}(\bm{x})$
$\displaystyle=$
$\displaystyle\int_{\mathscr{R}(\bm{X}_{L-m+1})}\underbrace{p_{L-m+1|L-m}(\bm{y},\bm{x})}_{\rho_{L-m}(\bm{y}-\bm{F}\left(\bm{x},\bm{w}_{L-m})\right)}q_{m-1}(\bm{y})d\bm{y},$
(86)

we define the following orthogonal projection operator111111The projection
operator (88) can be extended to vector-valued functions and conditional
expectations by defining an appropriate matrix-valued kernel
$\bm{K}(\bm{x},\bm{y})$.

$\displaystyle\mathcal{P}_{m}:L^{2}(\mathscr{R}(\bm{X}_{L-m})$
$\displaystyle)\mapsto L^{2}(\mathscr{R}(\bm{X}_{L-m}))$ (87) $\displaystyle
f\mapsto\int_{\mathscr{R}(\bm{X}_{L-m})}K_{L-m}(\bm{x},\bm{y})f(\bm{y})d\bm{y}.$
(88)

Since $\mathcal{P}_{m}$ is, by definition, an orthogonal projection we have
that $\mathcal{P}_{m}$ is idempotent ($\mathcal{P}_{m}^{2}=\mathcal{P}_{m}$),
bounded and self-adjoint relative to the inner product in
$L^{2}(\mathscr{R}(\bm{X}_{L-m}))$. These conditions imply that the kernel
$K_{L-m}(\bm{x},\bm{y})$ is a symmetric Hilbert-Schmidt kernel that satisfies
the reproducing kernel condition

$\int_{\mathscr{R}(\bm{X}_{L-m})}K_{L-m}(\bm{x},\bm{y})K_{L-m}(\bm{y},\bm{z})d\bm{y}=K_{L-m}(\bm{x},\bm{z}),\qquad\forall\bm{x},\bm{y}\in\mathscr{R}(\bm{X}_{L-m}).$
(89)

Note that the classical Mori’s projection [54, 53] can be written in the form
(88) if we set

$K_{L-m}(\bm{x},\bm{y})=\sum_{i=0}^{M}\eta^{m}_{i}(\bm{x})\eta^{m}_{i}(\bm{y}),$
(90)

where $\\{\eta^{m}_{0},\ldots,\eta^{m}_{M}\\}$ are orthonormal functions in
$L^{2}\left(\mathscr{R}(\bm{X}_{L-m})\right)$. In A we show that for general
noises $\\{\bm{\xi}_{m}\\}$, the range of $\bm{X}_{L-m}$ can vary across
layers. Therefore the choice of the orthonormal functions
$\\{\eta^{m}_{j}(\bm{x})\\}$ depends on the neural network layer “$m$”.

The projection operator $\mathcal{P}_{m}$ is said to be non-negative if for
all positive $v(\bm{x})\in L^{2}(\mathscr{R}(\bm{X}_{L-m}))$ ($v>0$) we have
that $\mathcal{P}_{m}v\geq 0$ [20]. Clearly, this implies that the kernel
$K_{L-m}(\bm{x},\bm{y})$ is non-negative in
$\mathscr{R}(\bm{X}_{L-m})\times\mathscr{R}(\bm{X}_{L-m})$ [12]. An example of
a kernel defining a non-negative orthogonal projection is

$K_{L-m}(\bm{x},\bm{y})=\eta^{m}(\bm{x})\eta^{m}(\bm{y}),\qquad\eta^{m}(\bm{x})\geq
0,\qquad\left\|\eta^{m}\right\|_{L^{2}(\mathscr{R}(\bm{X}_{L-m}))}=1.$ (91)

More generally, if $K_{L-m}(\bm{x},\bm{y})$ is any square-integrable symmetric
conditional probability density function on
$\mathscr{R}(\bm{X}_{L-m})\times\mathscr{R}(\bm{X}_{L-m})$, then
$\mathcal{P}_{m}$ is a non-negative orthogonal projection.

###  5.1 Analysis of the operator $\mathcal{Q}_{m}\mathcal{G}(m,m-1)$ and
fading memory property

The operator $\mathcal{Q}_{m}\mathcal{G}(m,m-1)$ plays a very important role
in the MZ equation (85). In fact, the memory operator $\Phi_{\mathcal{G}}$
defined in (83) is essentially an operator product involving multiple
$\mathcal{Q}_{m}\mathcal{G}(m,m-1)$. In this section, we aim at determining
conditions on
$Q_{m}\mathcal{G}(m,m-1)=(\mathcal{I}-\mathcal{P}_{m})\mathcal{G}(m,m-1)$,
e.g., noise level and distribution, such that

$\left\|\mathcal{Q}_{m}\mathcal{G}(m,m-1)\right\|=\sup_{\begin{subarray}{c}v\in
L^{2}(\mathscr{R}(\bm{X}_{L-m+1}))\\\ v\neq
0\end{subarray}}\frac{\left\|\mathcal{Q}_{m}\mathcal{G}(m,m-1)v\right\|_{L^{2}(\mathscr{R}(\bm{X}_{L-m}))}}{\left\|v\right\|_{L^{2}(\mathscr{R}(\bm{X}_{L-m+1}))}}<1.$
(92)

In this way, $\mathcal{Q}_{m}\mathcal{G}(m,m-1)$ becomes a contraction in
$L^{2}$, and therefore the MZ memory term in (85) decays with the number of
layers, while the first term (i.e., $\Phi_{G}(L,0)\bm{q}_{0}$) decays to zero.
Indeed, if (92) holds true, then the norm of memory operator
$\Phi_{\mathcal{G}}(n,m)$ defined in (83) decays with the number of
“$\mathcal{Q}\mathcal{G}$” operator products taken. For a neural network with
$L$ layers, we can formally write the operator
$\mathcal{Q}_{m}\mathcal{G}(m,m-1)$ as

$\displaystyle(\mathcal{Q}_{m}\mathcal{G}(m,m-1)v)(\bm{x})$
$\displaystyle=\mathcal{Q}_{m}\int_{\mathscr{R}(\bm{X}_{L-m+1})}\underbrace{\rho_{L-m}\left(\bm{y}-\bm{F}(\bm{x},\bm{w}_{L-m})\right)}_{p_{L-m+1|L-m}(\bm{y}|\bm{x})}v(\bm{y})d\bm{y}$
$\displaystyle=\int_{\mathscr{R}(\bm{X}_{L-m+1})}\gamma_{L-m}(\bm{y},\bm{x})v(\bm{y})d\bm{y},$
(93)

where

$\gamma_{L-m}(\bm{y},\bm{x})=\rho_{L-m}(\bm{y}-\bm{F}(\bm{x},\bm{w}_{L-m}))-\int_{\mathscr{R}(\bm{X}_{L-m})}K_{L-m}(\bm{x},\bm{z})\rho_{L-m}(\bm{y}-\bm{F}(\bm{z},\bm{w}_{L-m}))d\bm{z}.$
(94)

Recall that the conditional density
$p_{L-m+1|L-m}(\bm{y}|\bm{x})=\rho_{L-m}(\bm{y}-\bm{F}(\bm{x},\bm{w}_{L-m}))$
is defined on the set $\mathscr{B}_{L-m}$ (see Eq. (151)). In the zero noise
case, equation (93) reduces to

$(\mathcal{Q}_{m}\mathcal{G}(m,m-1)v)(\bm{x})=v\left(\bm{F}(\bm{x},\bm{w}_{L-m})\right)-\int_{\mathscr{R}(\bm{X}_{L-m})}K_{L-m}(\bm{x},\bm{z})v\left(\bm{F}(\bm{z},\bm{w}_{L-m})\right)d\bm{z}.$
(95)

Let us determine conditions on the projection kernel $K_{L-m}$ that guarantee
that $\mathcal{Q}_{m}\mathcal{G}(m,m-1)$ is a contraction independently of the
neural network weights. To this end, let us assume that $K_{L-m}$ is an
element of $L^{2}(\mathscr{R}(\bm{X}_{L-m})\times\mathscr{R}(\bm{X}_{L-m}))$
and expand it as121212As is well known, if $K(\bm{x},\bm{y})$ is a (symmetric)
bounded projection kernel satisfying (89) then $K$ is necessarily separable,
i.e., it can be written in the form (96).

$K_{L-m}(\bm{x},\bm{y})=c_{m}+\sum_{k=1}^{M}\eta^{m}_{i}(\bm{x})\eta^{m}_{i}(\bm{y}),$
(96)

where $c_{m}$ is a real number and $\eta^{m}_{i}$ are zero-mean orthonormal
basis functions, i.e.,

$\int_{\mathscr{R}(\bm{X}_{L-m})}\eta^{m}_{i}(\bm{x})d\bm{x}=0,\qquad\int_{\mathscr{R}(\bm{X}_{L-m})}\eta^{m}_{i}(\bm{x})\eta^{m}_{j}(\bm{x})d\bm{x}=\delta_{ij}.$
(97)

######  Lemma 1.

The kernel (96) satisfies the idempotency requirement (89) if and only if

$c_{m}=0\qquad\text{or}\qquad
c_{m}=\frac{1}{\lambda(\mathscr{R}(\bm{X}_{L-m}))},$ (98)

where $\lambda(\mathscr{R}(\bm{X}_{L-m}))$ denotes the Lebesgue measure of the
set $\mathscr{R}(\bm{X}_{L-m})$.

###### Proof.

By substituting (96) into (89) and taking into account (97) we obtain

$c_{m}^{2}\lambda(\mathscr{R}(\bm{X}_{L-m}))=c_{m},$ (99)

from which we obtain $c_{m}=0$ or
$c_{m}=1/\lambda(\mathscr{R}(\bm{X}_{L-m}))$.

∎

Clearly, if $\mathcal{G}(m,m-1)$ is itself a contraction and $\mathcal{Q}_{m}$
is an orthogonal projection, then the operator product
$\mathcal{Q}_{m}\mathcal{G}(m,m-1)$ is automatically a contraction. In the
following Lemma we compute a bound for the operator norm of
$\mathcal{Q}_{m}\mathcal{G}(m,m-1)$.

######  Proposition 2.

Let $\mathcal{Q}_{m}$ be an orthogonal projection. Define

$\left\|\rho_{L-m}\right\|_{L^{2}(\mathscr{R}(\bm{\xi}_{L-m}))}^{2}=\int_{\mathscr{R}(\bm{\xi}_{L-m})}\rho_{L-m}(\bm{x})^{2}d\bm{x}.$
(100)

Then under the same hypotheses of Proposition 17 in A we have that

$\left\|\mathcal{Q}_{m}\mathcal{G}(m,m-1)\right\|^{2}\leq\lambda(\Omega_{L-m})\left\|\rho_{L-m}\right\|^{2}_{L^{2}(\mathscr{R}(\bm{\xi}_{L-m}))},$
(101)

where $\lambda(\Omega_{L-m})$ is the Lebesgue measure of the set
$\Omega_{L-m}$ defined in (140). In particular, if $\mathcal{G}(m,m-1)$ is a
contraction then $\mathcal{Q}_{m}\mathcal{G}(m,m-1)$ is a contraction.

###### Proof.

The last statement in the Lemma is trivial. In fact, if $\mathcal{Q}_{m}$ is
an orthogonal projection then its operator norm is less or equal to one.
Hence,

$\left\|Q_{m}\mathcal{G}(m,m-1)\right\|^{2}\leq\underbrace{\left\|Q_{m}\right\|^{2}}_{\leq
1}\left\|\mathcal{G}(m,m-1)\right\|^{2}.$ (102)

Therefore if $\mathcal{G}(m,m-1)$ is a contraction and $\mathcal{Q}_{m}$ is an
orthogonal projection then $\mathcal{Q}_{m}\mathcal{G}(m,m-1)$ is a
contraction. The upper bound in (101) follows from Proposition 17 applied to

$\mathcal{G}(m,m-1)=\mathcal{M}(L-m+1,L-m).$

∎

The upper bound in (101) can be slightly improved using the definition of the
projection kernel $K_{L-m}$. This is stated in the following Theorem.

######  Theorem 3.

Suppose that the projection kernel $K_{L-m}$ in (88) has the form (96) with
$c_{m}=1/\lambda(\mathscr{R}(\bm{X}_{L-m}))$. Then the operator norm of
$\mathcal{Q}_{m}\mathcal{G}(m,m-1)$ can be bounded as

$\left\|\mathcal{Q}_{m}\mathcal{G}(m,m-1)\right\|^{2}\leq\lambda(\Omega_{L-m})\left(\left\|\rho_{L-m}\right\|_{L^{2}(\mathscr{R}(\bm{\xi}_{L-m}))}^{2}-\frac{1}{\lambda(\Omega_{L-m+1})}\right).$
(103)

The upper bound in (103) is independent of the neural network weights.

###### Proof.

The function $\gamma_{L-m}(\bm{y},\bm{x})$ defined in (94) is a Hilbert-
Schmidt kernel. Therefore, similarly to (157) we have

$\left\|\mathcal{Q}_{m}\mathcal{G}(m,m-1)\right\|^{2}\leq\left\|\gamma_{L-m}\right\|^{2}_{L^{2}(\mathscr{R}(\bm{X}_{L-m+1})\times\mathscr{R}(\bm{X}_{L-m}))}.$
(104)

The $L^{2}$ norm of $\gamma_{L-m}$ can be written as (see (94))

$\displaystyle\int_{\mathscr{R}(\bm{X}_{L-m+1})}\int_{\mathscr{R}(\bm{X}_{L-m})}\gamma_{L-m}(\bm{y},\bm{x})^{2}d\bm{y}d\bm{x}=\int_{\mathscr{R}(\bm{X}_{L-m+1})}\int_{\mathscr{R}(\bm{X}_{L-m})}\rho_{L-m}(\bm{y}-\bm{F}(\bm{x},\bm{w}_{L-m}))^{2}d\bm{y}d\bm{x}+$
$\displaystyle\int_{\mathscr{R}(\bm{X}_{L-m+1})}\int_{\mathscr{R}(\bm{X}_{L-m})}\left(\int_{\mathscr{R}(\bm{X}_{L-m})}K_{L-m}(\bm{x},\bm{z})\rho_{L-m}(\bm{y}-\bm{F}(\bm{z},\bm{w}_{L-m}))d\bm{z}\right)^{2}d\bm{y}d\bm{x}-$
$\displaystyle
2\int_{\mathscr{R}(\bm{X}_{L-m+1})}\int_{\mathscr{R}(\bm{X}_{L-m})}\rho_{L-m}(\bm{y}-\bm{F}(\bm{x},\bm{w}_{L-m}))\left(\int_{\mathscr{R}(\bm{X}_{L-m})}K_{L-m}(\bm{x},\bm{z})\rho_{L-m}(\bm{y}-\bm{F}(\bm{z},\bm{w}_{L-m}))d\bm{z}\right)d\bm{y}d\bm{x}.$
(105)

In equation (161) we show that the first term at the right hand side of (105)
can be written as131313Note that the PDF $\rho_{L-m}$ is nonzero only within
the set defined in (151).

$\displaystyle\int_{\mathscr{R}(\bm{X}_{L-m+1})}\int_{\mathscr{R}(\bm{X}_{L-m})}\rho_{L-m}(\bm{y}-\bm{F}(\bm{x},\bm{w}_{L-m}))^{2}d\bm{y}d\bm{x}=$
$\displaystyle\lambda(\mathscr{R}(\bm{X}_{L-m}))\left\|\rho_{L-m}\right\|_{L^{2}(\mathscr{R}(\bm{\xi}_{L-m}))}^{2}$
$\displaystyle\leq$
$\displaystyle\lambda(\Omega_{L-m})\left\|\rho_{L-m}\right\|_{L^{2}(\mathscr{R}(\bm{\xi}_{L-m}))}^{2}.$
(106)

The $L^{2}$ norm of $\rho_{L-m}$ is independent of the neural network weights
and it can be calculated analytically (if $\rho_{L-m}^{2}$ is easily
integrable), numerically, or bounded as in Lemma 18. A substitution of the
series expansion (96) in the second term at the right hand side of (105)
yields

$\displaystyle\int_{\mathscr{R}(\bm{X}_{L-m+1})}\int_{\mathscr{R}(\bm{X}_{L-m})}\left(\int_{\mathscr{R}(\bm{X}_{L-m})}K_{L-m}(\bm{x},\bm{z})\rho_{L-m}(\bm{y}-\bm{F}(\bm{z},\bm{w}_{L-m}))d\bm{z}\right)^{2}d\bm{y}d\bm{x}=$
$\displaystyle\frac{1}{\lambda(\mathscr{R}(\bm{X}_{L-m}))}\int_{\mathscr{R}(\bm{X}_{L-m+1})}\left(\int_{\mathscr{R}(\bm{X}_{L-m})}\rho_{L-m}(\bm{y}-\bm{F}(\bm{x},\bm{w}_{L-m}))d\bm{z}\right)^{2}d\bm{y}+$
$\displaystyle\sum_{k=1}^{M}\int_{\mathscr{R}(\bm{X}_{L-m+1})}\left(\int_{\mathscr{R}(\bm{X}_{L-m})}\rho_{L-m}(\bm{y}-\bm{F}(\bm{z},\bm{w}_{L-m}))\eta^{m}_{k}(\bm{z})d\bm{z}\right)^{2}d\bm{y}.$
(107)

Here, we used the fact that the basis functions $\eta^{m}_{k}(\bm{x})$ are
zero-mean and orthonormal in $\mathscr{R}(\bm{X}_{L-m})$ (see Eq. (97)).
Similarly, a substitution of the expansion (96) in the third term at the right
hand side of (105) yields

$\displaystyle
2\int_{\mathscr{R}(\bm{X}_{L-m+1})}\int_{\mathscr{R}(\bm{X}_{L-m})}\rho_{L-m}(\bm{y}-\bm{F}(\bm{x},\bm{w}_{L-m}))\left(\int_{\mathscr{R}(\bm{X}_{L-m})}K_{L-m}(\bm{x},\bm{z})\rho_{L-m}(\bm{y}-\bm{F}(\bm{z},\bm{w}_{L-m}))d\bm{z}\right)d\bm{y}d\bm{x}$
$\displaystyle=\frac{2}{\lambda(\mathscr{R}(\bm{X}_{L-m}))}\int_{\mathscr{R}(\bm{X}_{L-m+1})}\left(\int_{\mathscr{R}(\bm{X}_{L-m})}\rho_{L-m}(\bm{y}-\bm{F}(\bm{z},\bm{w}_{L-m}))d\bm{z}\right)^{2}d\bm{y}$
$\displaystyle+2\sum_{k=1}^{M}\int_{\mathscr{R}(\bm{X}_{L-m+1})}\left(\int_{\mathscr{R}(\bm{X}_{L-m})}\rho_{L-m}(\bm{y}-\bm{F}(\bm{z},\bm{w}_{L-m}))\eta^{m}_{k}(\bm{z})d\bm{z}\right)^{2}d\bm{y}.$
(108)

Combining (105)-(108) yields

$\displaystyle\int_{\mathscr{R}(\bm{X}_{L-m+1})}\int_{\mathscr{R}(\bm{X}_{L-m})}\gamma_{L-m}(\bm{y},\bm{x})^{2}d\bm{y}d\bm{x}=\lambda(\mathscr{R}(\bm{X}_{L-m}))\left\|\rho_{L-m}\right\|_{L^{2}(\mathscr{R}(\bm{\xi}_{L-m})}^{2}-$
$\displaystyle\frac{1}{\lambda(\mathscr{R}(\bm{X}_{L-m}))}\int_{\mathscr{R}(\bm{X}_{L-m+1})}\left(\int_{\mathscr{R}(\bm{X}_{L-m})}\rho_{L-m}(\bm{y}-\bm{F}(\bm{z},\bm{w}_{L-m}))d\bm{z}\right)^{2}d\bm{y}-$
$\displaystyle\sum_{k=1}^{M}\int_{\mathscr{R}(\bm{X}_{L-m+1})}\left(\int_{\mathscr{R}(\bm{X}_{L-m})}\rho_{L-m}(\bm{y}-\bm{F}(\bm{z},\bm{w}_{L-m}))\eta^{m}_{k}(\bm{z})d\bm{z}\right)^{2}d\bm{y}.$
(109)

By using the Cauchy-Schwarz inequality and well-known properties of
conditional PDFs we can bound the integral in the second term and the
integrals in the last summation as141414The inequality in (112) follows from
the Cauchy-Schwarz inequality. Specifically, let
$f(\bm{y})=\int_{\mathscr{R}(\bm{X}_{L-m})}\rho_{L-m}(\bm{y}-\bm{F}(\bm{z},\bm{w}_{L-m}))d\bm{z}.$
(110) Then
$\underbrace{\int_{\mathscr{R}(\bm{X}_{L-m+1})}1^{2}d\bm{y}}_{\lambda(\mathscr{R}(\bm{X}_{L-m+1}))}\int_{\mathscr{R}(\bm{X}_{L-m+1})}f(\bm{y})^{2}d\bm{y}\geq\left(\int_{\mathscr{R}(\bm{X}_{L-m+1})}1\cdot
f(\bm{y})d\bm{y}\right)^{2}=\lambda(\mathscr{R}(\bm{X}_{L-m}))^{2}.$ (111)

$\displaystyle\frac{\lambda(\mathscr{R}(\bm{X}_{L-m}))^{2}}{\lambda(\mathscr{R}(\bm{X}_{L-m+1}))}\leq\int_{\mathscr{R}(\bm{X}_{L-m+1})}\left(\int_{\mathscr{R}(\bm{X}_{L-m})}\rho_{L-m}(\bm{y}-\bm{F}(\bm{z},\bm{w}_{L-m}))d\bm{z}\right)^{2}d\bm{y},$
(112)

and

$\displaystyle\sum_{k=1}^{M}\int_{\mathscr{R}(\bm{X}_{L-m+1})}\left(\int_{\mathscr{R}(\bm{X}_{L-m})}\rho_{L-m}(\bm{y}-\bm{F}(\bm{z},\bm{w}_{L-m}))\eta^{m}_{k}(\bm{z})d\bm{z}\right)^{2}d\bm{y}\geq
0.$ (113)

By combining (109)-(113) we obtain

$\displaystyle\left\|\mathcal{Q}_{m}\mathcal{G}(m,m-1)\right\|^{2}\leq$
$\displaystyle\left\|\gamma_{L-m}\right\|^{2}_{L^{2}(\mathscr{R}(\bm{X}_{L-m+1})\times\mathscr{R}(\bm{X}_{L-m}))}$
$\displaystyle\leq$
$\displaystyle\lambda(\mathscr{R}(\bm{X}_{L-m}))\left\|\rho_{L-m}\right\|_{L^{2}(\mathscr{R}(\bm{\xi}_{L-m}))}^{2}-\frac{\lambda(\mathscr{R}(\bm{X}_{L-m}))}{\lambda(\mathscr{R}(\bm{X}_{L-m+1}))}$
(114) $\displaystyle\leq$
$\displaystyle\lambda(\Omega_{L-m})\left(\left\|\rho_{L-m}\right\|_{L^{2}(\mathscr{R}(\bm{\xi}_{L-m}))}^{2}-\frac{1}{\lambda(\Omega_{L-m+1})}\right),$
(115)

which proves the Theorem.

∎

Remark: The last two terms in (109) represent the $L^{2}$ norm of the
projection of $\rho_{L-m}$ onto the orthonormal basis
$\\{\lambda(\mathscr{R}(\bm{X}_{L-m}))^{-1/2},\eta^{m}_{1},\ldots,\eta^{m}_{M}\\}$.
As a matter of fact, if we assume that
$\rho_{L-m}(\bm{y}-\bm{F}(\bm{x},\bm{w}_{L-m})$ is in
$L^{2}(\mathscr{R}(\bm{X}_{L-m+1})\times\mathscr{R}(\bm{X}_{L-m}))$, then by
using Parseval’s identity we can write (105) as

$\displaystyle\int_{\mathscr{R}(\bm{X}_{L-m+1})}\int_{\mathscr{R}(\bm{X}_{L-m})}\gamma_{L-m}(\bm{y},\bm{x})^{2}d\bm{y}d\bm{x}=$
$\displaystyle\sum_{k=M+1}^{\infty}\int_{\mathscr{R}(\bm{X}_{L-m+1})}\left(\int_{\mathscr{R}(\bm{X}_{L-m})}\rho_{L-m}(\bm{y}-\bm{F}(\bm{z},\bm{w}_{L-m}))\eta^{m}_{k}(\bm{z})d\bm{z}\right)^{2}d\bm{y},$
(116)

where $\\{\eta_{M+1},\eta_{M+2}\ldots\\}$ is an orthonormal basis for the
orthogonal complement (in $L^{2}(\mathscr{R}(\bm{X}_{L-m}))$) of the space
spanned by the basis
$\\{\lambda(\mathscr{R}(\bm{X}_{L-m}))^{-1/2},\eta^{m}_{1},\ldots,\eta^{m}_{M}\\}$.
This allows us to bound (113) from below (with a nonzero bound). Such lower
bound depends on the basis functions $\eta^{m}_{k}$, on the weights
$\bm{w}_{L-m}$ as well as on the choice of the activation function $\bm{F}$.
In other words, the bound (103) can be improved, if we provide information on
$\eta^{m}_{k}$ and the activation function $\bm{F}$. Note also that the bound
(103) is formulated in terms of the Lebesgue measure of $\Omega_{L-m}$, i.e.,
$\lambda(\Omega_{L-m})$, rather than
$\lambda\left(\mathscr{R}(\bm{X}_{L-m})\right)$ as in equation (114). The
reason is that $\lambda(\Omega_{L-m})$ depends only on the range of the noise
(see definition (140) and Lemma 8), while
$\lambda\left(\mathscr{R}(\bm{X}_{L-m})\right)$ depends on the range of the
noise, on the weights of the layer $L-m$, and on the range of
$\bm{X}_{L-m+1}$.

######  Lemma 4.

Consider the projection kernel (96) with
$c_{m}=1/\lambda(\mathscr{R}(\bm{X}_{L-m}))$ and let $\kappa\geq 0$. If

$\left\|\rho_{L-m}\right\|_{L^{2}(\mathscr{R}(\bm{\xi}_{L-m}))}^{2}\leq\frac{\kappa}{\lambda(\Omega_{L-m})}+\frac{1}{\lambda(\Omega_{L-m+1})},$
(117)

then

$\left\|\mathcal{Q}_{m}\mathcal{G}(m,m-1)\right\|^{2}\leq\kappa.$ (118)

In particular, if $0\leq\kappa<1$ then $\mathcal{Q}_{m}\mathcal{G}(m,m-1)$ is
a contraction.

###### Proof.

The proof follows from equation (103).

∎

The upper bound in (117) represents a slight improvement of the bound we
obtained in A, Lemma 19. It is worthwhile noting that by using (167) we can
write (103) in terms of the max and min of the PDF $\rho_{L-m}$ in the domain
$\mathscr{R}(\bm{\xi}_{L-m})$.

###  5.2 Contractions induced by uniform random noise

Let $\\{\bm{\xi}_{0},\bm{\xi}_{1},,\ldots\\}$ be a discrete Markov process
driving the neural network model (6). As a concrete example to examine the
property in Lemma 4, suppose that each $\bm{\xi}_{i}$ is a random vector with
i.i.d. uniform components supported in $[-b_{i},b_{i}]$ ($b_{i}>0$). In this
assumption, the $L^{2}(\mathscr{R}(\bm{\xi}_{L-m}))$ norm of $\rho_{L-m}$
appearing in Theorem 3 and Lemma 4 can be computed analytically as

$\left\|\rho_{L-m}\right\|_{L^{2}(\mathscr{R}(\bm{\xi}_{L-m}))}^{2}=\frac{1}{\lambda(\mathscr{R}(\bm{\xi}_{L-m}))}=\frac{1}{(2b_{L-m})^{N}}.$
(119)

where $N$ is the number of neurons in each layer. By using Lemma 11 in A we
also have that $\lambda(\Omega_{L-m})=2^{N}(1+b_{L-m-1})^{N}.$ This allows us
to write the inequality (117) as

$\left(\frac{1+b_{L-m-1}}{b_{L-m}}\right)^{N}\leq\kappa+\left(\frac{1+b_{L-m-1}}{1+b_{L-m}}\right)^{N}.$
(120)

By setting $n=L-m$, we obtain

$\left(\frac{1+b_{n-1}}{b_{n}}\right)^{N}\leq\kappa+\left(\frac{1+b_{n-1}}{1+b_{n}}\right)^{N},$
(121)

i.e.,

$\frac{b_{n}(b_{n}+1)}{\left[(b_{n}+1)^{N}-b_{n}^{N}\right]^{1/N}}\geq\frac{b_{n-1}+1}{\kappa^{1/N}}\qquad
n=1,\ldots,L-1.$ (122)

A lower bound for the coefficient $b_{0}$ can be determined by using
Proposition 20 in A, i.e.,

$b_{0}\geq\frac{1}{2}\left(\frac{\lambda(\Omega_{0})}{\kappa}\right)^{1/N}.$
(123)

With a lower bound for $b_{0}$ available, we can compute a lower bound for
each $b_{n}$ ($n=1,2,\ldots$) by solving the recursion (122) with an equality
sign.

$\kappa=0.2$ $\kappa=0.2$

$\kappa=10^{-4}$ $\kappa=10^{-4}$

Figure 4: Lower bound on the noise amplitude (122) versus the number of
neurons $N$ for $\lambda(\Omega_{0})=1$ and different contraction factors
$\kappa$. With these values of $b_{n}$ the operator
$\mathcal{Q}_{L-n+1}\mathcal{G}(L-n+1,L-n)$ is a contraction satisfying
$\left\|\mathcal{Q}_{L-n+1}\mathcal{G}(L-n+1,L-n)\right\|^{2}\leq\kappa$
regardless of the weights in the neural network. A comparison between these
plots and the plots in Figure 8 is provided in Figure 5.

It is interesting to compare the sequence of lower bounds on $b_{n}$ generated
by the recursion (122) with the one generated by (188) (see also Proposition
2). This is done in Figure 4, where we plot the lower bound on $b_{n}$ we
obtained by solving the recursion (122) numerically for different contraction
factors $\kappa$151515To compute the lower bounds of $b_{n}$ we solved the
recursion (122) numerically (with an equality sign) for $b_{n}$ using the
Newton method. To improve numerical accuracy we wrote the left hand side of
(122) in the equivalent form
$\frac{b_{n}(b_{n}+1)}{\left[(b_{n}+1)^{N}-b_{n}^{N}\right]^{1/N}}=\displaystyle\frac{b_{n}}{\left[1-\left(\displaystyle\frac{1}{b_{n}}+1\right)^{-N}\right]^{1/N}}.$
(124) . The plots in Figure 4 and Figure 8 show that dependence of $b_{n}$ on
the number of neurons is indeed different if we consider the operator
$\mathcal{G}(L-n+1,L-n)$ alone or its projected version
$\mathcal{Q}_{L-n+1}\mathcal{G}(L-n+1,L-n)$. In particular, in Figure 4 we
notice the existence of a hump that is not present in Figure 8. Such a hump
allows for smaller values of $b_{n}$ relative to the case where we consider
the operator $\mathcal{G}$ alone. Moreover, the asymptotic value of $b_{n}$
for $N\rightarrow\infty$ (infinite number of neurons) is the same in the plots
in Figure 8 and Figure 4. This is demonstrated numerically in Figure 5 and
theoretically in the following Lemma.

$\kappa=0.2$ $\kappa=10^{-4}$

Figure 5: Comparison between the sequences of lower bounds on $b_{n}$
generated by the recursion (122) (operator $\mathcal{Q}\mathcal{G}$, black
lines) and the recursion (188) (operator $\mathcal{G}$, blue lines) for
$\lambda(\Omega_{0})=1$ and different contraction factors $\kappa$.

######  Lemma 5.

The asymptotic result of Lemma 21 holds also for the sequence of lower bounds
generated by (122), for any contraction factor $\kappa\geq 0$ and any
$\Omega_{0}$ (domain of the neural network input). In other words, the lower
bounds generated by (122) satisfy the condition

$\lim_{N\rightarrow\infty}b_{n}=\frac{1}{2}+n\qquad n=0,\ldots,L-1$ (125)

independently of the contraction factor $\kappa$ and $\Omega_{0}$.

###### Proof.

The proof follows by substituting the identity

$\lim_{N\rightarrow\infty}\left[(b_{n}+1)^{N}-b_{n}^{N}\right]^{1/N}=b_{n}+1$
(126)

into (122). This yields (183), which is the starting point of Lemma 21.

∎

Remark: In the region between the back and blue curvesin Figure 5, we have
that the operator $\mathcal{Q}_{L-n+1}\mathcal{G}(L-n+1,L-n)$ is certainly a
contraction while $\mathcal{G}(L-n+1,L-n)$ maybe a contraction or an
expansion. On the other hand, if we fix the number of neurons $N$ in each
layer, say $N=50$, then the uniform noise amplitude $b_{n}$ that induces a
certain contraction factor $\kappa$ in $\mathcal{G}(m+1,m)$ yields a much
larger contraction factor in $\mathcal{Q}_{m+1}\mathcal{G}(m+1,n)$. For
instance, if we choose $b_{n}$ as in equation (188) then we know that
$\mathcal{G}(L-n+1,L-n)$ is contractive, with contraction factor $\kappa$.
Corresponding to such sequence of $b_{n}$ we have that
$\mathcal{Q}_{L-n+1}\mathcal{G}(L-n+1,L-n)$ has contraction factors

$\kappa_{0}\geq\frac{\lambda(\Omega_{0})}{2^{N}}\left[\frac{1}{b^{N}_{0}}-\frac{1}{(b_{0}+1)^{N}}\right],\qquad\kappa_{n}\geq\left(\frac{1+b_{n-1}}{b_{n}}\right)^{N}-\left(\frac{1+b_{n-1}}{1+b_{n}}\right)^{N}.$
(127)

The lower bound of such contraction factors (i.e., the right hand side of
(127)) is plotted in Figure 6. In summary, if $\mathcal{G}(L-n+1,L-n)$ is a
contraction with contraction factor $\kappa<1$ across the layers of the neural
network, then the operator $\mathcal{Q}_{L-n+1}\mathcal{G}(L-n+1,L-n)$
contracts even more, i.e., it has a contraction factor $\kappa_{n}<\kappa$.

$\kappa=0.2$ $\kappa=10^{-4}$

Figure 6: Lower bound on the contraction factors $\kappa_{n}$ defined in (127)
versus the number of neurons. Here we choose $b_{n}$ as in (188), i.e., in a
way that guarantees that all $\mathcal{G}(L-n+1,L-n)$ are contractions with
the same constant contraction factor $\kappa$. With the factors $\kappa_{n}$
shown above we have
$\|\mathcal{Q}_{L-n+1}\mathcal{G}(L-n+1,L-n)\|^{2}\leq\kappa_{n}$. Hence, if
$\mathcal{G}(L-n+1,L-n)$ is a contraction with constant contraction factor
$\kappa<1$ then $\mathcal{Q}_{L-n+1}\mathcal{G}(L-n+1,L-n)$ contracts even
more, i.e., it has a contraction factor $\kappa_{n}<\kappa$.

###  5.3 Decay property of the neural network memory operator

Now we discuss the implications of the contraction property of
$\mathcal{Q}_{m}\mathcal{G}(m,m-1)$ on the MZ equations. It is straightforward
to show that if Proposition 2 or Lemma 4 hold true then the MZ memory and
noise terms in (82) decay with the number of layers. The key property is
summarized in the following Theorem.

######  Theorem 6.

If the conditions of Lemma 4 are satisfied, then the MZ memory operator (83)
decays with the number of layers in the neural network, i.e.,

$\left\|\Phi_{\mathcal{G}}(n,m+1)\right\|^{2}\geq\frac{1}{\kappa}\left\|\Phi_{\mathcal{G}}(n,m)\right\|^{2}\qquad\forall
n\geq m+1,\qquad 0<\kappa<1.$ (128)

Moreover,

$\left\|\Phi_{\mathcal{G}}(n,0)\right\|^{2}\leq\kappa^{n},$ (129)

i.e., the memory operator $\Phi_{\mathcal{G}}(n,0)$ decays exponentially fast
with the number of layers.

###### Proof.

The proof follows from
$\left\|Q_{m+1}\mathcal{G}(m+1,m)\right\|^{2}\leq\kappa$ and equa (83). In
fact, for all $n\geq m+1$

$\left\|\Phi_{\mathcal{G}}(n,m)\right\|^{2}=\left\|\Phi_{\mathcal{G}}(n,m+1)\mathcal{Q}_{m+1}\mathcal{G}(m+1,m)\right\|^{2}\leq\kappa\left\|\Phi_{\mathcal{G}}(n,m+1)\right\|^{2}.$
(130)

∎

This result can be used to approximate the MZ equation of a neural net with a
large number of layers to an equivalent one involving only a few layers. Note
that the fading memory property (129) holds for the conditional expectation of
any phase space function, including conditional mean and conditional second-
order moment of the neural network vector of state.

###  5.4 Reducing deep neural networks to shallow neural networks

Consider the MZ equation (85), hereafter rewritten for convenience

$\bm{q}_{L}=\Phi_{\mathcal{G}}(L,0)\bm{q}_{0}+\sum_{m=0}^{L-1}\mathcal{G}(L,m+1)\mathcal{P}_{m+1}\mathcal{G}(m+1,m)\Phi_{\mathcal{G}}(m,0)\bm{q}_{0}.$
(131)

We have seen in Lemma 4 that the memory term $\Phi_{\mathcal{G}}(L,0)$ decays
exponentially fast with the number of layers if the operator
$\mathcal{Q}_{m}\mathcal{G}(m,m-1)$ is a contraction. Specifically, we proved
that

$\left\|\Phi_{\mathcal{G}}(L,0)\right\|^{2}\leq\kappa^{L}\qquad
0\leq\kappa<1,$ (132)

where $\kappa$ is a parameter that can be made as small as we like. If $L$ is
sufficiently large and $\kappa$ is sufficiently small we have that the first
term at the right hand side of (131) is negligible. Regarding the other terms
in the summation, we show hereafter that they can be controlled with $\kappa$
independently of the neural network weights, and that they decay with the
layer index $m$.

######  Proposition 7.

Under the same assumptions of Theorem 14 in A we have that

$\left\|\mathcal{G}(L,m+1)\mathcal{P}_{m+1}\mathcal{G}(m+1,m)\Phi_{\mathcal{G}}(m,0)\right\|^{2}\leq
B^{L}\left(\frac{\kappa}{B}\right)^{m},$ (133)

where $\kappa$ is defined in Lemma 3 and

$B\geq\sup_{i=1,\ldots,L-1}\left\|\mathcal{G}(i+1,i)\right\|^{2}\qquad(B<\infty)$
(134)

is independent of the neural network weights161616Note that if  $0<\kappa<B$
(135) then the upper bound in (133) decays with $m$. Also, if $B>1$ and
$k=B^{1-L}$ then the upper bound in (133) is exactly $B^{1-mL}$..

###### Proof.

By using standard operator norm inequalities and Theorem 6 we immediately
obtain

$\displaystyle\left\|\mathcal{G}(L,m+1)\mathcal{P}\mathcal{G}(m+1,m)\Phi_{\mathcal{G}}(m,0)\right\|^{2}\leq$
$\displaystyle\left\|\mathcal{G}(L,L-1)\cdots\mathcal{G}(m+2,m+1)\right\|^{2}\left\|\mathcal{G}(m+1,m)\right\|^{2}\kappa^{m}$
$\displaystyle\leq$
$\displaystyle\left(\sup_{i=m,\ldots,L-1}\left\|\mathcal{G}(i+1,i)\right\|^{2}\right)^{L-m}\kappa^{m}$
$\displaystyle\leq$ $\displaystyle B^{L-m}\kappa^{m},$ (136)

where $B$ is defined in (134). In Proposition 17 we proved that $B$ exists, is
finite, and it does not depend on the weights of the neural network.

∎

##  6 Summary

We developed a new formulation of deep learning based on the Mori-Zwanzig (MZ)
projection operator formalism of irreversible statistical mechanics. The new
formulation provides new insights on how information propagates through deep
neural networks in terms formally exact linear operator equations, and it
introduces a new important concept, i.e., the memory of the neural network,
which plays a fundamental role in low-dimensional modeling and
parameterization of the network (see, e.g., [25]). By using the theory of
contraction mappings, we developed sufficient conditions for the memory of the
neural network to decay with the number of layers. This allowed us to
rigorously transform deep networks into shallow ones, e.g., by reducing the
number of neurons per layer (using projections), or the by reducing the total
number of layers (using the decay property of the memory operator). We
developed most of our analysis for MZ equations involving conditional
expectations, i.e., Eqs. (82)-(85). However, by using the well-known duality
between PDF dynamics and conditional expectation dynamics [9], it is
straightforward to derive similar analytic results for MZ equations involving
PDFs, i.e., Eqs. (77)-(78). The mathematical techniques we developed in this
paper can be generalized to other types of stochastic neural network models.

An important open question is the development of effective approximation
methods for the MZ memory operator and the noise term. With a specific
reference to equation (66), it can be shown that the noise term has zero mean
if the projection operator is chosen appropriately. It would then be of
practical interest to study the fluctuation and concentration properties of
such zero-mean noise term with respect to the size of the neural network.
Other approximation methods for the MZ memory and noise terms can be built on
the backbone of continuous-time approximation methods, e.g., based on
functional analysis [53, 54, 29], combinatorics [55], data-driven methods [3,
40, 30, 33], Markovian embedding techniques [21, 18, 32, 25, 7], or
projections based on reproducing kernel Hilbert or Banach spaces [1, 39, 52].

## Acknowledgements

Dr. Venturi was partially supported by the U.S. Air Force Office of Scientific
Research grant FA9550-20-1-0174, and by the U.S. Army Research Office grant
W911NF1810309. Dr. Li was supported by the NSF grant DMS-1953120.

## Data availability statement

The data that support the findings of this study are available from the
corresponding author upon request.

##  Appendix A Functional setting

Let $(\mathcal{S},\mathcal{F},\mathscr{P})$ be a probability space. Consider
the neural network model (see Figure 2)

$\bm{X}_{1}=\bm{F}_{0}(\bm{X}_{0},\bm{w}_{0})+\bm{\xi}_{0}\qquad\bm{X}_{n+1}=\bm{F}(\bm{X}_{n},\bm{w}_{n})+\bm{\xi}_{n}\qquad
n=1,\ldots,L-1,$ (137)

where $\\{\bm{\xi}_{0},\ldots,\bm{\xi}_{L-1}\\}$ is a discrete (vector-valued)
Markov process. Suppose we are interested in using the model (137) to
approximate a multivariate function $f(\bm{x})$. This is usually done by
taking a linear combination of the network output, e.g., Eq. (7). In this
setting, the neural network can be thought of as a process of constructing an
adaptive basis by function composition. Here we consider the case where the
function we are approximating is defined on a compact subset $\Omega_{0}$ of
$\mathbb{R}^{d}$. This means that the input vector of the neural network, i.e.
$\bm{X}_{0}$, is and element of $\Omega_{0}$. We assume that the following
conditions are satisfied

  1. 1.

$\bm{X}_{0}\in\Omega_{0}\subseteq\mathbb{R}^{d}$ ($\Omega$ compact),
$\bm{X}_{n}\in\mathbb{R}^{N}$ for $n=1,\ldots,L-1$;

  2. 2.

The image of $\bm{F}_{0}$ and $\bm{F}$ is the hyper-cube $[-1,1]^{N}$.

For example, if $\bm{F}$ in (137) is of the form

$\bm{F}(\bm{x},\bm{w})=\tanh(\bm{W}\bm{x}+\bm{b})\qquad\bm{w}=\\{\bm{W},\bm{b}\\},$
(138)

then conditions 1. and 2. imply that $\bm{W}_{0}\in M_{N\times d}(\mathbb{R})$
and $\bm{W}_{n}\in M_{N\times N}(\mathbb{R})$ for $n=1,\ldots L-1$, while the
biases are $\bm{b}_{n}\in M_{N\times 1}(\mathbb{R})$ for $n=0,\ldots L-1$. The
random vectors $\\{\bm{\xi}_{0},\ldots,\bm{\xi}_{L-1}\\}$ added to the output
of each layer make $\\{\bm{X}_{1},\ldots,\bm{X}_{K}\\}$ a discrete Markov
process (each $\bm{X}_{i}$ is a random vector). The range of $\bm{X}_{n+1}$
depends essentially on the range of $\bm{\xi}_{n}$, as the image of $\bm{F}$
is the hyper-cube $[-1,1]^{N}$ (see condition 2. above). Let us
define171717The notation $[-1,1]^{N}$ denotes a Cartesian product of $N$ one-
dimensional domains $[-1,1]$, i.e.,
$[-1,1]^{N}=\bigtimes_{k=1}^{N}[-1,1]=\underbrace{[-1,1]\times[-1,1]\times\cdots\times[-1,1]}_{\text{$N$
times}}.$ (139)

$\displaystyle\Omega_{n+1}=$
$\displaystyle[-1,1]^{N}+\mathscr{R}(\bm{\xi}_{n})$ $\displaystyle=$
$\displaystyle\\{\bm{c}\in\mathbb{R}^{N}\,:\,\bm{c}=\bm{a}+\bm{b}\quad\bm{a}\in[-1,1]^{N},\,\bm{b}\in\mathscr{R}(\bm{\xi}_{n})\\},$
(140)

where $\mathscr{R}(\bm{\xi}_{n})$ denotes the range of the random vector
$\bm{\xi}_{n}$, i.e.,

$\mathscr{R}(\bm{\xi}_{n})=\\{\bm{\xi}_{n}(\omega)\in\mathbb{R}^{N}:\omega\in\mathcal{S}\\}.$
(141)

Clearly, the range of the random vector $\bm{X}_{n+1}$ is a subset181818We
emphasize that if we are given a specific form of the activation function
$\bm{F}$ together with suitable bounds on the neural network weights and
biases $\\{\bm{W},\bm{b}\\}$ then we can easily identify a domain that is
smaller than $\Omega_{n}$, and that still contains $\mathscr{R}(\bm{X}_{n})$.
This allows us to construct a tighter bound for
$\lambda(\mathscr{R}(\bm{X}_{n+1})$ in Lemma 8, which depends on the
activation function and on the bounds we set on neural network weights and
biases. of $\Omega_{n+1}$, i.e.,
$\mathscr{R}(\bm{X}_{n+1})\subseteq\Omega_{n+1}$. This implies the following
lemma.

######  Lemma 8.

Let $\lambda(\Omega_{n+1})$ the Lebesgue measure of the set (140). The
Lebesgue measure of the range of $\bm{X}_{n+1}$ satisfies

$\lambda(\mathscr{R}(\bm{X}_{n+1}))\leq\lambda(\Omega_{n+1}).$ (142)

###### Proof.

The proof follows from the inclusion
$\mathscr{R}(\bm{X}_{n+1})\subseteq\Omega_{n+1}$. ∎

The $L^{\infty}$ norm of the random vector $\bm{\xi}$ is defined as the
largest value of $r\geq 0$ that yields a nonzero probability on the event
$\\{\omega\in\mathcal{S}:\,\left\|\bm{\xi}(\omega)\right\|_{\infty}>r\\}\in\mathcal{F}$,
i.e.,

$\left\|\bm{\xi}\right\|_{\infty}=\sup_{r\in\mathbb{R}}\\{\,\mathscr{P}(\\{\omega\in\mathcal{S}:\,\left\|\bm{\xi}(\omega)\right\|_{\infty}>r\\})>0\\}.$
(143)

This definition allows us to bound the Lebesgue measure of $\Omega_{n+1}$ as
follows.

######  Proposition 9.

The Lebesgue measure of the set $\Omega_{n+1}$ defined in (140) can be bounded
as

$\lambda(\Omega_{n+1})\leq\left(\sqrt{N}+\left\|\bm{\xi}_{n}\right\|_{\infty}\right)^{N}\frac{\pi^{N/2}}{\Gamma(1+N/2)},$
(144)

where $N$ is the number of neurons and $\Gamma(\cdot)$ is the Gamma function.

###### Proof.

As is well known, the length of the diagonal of the hypercube $[-1,1]^{N}$ is
$\sqrt{N}$. Hence, $\sqrt{N}+\left\|\bm{\xi}_{n}\right\|_{\infty}$ is the
radius of a ball that encloses all elements of $\Omega_{n+1}$. The Lebesgue
measure of such ball is obtained by multiplying the Lebesgue measure of the
unit ball in $\mathbb{R}^{N}$, i.e., $\pi^{N/2}/\Gamma(1+N/2)$ by the scaling
factor $\left(\sqrt{N}+\left\|\bm{\xi}_{n}\right\|_{\infty}\right)^{N}$.

∎

######  Lemma 10.

If $\mathscr{R}(\bm{\xi}_{n})$ is bounded then $\mathscr{R}(\bm{X}_{n+1})$ is
bounded.

###### Proof.

The image of the activation function $\bm{F}$ is a bounded set. If
$\mathscr{R}(\bm{\xi}_{n})$ is bounded then $\Omega_{n+1}$ in (140) is
bounded. $\mathscr{R}(\bm{X}_{n+1})\subseteq\Omega_{n+1}$ and therefore
$\mathscr{R}(\bm{X}_{n+1})$ is bounded.

∎

Clearly, if $\\{\bm{\xi}_{0},\ldots,\bm{\xi}_{L-1}\\}$ are i.i.d. random
variables then there exists a domain $V=\Omega_{1}=\cdots=\Omega_{L}$ such
that

$\mathscr{R}(\bm{\xi}_{n})\subseteq\mathscr{R}(\bm{X}_{n+1})\subseteq
V\qquad\forall n=0,\ldots,L-1.$ (145)

In fact, if $\\{\bm{\xi}_{0},\ldots,\bm{\xi}_{L-1}\\}$ are i.i.d. random
variables then we have

$\mathscr{R}(\bm{\xi}_{0})=\mathscr{R}(\bm{\xi}_{1})=\cdots=\mathscr{R}(\bm{\xi}_{L-1}),$
(146)

which implies that all of $\Omega_{i}$ defined in (140) are the same. If the
range of each random vector $\bm{\xi}_{n}$ is a tensor product of one-
dimensional domain, e.g., if the components of $\bm{\xi}_{n}$ are
statistically independent, then $V=\Omega_{1}=\cdots=\Omega_{L}$ becomes
particularly simple, i.e., a hypercube.

######  Lemma 11.

Let $\\{\bm{\xi}_{0},\ldots,\bm{\xi}_{L-1}\\}$ be i.i.d. random variables with
bounded range and suppose that each $\bm{\xi}_{k}$ has statistically
independent components with range $[a,b]$. Then all domains
$\\{\Omega_{1},\ldots,\Omega_{L}\\}$ defined in equation (140) are the same,
and they are equivalent to

$V=\bigtimes_{k=1}^{N}[-1+a,1+b].$ (147)

$V$ includes the range of all random vectors $\bm{X}_{n}$ ($n=1,\ldots,L$) and
has Lebesgue measure

$\lambda(V)=(2+b-a)^{N}.$ (148)

###### Proof.

The proof is trivial and therefore omitted.

∎

Remark: It is worth noticing that if each $\bm{\xi}_{k}$ is a uniformly
distributed random vector with statistically independent components in
$[-1,1]$, then for $N=10$ (number of neurons) the upper bound in (144) is
$3.98\times 10^{6}$ while the exact result (148) gives $1.05\times 10^{6}$.
Hence the estimate (144) is quite sharp in the case of uniform random vectors.

#### Boundedness of composition and transfer operators

Lemma 10 states that if we perturb the output of the $n$-th layer of a neural
network by a random vector $\bm{\xi}_{n}$ with finite range then we obtain a
random vector $\bm{X}_{n+1}$ with finite range. In this hypothesis it is
straightforward to show that that the composition and transfer operators
defined in (22) and (16) are bounded. We have seen that these operators can be
written as

$\mathcal{M}(n,n+1)v=\int_{\mathscr{R}(\bm{X}_{n+1})}v(\bm{y})p_{n+1|n}(\bm{y}|\bm{x})d\bm{y},\qquad\mathcal{N}(n+1,n)v=\int_{\mathscr{R}(\bm{X}_{n})}p_{n+1|n}(\bm{x}|\bm{y})v(\bm{y})d\bm{y},$
(149)

where $p_{n+1|n}(\bm{y}|\bm{x})=\rho_{n}(\bm{y}-\bm{F}(\bm{x},\bm{w}_{n}))$ is
the conditional transition density of $\bm{X}_{n+1}$ given $\bm{X}_{n}$, and
$\rho_{n}$ is the joint PDF of the random vector $\bm{\xi}_{n}$. The
conditional transition density $p_{n+1|n}(\bm{y}|\bm{x})$ is always non-
negative, i.e.,

$p_{n+1|n}(\bm{y}|\bm{x})\geq
0\qquad\forall\bm{y}\in\mathscr{R}(\bm{X}_{n+1}),\quad\forall\bm{x}\in\mathscr{R}(\bm{X}_{n}).$
(150)

Moreover, the conditional density $p_{n+1|n}$ is defined on the set

$\mathscr{B}_{n}=\\{(\bm{x},\bm{y})\in\mathscr{R}(\bm{X}_{n})\times\mathscr{R}(\bm{X}_{n+1}):\,(\bm{y}-\bm{F}(\bm{x},\bm{w}_{n}))\in\mathscr{R}(\bm{\xi}_{n})\\}.$
(151)

It is also important to emphasize that $\bm{y}\in\mathscr{R}(\bm{X}_{n+1})$
and $\bm{x}\in\mathscr{R}(\bm{X}_{n})$. Both $\mathscr{R}(\bm{X}_{n+1})$ and
$\mathscr{R}(\bm{X}_{n})$ depend on $\Omega_{0}$ (domain of the input), the
neural network weights, and the noise amplitude. Thanks to Lemma 8, we have
that

$\mathscr{B}_{n}\subseteq\Omega_{n}\times\Omega_{n+1}.$ (152)

The Lebesgue measure of $\mathscr{B}_{n}$ can be calculated as follows.

######  Lemma 12.

The Lebesgue measure of the set $\mathscr{B}_{n}$ defined in (151) is equal to
the product of the measure of $\lambda(\mathscr{R}(\bm{X}_{n}))$ and the
measure of $\mathscr{R}(\bm{\xi}_{n})$, i.e.,

$\lambda(\mathscr{B}_{n})=\lambda(\mathscr{R}(\bm{X}_{n}))\lambda(\mathscr{R}(\bm{\xi}_{n})).$
(153)

Moreover, $\lambda(\mathscr{B}_{n})$ is bounded by
$\lambda(\mathscr{R}(\Omega_{n}))\lambda(\mathscr{R}(\bm{\xi}_{n}))$, which is
independent of the neural network weights.

###### Proof.

Let $\chi_{n}$ be the indicator function of the set
$\mathscr{R}(\bm{\xi}_{n})$, $\bm{y}\in\mathscr{R}(\bm{X}_{n+1})$ and
$\bm{x}\in\mathscr{R}(\bm{X}_{n})$. Then

$\displaystyle\lambda(\mathscr{B}_{n})=$
$\displaystyle\int_{\mathscr{R}(\bm{X}_{n+1})}\int_{\mathscr{R}(\bm{X}_{n})}\chi_{n}(\bm{y}-\bm{F}(\bm{x},\bm{w}_{n}))d\bm{x}d\bm{y}$
$\displaystyle=$
$\displaystyle\lambda(\mathscr{R}(\bm{\xi}_{n}))\int_{\mathscr{R}(\bm{X}_{n})}d\bm{x}$
$\displaystyle=$
$\displaystyle\lambda(\mathscr{R}(\bm{X}_{n}))\lambda(\mathscr{R}(\bm{\xi}_{n})).$
(154)

By using Lemma 8 we conclude that $\lambda(\mathscr{B}_{n})$ is bounded from
above by
$\lambda(\mathscr{R}(\Omega_{L-m}))\lambda(\mathscr{R}(\bm{\xi}_{n}))$, which
is independent of the neural network weights.

∎

Remark: The result (153) has a straightforward geometrical interpretation in
two dimensions. Pick a ruler of length $r=\lambda(\mathscr{R}(\xi_{n}))$ with
endpoints that can leave markings if we slide it on a rectangular table with
side lengths $s_{b}=\lambda(\mathscr{R}(X_{n+1}))$ (horizontal side)
$s_{h}=\lambda(\mathscr{R}(X_{n}))$ (vertical side). Then slide the ruler from
the top to the bottom of the table, while keeping it horizontal, i.e.,
parallel to the horizontal sides of the table (see Figure 7). The area of the
domain defined by the two curves drawn by the endpoints of the ruler is always
$r\times s_{h}$ independently of the way we slide the ruler laterally –
provided the ruler never gets out of the table.

######  Lemma 13.

If the range of $\bm{\xi}_{n-1}$ is a bounded subset of $\mathbb{R}^{N}$ then
the transition density $p_{n+1|n}(\bm{y}|\bm{x})$ is an element of
$L^{1}\left(\mathscr{R}(\bm{X}_{n+1})\times\mathscr{R}(\bm{X}_{n})\right)$.

###### Proof.

Note that

$\int_{\mathscr{R}(\bm{X}_{n+1})}\int_{\mathscr{R}(\bm{X}_{n})}p_{n+1|n}(\bm{y}|\bm{x})d\bm{y}d\bm{x}=\lambda\left(\mathscr{R}(\bm{X}_{n})\right)\leq\lambda(\Omega_{n}).$
(155)

The Lebesgue measure $\lambda(\Omega_{n})$ can be bounded as (see Proposition
9)

$\lambda(\Omega_{n})\leq\left(\sqrt{N}+\left\|\bm{\xi}_{n-1}\right\|_{\infty}\right)^{N}\frac{\pi^{N/2}}{\Gamma(1+N/2)}.$
(156)

Since the range of $\bm{\xi}_{n-1}$ is bounded by hypothesis we have that
there exists a finite real number $M>0$ such that
$\left\|\bm{\xi}_{n-1}\right\|_{\infty}\leq M$. This implies that the integral
in (155) is finite, i.e., that the transition kernel
$p_{n+1|n}(\bm{y}|\bm{x})$ is in
$L^{1}\left(\mathscr{R}(\bm{X}_{n+1})\times\mathscr{R}(\bm{X}_{n})\right)$.

∎

######  Theorem 14.

Let $C_{\bm{\xi}_{n}}(\bm{x})$ be the cumulative distribution function
$\bm{\xi}_{n}$. If $C_{\bm{\xi}_{n}}(\bm{x})$ is Lipschitz continuous on
$\mathscr{R}(\bm{\xi}_{n})$ and the partial derivatives $\partial
C_{\bm{\xi}_{n}}/\partial x_{k}$ ($k=1,\ldots,N$) are Lipschitz continuous in
$x_{1}$, $x_{2}$, …, $x_{N}$, respectively, then the joint probability density
function of $\bm{\xi}_{n}$ is bounded on $\mathscr{R}(\bm{\xi}_{n})$.

###### Proof.

By using Rademacher’s theorem we have that if $C_{\bm{\xi}_{n}}(\bm{x})$ is
Lipschitz on $\mathscr{R}(\bm{\xi}_{n})$ then it is differentiable almost
everywhere on $\mathscr{R}(\bm{\xi}_{n})$ (except on a set with zero Lebesgue
measure). Therefore the partial derivatives $\partial
C_{\bm{\xi}_{n}}/\partial x_{k}$ exist almost everywhere on
$\mathscr{R}(\bm{\xi}_{n})$. If, in addition, we assume that $\partial
C_{\bm{\xi}_{n}}/\partial x_{k}$ are Lipschitz continuous with respect to
$x_{k}$ (for all $k=1,\ldots,N$) then by applying [34, Theorem 9] recursively
we conclude that the joint probability density function of $\bm{\xi}_{n}$ is
bounded.

∎

######  Lemma 15.

Under the same assumptions of Theorem 14 we have that the conditional PDF
$p_{n+1|n}(\bm{y}|\bm{x})=\rho_{n}(\bm{y}-\bm{F}(\bm{x},\bm{w}))$ is bounded
on $\mathscr{R}(\bm{X}_{n+1})\times\mathscr{R}(\bm{X}_{n})$.

###### Proof.

Theorem 14 states that $\rho_{n}$ is a bounded function. This implies that the
conditional density
$p_{n+1|n}(\bm{y}|\bm{x})=\rho_{n}(\bm{y}-\bm{F}(\bm{x},\bm{w}))$ is bounded
on $\mathscr{R}(\bm{X}_{n+1})\times\mathscr{R}(\bm{X}_{n})$.

∎

######  Proposition 16.

Let $\mathscr{R}(\bm{\xi}_{n})$ and $\mathscr{R}(\bm{\xi}_{n-1})$ be bounded
subsets of $\mathbb{R}^{N}$. Then, under the same assumptions of Theorem 14,
we have that the composition and the transfer operators defined in (149) are
bounded in $L^{2}$.

###### Proof.

Let us first prove that $\mathcal{M}(n,n+1)$ is a bounded linear operator from
$L^{2}(\mathscr{R}(\bm{X}_{n+1}))$ into $L^{2}(\mathscr{R}(\bm{X}_{n}))$. To
this end, note that

$\displaystyle\left\|\mathcal{M}(n,n+1)v\right\|^{2}_{L^{2}(\mathscr{R}(\bm{X}_{n}))}=$
$\displaystyle\int_{\mathscr{R}(\bm{X}_{n})}\left|\int_{\mathscr{R}(\bm{X}_{n+1})}v(\bm{y})p_{n+1|n}(\bm{y}|\bm{x})d\bm{y}\right|^{2}d\bm{x}$
$\displaystyle\leq$
$\displaystyle\left\|v\right\|^{2}_{L^{2}(\mathscr{R}(\bm{X}_{n+1}))}\underbrace{\int_{\mathscr{R}(\bm{X}_{n})}\int_{\mathscr{R}(\bm{X}_{n+1})}p_{n+1|n}(\bm{y}|\bm{x})^{2}d\bm{y}d\bm{x}}_{K_{n}}$
$\displaystyle=$ $\displaystyle
K_{n}\left\|v\right\|^{2}_{L^{2}(\mathscr{R}(\bm{X}_{n+1}))}.$ (157)

Clearly, $K_{n}<\infty$. In fact, if $\mathscr{R}(\bm{\xi}_{n})$ and
$\mathscr{R}(\bm{\xi}_{n-1})$ are bounded then $\mathscr{R}(\bm{X}_{n+1})$ and
$\mathscr{R}(\bm{X}_{n})$ are bounded. Moreover, thanks to Lemma 15 we have
that $p_{n+1|n}(\bm{y}|\bm{x})$ is bounded on
$\mathscr{R}(\bm{X}_{n+1})\times\mathscr{R}(\bm{X}_{n})$. Hence, $K_{n}$ is
the integral of the square of a bounded function defined on a bounded domain,
and therefore it is finite. By following the same steps it is straightforward
to show that the transfer operator $\mathcal{N}$ is a bounded linear operator.
Alternatively, simply recall that $\mathcal{N}$ is the adjoint of
$\mathcal{M}$, and the adjoint of a bounded linear operator is bounded.
Specifically we have,

$\left\|\mathcal{N}(n+1,n)p\right\|^{2}_{L^{2}(\mathscr{R}(\bm{X}_{n+1}))}\leq
K_{n}\left\|p\right\|^{2}_{L^{2}(\mathscr{R}(\bm{X}_{n}))}.$ (158)

∎

Remark: The integrals

$K_{n}=\int_{\mathscr{R}(\bm{X}_{n})}\int_{\mathscr{R}(\bm{X}_{n+1})}p_{n+1|n}(\bm{y}|\bm{x})^{2}d\bm{y}d\bm{x}$
(159)

can be computed by noting that

$p_{n+1|n}(\bm{y}|\bm{x})=\rho_{n}(\bm{y}-\bm{F}(\bm{x},\bm{w}))$ (160)

is essentially a shift of the PDF $\rho_{n}$ by a quantity
$\bm{F}(\bm{x},\bm{w})$ that depends on $\bm{x}$ and $\bm{w}$ (see, e.g.,
Figure 7). Such a shift does not influence the integral with respect to
$\bm{y}$, meaning that the integral of $p_{n+1|n}(\bm{y}|\bm{x})$ or
$p_{n+1|n}(\bm{y}|\bm{x})^{2}$ with respect to $\bm{y}$ is the same for all
$\bm{x}$. Hence, by changing variables we have that the integral (159) is
equivalent to

$K_{n}=\lambda(\mathscr{R}(\bm{X}_{n}))\int_{\mathscr{R}(\bm{\xi}_{n})}\rho_{n}(\bm{x})^{2}d\bm{x},$
(161)

where $\lambda(\mathscr{R}(\bm{X}_{n}))$ is the Lebesgue measure of
$\mathscr{R}(\bm{X}_{n})$, and $\mathscr{R}(\bm{\xi}_{n})$ is the range of
$\bm{\xi}_{n}$. Note that $K_{n}$ depends on the neural net weights only
through the Lebesgue measure of $\mathscr{R}(\bm{X}_{n})$. Clearly, since the
set $\Omega_{n}$ includes $\mathscr{R}(\bm{X}_{n})$ we have by Lemma 8 that
$\lambda(\mathscr{R}(\bm{X}_{n}))\leq\lambda(\Omega_{n})$. This implies that

$K_{n}\leq\lambda(\Omega_{n})\int_{\mathscr{R}(\bm{\xi}_{n})}\rho_{n}(\bm{x})^{2}d\bm{x}.$
(162)

The upper bound here does not depend on the neural network weights. The
following lemma summarizes all these remarks.

######  Proposition 17.

Under the same assumptions of Theorem 14, we have that the composition and the
transfer operators defined in (149) can be bounded as

$\left\|\mathcal{M}(n,n+1)\right\|^{2}\leq
K_{n},\qquad\left\|\mathcal{N}(n+1,n)\right\|^{2}\leq K_{n},$ (163)

where

$K_{n}=\lambda(\mathscr{R}(\bm{X}_{n}))\int_{\mathscr{R}(\bm{\xi}_{n})}\rho_{n}(\bm{x})^{2}d\bm{x}.$
(164)

Moreover, $K_{n}$ can be bounded as

$K_{n}\leq\lambda(\Omega_{n})\int_{\mathscr{R}(\bm{\xi}_{n})}\rho_{n}(\bm{x})^{2}d\bm{x},$
(165)

where $\Omega_{n}$ is defined in (140) and $\rho_{n}$ is the PDF of
$\bm{\xi}_{n}$. The upper bound in (165) does not depend on the neural network
weights and biases.

Under additional assumptions on the PDF $\rho_{n}(\bm{x})$ it is also possible
to bound the integrals at the right hand side of (164) and (165). Specifically
we have the following sharp bound.

######  Lemma 18.

Let

$s_{n}=\inf_{\bm{x}\in\mathscr{R}(\bm{\xi}_{n})}{\rho_{n}}(\bm{x}),\qquad
S_{n}=\sup_{\bm{x}\in\mathscr{R}(\bm{\xi}_{n})}{\rho_{n}}(\bm{x}).$ (166)

If $s_{n}>0$ then under the same assumptions of Theorem 14 we have that

$\int_{\mathscr{R}(\bm{\xi}_{n})}\rho_{n}(\bm{x})^{2}d\bm{x}\leq\frac{1}{\lambda(\mathscr{R}(\bm{\xi}_{n}))}\frac{(S_{n}+s_{n})^{2}}{4S_{n}s_{n}}.$
(167)

###### Proof.

First we notice that if the random vector $\bm{\xi}_{n}$ satisfies the
assumptions of Theorem 14 then the upper bound $S_{n}$ is finite. By using the
definition (166) we have

$(\rho_{n}(\bm{x})-s_{n})(S_{n}-\rho_{n}(\bm{x}))\geq 0\quad\text{for
all}\quad\bm{x}\in\mathscr{R}(\bm{\xi}_{n}).$ (168)

This implies

$\int_{\mathscr{R}(\bm{\xi}_{n})}\rho_{n}(\bm{x})^{2}d\bm{x}\leq(S_{n}+s_{n})-S_{n}s_{n}\lambda(\mathscr{R}(\bm{\xi}_{n})),$
(169)

where we used the fact that the PDF $\rho_{n}$ integrates to one over
$\mathscr{R}(\bm{\xi}_{n})$. Next, define

$R_{n}=\frac{1}{\lambda(\mathscr{R}(\bm{\xi}_{n}))}\frac{(S_{n}+s_{n})^{2}}{4S_{n}s_{n}}.$
(170)

Clearly,

$\displaystyle
R_{n}\left(1-\frac{2S_{n}s_{n}}{s_{n}+S_{n}}\lambda(\mathscr{R}(\bm{\xi}_{n}))\right)^{2}=R_{n}-(S_{n}+s_{n})+S_{n}s_{n}\lambda(\mathscr{R}(\bm{\xi}_{n}))\geq
0$ (171)

which implies that

$(S_{n}+s_{n})-S_{n}s_{n}\lambda(\mathscr{R}(\bm{\xi}_{n}))\leq R_{n}.$ (172)

A substitution of (172) into (169) yields (167).

∎

An example: Let $X_{0}\in\Omega_{0}=[-1,1]$ and consider

$X_{1}=\tanh(X_{0}+3)+\xi_{0},\qquad X_{2}=\tanh(2X_{1}-1)+\xi_{1},$ (173)

where $\xi_{0}$ and $\xi_{1}$ are uniform random variables with range
$\mathscr{R}(\xi_{0})=\mathscr{R}(\xi_{1})=[-2,2]$. In this setting,

$\displaystyle\mathscr{R}(X_{1})=$ $\displaystyle[\tanh(2)-2,\tanh(4)+2],$
$\displaystyle\mathscr{R}(X_{2})=$
$\displaystyle[\tanh(2\tanh(2)-5)-2,\tanh(2\tanh(4)+3)+2].$

The conditional density of $X_{1}$ given $X_{0}$ is given by

$p_{1|0}(x_{1}|x_{0})=\begin{cases}\displaystyle\frac{1}{4}&\text{if
$\left|x_{1}-\tanh(x_{0}+3)\right|\leq 2$}\\\ 0&\text{otherwise}\end{cases}$
(174)

This function is plotted in Figure 7 together with the domain
$\mathscr{R}(X_{1})\times\mathscr{R}(X_{0})$ (interior of the rectangle
delimited by dashed red lines).

Figure 7: Conditional probability density function $p_{1|0}(x_{1}|x_{0})$
defined in equation (174). The domain
$\mathscr{R}(X_{1})\times\mathscr{R}(X_{0})$ is the interior of the rectangle
delimited by dashed red lines.

Clearly, the integral of the conditional PDF (174) is

$\int_{\mathscr{R}(X_{0})}\int_{\mathscr{R}(X_{1})}p_{1|0}(x_{1}|x_{0})dx_{1}dx_{0}=\lambda(\mathscr{R}(X_{0}))=2,$
(175)

where $\lambda(\mathscr{R}(X_{0}))$ is the Lebesgue measure of
$\mathscr{R}(X_{0})=[-1,1]$. The $L^{2}$ norm of the operators $\mathcal{N}$
and $\mathcal{M}$ is bounded by191919For uniformly distributed random
variables we have that
$\int_{\mathscr{R}(\bm{\xi}_{n})}\rho_{n}(\bm{x})^{2}d\bm{x}=\frac{1}{\lambda(\mathscr{R}(\bm{\xi}_{n}))}.$
(176) Therefore equation (164) yields
$K_{n}=\frac{\lambda\left(\mathscr{R}(\bm{X}_{n})\right)}{\lambda\left(\mathscr{R}(\bm{\xi}_{n})\right)}\leq\frac{\lambda\left(\Omega_{n}\right)}{\lambda\left(\mathscr{R}(\bm{\xi}_{n})\right)}.$
(177) Depending on the ratio between the Lebesgue measure of
$\mathscr{R}(\bm{X}_{n})$ and $\mathscr{R}(\bm{\xi}_{n})$ one can have $K_{n}$
smaller or larger than 1.

$\displaystyle K_{0}=$
$\displaystyle\int_{\mathscr{R}(X_{0})}\int_{\mathscr{R}(X_{1})}p_{1|0}(x_{1}|x_{0})^{2}dx_{1}dx_{0}=\frac{\lambda(\mathscr{R}(X_{0}))}{\lambda(\mathscr{R}(\xi_{0}))}=\frac{1}{2}.$
(178)

Hence, both operators $\mathcal{N}(1,0)$ and $\mathcal{M}(0,1)$ are
contractions (Proposition 17). On the other hand,

$\displaystyle
K_{1}=\frac{\lambda(\mathscr{R}(X_{1}))}{\lambda(\mathscr{R}(\xi_{1}))}=1+\frac{\tan(4)-\tan(2)}{4}>1.$
(179)

Next, define $V$ as in Lemma 11, i.e., $V=[-3,3]$. Clearly, both
$\mathscr{R}(X_{0})$ and $\mathscr{R}(X_{1})$ are subsets of $V$. If we
integrate the conditional PDF shown in Figure 7 in $V\times V$ we obtain

$\int_{V}\int_{V}p_{1|0}(x_{1}|x_{0})^{2}dx_{1}dx_{0}=\frac{\lambda(V)}{\lambda(\mathscr{R}(\xi_{0}))}=\frac{3}{2}.$
(180)

### Random noise can induce contractions

In this section we prove a result on neural networks perturbed by random noise
of increasing amplitude which states that it is possible to make both
operators $\mathcal{N}$ and $\mathcal{M}$ in (149) contractions202020An linear
operator is called a contraction if its operator norm is smaller than one. if
the noise is properly chosen. To this end, we begin with the following lemma.

######  Lemma 19.

Let

$\left\|\rho_{n}\right\|^{2}_{L^{2}(\mathscr{R}(\bm{\xi}_{n}))}=\int_{\mathscr{R}(\bm{\xi}_{n})}\rho_{n}(\bm{x})^{2}d\bm{x}.$
(181)

If

$\left\|\rho_{n}\right\|^{2}_{L^{2}(\mathscr{R}(\bm{\xi}_{n}))}\leq\frac{\kappa}{\lambda(\Omega_{n})}\qquad
0\leq\kappa<1$ (182)

then $\mathcal{M}(n,n+1)$ and $\mathcal{N}(n+1,n)$ are contractions. The
condition (182) is independent of the neural network weights.

###### Proof.

The proof follows from equation (165).

∎

Hereafter we specialize Lemma 19 to neural network perturbed by uniformly
distributed random noise.

######  Proposition 20.

Let $\\{\bm{\xi}_{0},\ldots,\bm{\xi}_{L-1}\\}$ be independent random vectors.
Suppose that the components of each $\bm{\xi}_{n}$ are zero-mean i.i.d.
uniform random variables with range $[-b_{n},b_{n}]$ ($b_{n}>0$). If

$b_{0}\geq\frac{1}{2}\left(\frac{\lambda(\Omega_{0})}{\kappa}\right)^{1/N}\qquad\text{and}\qquad
b_{n}\geq\frac{b_{n-1}+1}{\kappa^{1/N}}\qquad n=1,\ldots,L-1,$ (183)

where $\Omega_{0}$ is the domain of the neural network input, $0\leq\kappa<1$,
and $N$ is the number of neurons in each layer, then both operators
$\mathcal{M}(n,n+1)$ and $\mathcal{N}(n+1,n)$ defined in (149) are
contractions for all $n=0,\ldots,L-1$, i.e., their norm can be bounded by a
constant $K_{n}\leq\kappa$, independently of the weights of the neural
network.

###### Proof.

If $\bm{\xi}_{n}$ is uniformly distributed then from (164) we have that

$K_{n}=\frac{\lambda\left(\mathscr{R}(\bm{X}_{n})\right)}{\lambda\left(\mathscr{R}(\bm{\xi}_{n})\right)}.$
(184)

By using Lemma 11 we can bound $K_{n}$ as

$K_{n}\leq\left(\frac{1+b_{n-1}}{b_{n}}\right)^{N},$ (185)

where $N$ is the number of neurons in each layer of the neural network.
Therefore, if $b_{n}\geq(b_{n-1}+1)/\kappa^{1/N}$ ($n=1,\ldots,L-1$) we have
that $K_{n}$ is bounded by a quantity $\kappa$ smaller than one. Regarding
$b_{0}$, we notice that

$K_{0}=\frac{\lambda\left(\mathscr{R}(\bm{X}_{0})\right)}{\lambda\left(\mathscr{R}(\bm{\xi}_{0})\right)}=\frac{\lambda(\Omega_{0})}{(2b_{0})^{N}},$
(186)

where $\Omega_{0}$ is the domain of the neural network input. Hence, if
$b_{0}$ satisfies (184) then $K_{0}\leq\kappa$.

∎

One consequence of Proposition 20 is that the $L^{2}$ norm of the neural
network output decays with both the number of layers and the number of neurons
if the noise amplitude from one layer to the next increases as in (183). For
example, if we represent the input-output map as a sequence of conditional
expectations (see (25)), and set $u(\bm{x})=\bm{\alpha}^{T}\bm{x}$ (linear
output) then we have

$q_{0}(\bm{x})=\mathcal{M}(0,1)\mathcal{M}(1,2)\cdots\mathcal{M}(L-1,L)(\bm{\alpha}^{T}\bm{x}).$
(187)

By iterating the inequalities (183) in Proposition 20 we find that

$b_{n}\geq\frac{1}{2\kappa^{n/N}}\left(\frac{\lambda(\Omega_{0})}{\kappa}\right)^{1/N}+\sum_{k=1}^{n}\frac{1}{\kappa^{k/N}}\qquad
n=0,\ldots,L-1,$ (188)

In Figure 8 we plot the lower bound at the right hand side of (188) for
$\kappa=0.2$ and $\kappa=10^{-4}$ as a function of the number of neurons
($N$).

$\kappa=0.2$ $\kappa=0.2$

$\kappa=10^{-4}$ $\kappa=10^{-4}$

Figure 8: Lower bound on the coefficients $b_{n}$ defined in (188) for
$\lambda(\Omega_{0})=1$ as a function of the number of neurons $N$ and number
of layers of the neural network. With such values of $b_{n}$ the operator
$\mathcal{G}(L-n+1,L-n)$ is a contraction satisfying
$\left\|\mathcal{G}(L-n+1,L-n)\right\|^{2}\leq\kappa$. Shown are results for
$\kappa=0.2$ and $\kappa=10^{-4}$ (contraction index).

With $b_{n}$ given in (188) we have that the operator norms of
$\mathcal{M}(n,n+1)$ and $\mathcal{N}(n+1,n)$ ($n=0,\ldots,L-1$) are bounded
exactly by $\kappa$ (see Lemma 19). Hence, by taking the $L^{2}$ norm of
(187), and recalling that $\left\|\mathcal{M}(n,n+1)\right\|^{2}\leq\kappa$ we
obtain

$\displaystyle\left\|q_{0}\right\|^{2}_{L^{2}(\Omega_{0})}\leq$ $\displaystyle
Z^{2}\left\|\bm{\alpha}\right\|^{2}_{2}\kappa^{L},$ (189)

where212121In equation (189) we used the Cauchy-Schwarz inequality
$\left\|\bm{\alpha}^{T}\bm{x}\right\|^{2}_{L^{2}(\mathscr{R}(\bm{X}_{L}))}\leq
Z^{2}\left\|\bm{\alpha}\right\|_{2}^{2}.$ (190)

$Z^{2}=\sum_{k=1}^{N}\int_{\mathscr{R}(\bm{X}_{L})}x_{k}^{2}d\bm{x}\qquad\text{and}\qquad\left\|\bm{\alpha}\right\|^{2}_{2}=\sum_{k=1}^{N}\alpha_{k}^{2}.$
(191)

The inequality (189) shows that the 2-norm of the vector of weights
$\bm{\alpha}$ must increase exponentially fast with the number of layers $L$
if we chose the noise amplitude as in (188). As shown in the following Lemma,
the growth rate of $b_{n}$ that guarantees that both $\mathcal{M}$ and
$\mathcal{N}$ are contractions is linear (asymptotically with the number of
neurons).

######  Lemma 21.

Consider a neural network satisfying the hypotheses of Proposition 20. Then,
in the limit of an infinite number of neurons ($N\rightarrow\infty$), the
noise amplitude (188) satisfies

$\lim_{N\rightarrow\infty}b_{n}=\frac{1}{2}+n,$ (192)

independently of the contraction factor $\kappa$ and the domain $\Omega_{0}$.
This means that for a finite number of neurons the noise amplitude $b_{n}$
that guarantees that $\|\mathcal{M}(n,n+1)\|\leq\kappa$ is bounded from below
($\kappa<1$) or from above ($\kappa>1$) by a function that increases linearly
with the number of layers.

###### Proof.

The proof follows by taking the limit of (188) for $N\rightarrow\infty$. ∎

An example: Set $k=10^{-4}$, $L=4$ (four layers) and $N=10$ neurons per layer.
The factor $\kappa^{L}$ in (189) is $10^{-16}$. If we are interested in
representing a $d$-dimensional function in the unit cube
$\Omega_{0}=[0,1]^{d}$ then we have222222In fact,
$\sum_{i=1}^{N}\int_{\mathscr{R}(\bm{X}_{4})}x_{i}^{2}d\bm{x}\leq\sum_{i=1}^{N}\int_{\Omega_{4}}x_{i}^{2}d\bm{x}=\sum_{i=1}^{N}\underbrace{\int_{-1-b_{3}}^{1+b_{3}}\cdots\int_{-1-b_{3}}^{1+b_{3}}}_{\text{$N$
times}}x_{i}^{2}d\bm{x}=\frac{2^{N}N(1+b_{3})^{N+2}}{3}.$ (193)

$\lambda(\Omega_{0})=1\qquad Z^{2}\leq N\frac{2^{N}(1+b_{3})^{N+2}}{3},$ (194)

where $b_{3}=3.684$ (see Figure 8 for $N=10$). Hence, the norm of the output
(189) is bounded by

$\displaystyle\left\|q_{L}\right\|_{L^{2}(\Omega_{0})}\leq$ $\displaystyle
C\left\|\bm{\alpha}\right\|_{2}\qquad
C=10^{-16}\sqrt{10\frac{2^{10}(1+b_{3})^{12}}{3}}\simeq 6.17\times 10^{-11}.$
(195)

This means that the 2-norm of the output coefficients $\bm{\alpha}$ has to be
of the order of $10^{11}$ to represent, e.g., a two-dimensional function of
norm about one on the square $\Omega_{0}=[0,1]^{2}$.

##  Appendix B Training over weights versus training over noise

By adding random noise to the output of each layer in a neural network we are
essentially adding an infinite number of degrees of freedom to our system.
This allows us to rethink the process of training the neural network from a
probabilistic perspective. In particular, instead of optimizing a performance
metric232323In a supervised learning setting the neural network weights are
usually determined by minimizing a dissimilarity measure between the output of
the network and a target function. Such measure may be an entropy measure, the
Wasserstein distance, the Kullback–Leibler divergence, or other measures
defined by classical $L^{p}$ norms. relative to the neural network weights
$\bm{w}=\\{\bm{w}_{0},\bm{w}_{1},\ldots,\bm{w}_{L-1}\\}$ for fixed noise, we
can now optimize the transition density242424 In a deterministic setting, the
transition density for a neural network model of the form
$\bm{X}_{n+1}=\bm{F}(\bm{X}_{n},\bm{w}_{n})$ is simply
$p_{n+1|n}(\bm{x}_{n+1}|\bm{x}_{n})=\delta\left(\bm{x}_{n+1}-\bm{F}(\bm{x}_{n},\bm{w}_{n})\right),$
(196) where $\delta(\cdot)$ is the Dirac delta function. Such density does not
have any degree of freedom other than $\bm{w}_{n}$. On the other hand, in a
stochastic setting we are free to choose the PDF of $\bm{\xi}_{n}$. For a
neural network model of the form
$\bm{X}_{n+1}=\bm{F}(\bm{X}_{n},\bm{w}_{n})+\bm{\xi}_{n}$ the transition
density has the form
$p_{n+1|n}(\bm{x}_{n+1}|\bm{x}_{n})=\rho_{n}\left(\bm{x}_{n+1}-\bm{F}(\bm{x}_{n},\bm{w}_{n})\right),$
(197) where $\rho_{n}(\bm{\xi})$ is the PDF of $\bm{\xi}_{n}$. We are clearly
free to choose the functional form of $\rho_{n}$.
$p_{n+1|n}(\bm{x}_{n+1}|\bm{x}_{n})$. Clearly, such transition density is
connected to the neural network weights, e.g., by equation (39). Hence, if we
prescribe the PDF of the random noise, i.e., $\rho_{n}(\cdot)$ in (39), then
the transition density $p_{n+1|n}$ is uniquely determined by the functional
form of the activation function $\bm{F}$, and by the weights $\bm{w}_{n}$. On
the other hand, we can optimize $\rho_{n}$ (probability density of the random
noise $\bm{\xi}_{n}$) while keeping the weights $\bm{w}_{n}$ fixed. As we
shall see hereafter, this process opens the possibility approximately encode
and decode secret messages in a fully trained neural network using random
noise.

### Encoding secret messages using random noise

An interesting question is whether random noise added to the output of each
layer in the neural network can enhance features of the output, or allow us to
encode/decode secret signals in the network. The interaction between random
noise and the nonlinear dynamics modeled by the network can yield indeed many
surprising results. For example, in stochastic resonance [36, 50] it is well
known that random noise added to a properly tuned bi-stable system can induce
a peak in the Fourier power spectrum of the output - hence effectively
amplifying the signal. Similarly, random noise added to a neural network can
have remarkable effects. In particular, it allows us to re-purpose (to some
extent) a previously trained network by hiding a secret signal in it, which
can be approximately encoded and decoded by using random noise. To this end,
it is sufficient to add random noise with appropriately selected probability
density function to the output of some (or all) layers, and then train the
expectation of the output over the PDF of the noise.

To describe this method, suppose that we are given a fully trained
deterministic neural network with only two layers, and weights chosen to
represent an input-output map defined on some domain
$\Omega\subseteq\mathbb{R}^{d}$. In the absence of noise we can write the
output of the neural network as

$q_{2}(\bm{x})=\bm{\alpha}^{T}\bm{F}(\bm{F}_{0}(\bm{x},\bm{w}_{0}),\bm{w}_{1})$
(198)

where $\\{\bm{\alpha},\bm{w}_{0},\bm{w}_{1}\\}$ can be optimized to minimize
the distance between $q_{2}(\bm{x})$ and a given target function $f(\bm{x})$
($\bm{x}\in\Omega$) . Injecting noise $\bm{\xi}_{0}$ in the output of the
first layer yields the input-output map

$h_{2}(\bm{x})=\bm{\alpha}^{T}\int\bm{F}(\bm{y}+\bm{F}_{0}(\bm{x},\bm{w}_{0}),\bm{w}_{1})\rho_{0}\left(\bm{y}\right)d\bm{y},$
(199)

where $\\{\bm{\alpha},\bm{w}_{1},\bm{w}_{0}\\}$ here are fixed, and $\rho_{0}$
is the PDF of $\bm{\xi}_{0}$. Equation (199) resembles a Fredholm integral
equation of the first kind. In fact, it can be written as

$h_{2}(\bm{x})=\int\kappa_{2}(\bm{x},\bm{y})\rho_{0}\left(\bm{y}\right)d\bm{y},$
(200)

where

$\kappa_{2}(\bm{x},\bm{y})=\bm{\alpha}^{T}\bm{F}(\bm{y}+\bm{F}_{0}(\bm{x},\bm{w}_{0}),\bm{w}_{1}).$
(201)

However, differently from standard Fredholm equations of the first kind, here
we have $\bm{x}\in\Omega\subseteq\mathbb{R}^{d}$ while
$\bm{y}\in\mathbb{R}^{N}$, i.e., the integral operator with kernel
$\kappa_{2}$ maps functions in $N$ variables into functions in $d$ variables.
We are interested in finding a PDF $\rho_{0}(\bm{y})$ that solves (199) for a
given function $h_{2}(\bm{x})$. In other words, we are re-purposing the neural
network (198) with output $q_{2}(\bm{x})\simeq f(\bm{x})$ to approximate now a
different function $h_{2}(\bm{x})\simeq v(\bm{x})$, without modifying the
weights $\\{\bm{\alpha},\bm{w}_{1},\bm{w}_{0}\\}$ but rather simply adding
noise $\bm{\xi}_{0}$ and averaging the output over the PDF $\rho_{0}$ of the
noise (Eq. (200)). Equation (200) is unfortunately ill-posed in the space of
probability distributions. In other words, for a given kernel $\kappa_{2}$ and
a given target $q_{2}$ there is (in general) no PDF $\rho_{0}$ that satisfies
(200) exactly. However, one can proceed by optimization. For instance,
$\rho_{0}$ can be determined by solving the constrained least squares
problem252525The optimization problem (202) is a quadratic program with linear
constraints if we represent $\rho_{0}$ in the span of a basis made of positive
functions, e.g., Gaussian kernels [2].

$\rho_{0}=\operatorname*{\arg\\!\min}_{\rho}\left\|h_{2}(\bm{x})-\int\kappa_{2}(\bm{x},\bm{y})\rho(\bm{y})d\bm{y}\right\|_{L^{2}(\Omega)}\qquad\text{subject
to}\qquad\|\rho\|_{L^{1}(\mathbb{R}^{N})}=1\qquad\rho\geq 0.$ (202)

#### Encoding secret messages using random shifts

A related setting is one in which we introduce a random shift in the input of
the neural network rather than adding noise262626Note that if we do not have
access to the output of each layer in the neural network, we can just
introduce random shifts in the input. This is the setting one could use to
encode secret messages in pre-trained neural networks in which the user is
allowed to modify the input and just observe the output., i.e.,
$\bm{x}\rightarrow\bm{x}-\bm{\xi}$, and average the output over the PDF $\rho$
of random vector $\bm{\xi}$. This yields

$h_{2}(\bm{x})=\int q_{2}(\bm{x}-\bm{y})\rho\left(\bm{y}\right)d\bm{y},$ (203)

where the convolution kernel

$q_{2}(\bm{x})=\bm{\alpha}^{T}\bm{F}(\bm{F}_{0}(\bm{x},\bm{w}_{0}),\bm{w}_{1}).$
(204)

represents the output of the two layer network without noise trained on a
given function $f(\bm{x})$, i.e., the same function as in (198)). Clearly,
equation (199) is the expectation of the noiseless neural network output
$q_{2}(\bm{x})$ under a random shift with PDF $\rho(\bm{y})$, i.e.,

$h_{2}(\bm{x})=\mathbb{E}_{\bm{y}\sim\rho(\bm{y})}\left[q_{2}(\bm{x}-\bm{y})\right]$
(205)

To solve (203) approximately in the space of PDF one can proceed by
optimization, i.e., solve the quadratic problem (202) with
$\kappa_{2}(\bm{x},\bm{y})$ replaced by $q_{2}(\bm{x}-\bm{y})$.

Remark: The integral equation (203) can, in principle, be solved analytically
in $\mathbb{R}^{d}$ (or in any $d$-dimensional torus) by using Fourier
transforms (provided the Fourier transform of $\ell_{2}$ is nonzero). This
yields272727In equation (206) we assumed that
$\mathcal{F}[q_{2}(\bm{x})](\bm{\zeta})\neq 0$ for all
$\zeta\in\mathbb{R}^{d}$.

$\rho\left(\bm{y}\right)=\int\frac{\mathcal{F}[h_{2}(\bm{x})](\bm{\zeta})}{\mathcal{F}[q_{2}(\bm{x})](\bm{\zeta})}e^{2\pi
i\bm{\zeta}\cdot\bm{x}}d\bm{\zeta},$ (206)

where $\mathcal{F}[\cdot]$ denotes the multivariate Fourier transform operator

$\mathcal{F}[f(\bm{x})](\bm{\zeta})=\int f(\bm{x})e^{-2\pi
i\bm{x}\cdot\bm{\zeta}}d\bm{x}.$ (207)

However, the function $\rho(\bm{y})$ defined in (206) is in general not a PDF,
except in those cases where both $q_{2}$ and $h_{2}$ are PDFs, and the PDF of
$h_{2}$ is essentially the PDF of the sum of two independent random variables
[38].

## References

  * [1] F. Bartolucci, E. De Vito, L. Rosasco, and S. Vigogna.  Understanding neural networks with reproducing kernel Banach spaces.  ArXiv:2109.09710, pages 1–42, 2021.
  * [2] Z. I. Botev, J. F. Grotowski, and D. P. Kroese.  Kernel density estimation via diffusion.  Ann. Stat., 38(5):2916–2957, 2010.
  * [3] C. Brennan and D. Venturi.  Data-driven closures for stochastic dynamical systems.  J. Comp. Phys., 372:281–298, 2018.
  * [4] M. Chen, X. Li, and C. Liu.  Computation of the memory functions in the generalized langevin models for collective dynamics of macromolecules.  J. Chem. Phys, 141(6):064112, 2014.
  * [5] H. Cho, D. Venturi, and G. E. Karniadakis.  Statistical analysis and simulation of random shocks in Burgers equation.  Proc. R. Soc. A, 2171(470):1–21, 2014.
  * [6] A. J. Chorin, O. H. Hald, and R. Kupferman.  Optimal prediction and the Mori–Zwanzig representation of irreversible processes.  Proc. Natl. Acad. Sci., 97(7):2968–2973, 2000.
  * [7] W. Chu and X. Li.  The Mori–Zwanzig formalism for the derivation of a fluctuating heat conduction model from molecular dynamics.  Comm. Math. Sci., 17(2), 2019.
  * [8] G. Ciccotti and J.-P. Ryckaert.  On the derivation of the generalized Langevin equation for interacting Brownian particles.  J. Stat. Phys., 26(1):73–82, 1981.
  * [9] J. M. Dominy and D. Venturi.  Duality and conditional expectations in the Nakajima-Mori-Zwanzig formulation.  J. Math. Phys., 58(8):082701, 2017.
  * [10] W. E.  A proposal on machine learning via dynamical systems.  Commun. Math. Stat., 5:1–10, 2017.
  * [11] W. E, J. Han, and Q. Li.  A mean-field optimal control formulation of deep learning.  Res. Math. Sci., 6:1–41, 2019.
  * [12] S. Gibert and A. Mukherjea.  Nonnegative idempotent kernels.  J. Math. Anal. Appl., 135(1):326–341, 1988.
  * [13] L. Gonon, L. Grigoryeva, and J.-P. Ortega.  Risk bounds for reservoir computing.  JMLR, 21(240):1–61, 2020.
  * [14] J. Harlim, S. W. Jiang, S. Liang, and H. Yang.  Machine learning for prediction with missing dynamics.  J. Comput. Phys., 428:109922, 2021.
  * [15] K. He, X. Zhang, S. Ren, and J. Sun.  Deep residual learning for image recognition.  In Proc. IEEE Comput. Soc. Conf. Comput. Vis. Pattern Recognit., pages 770–778, 2016.
  * [16] K. He, X. Zhang, S. Ren, and J. Sun.  Identity mappings in deep residual networks.  In ECCV, pages 630–645. Springer, 2016.
  * [17] C. Hijón, P. Español, E. Vanden-Eijnden, and R. Delgado-Buscalioni.  Mori–Zwanzig formalism as a practical computational tool.  Faraday discussions, 144:301–322, 2010.
  * [18] C. Hijón, M. Serrano, and P. Español.  Markovian approximation in a coarse-grained description of atomic systems.  J. Chem. Phys., 125:204101, 2006.
  * [19] S. Izvekov and G. A. Voth.  Modeling real dynamics in the coarse-grained representation of condensed phase systems.  J. Chem. Phys., 125:151101–151104, 2006.
  * [20] G. J. O. Jameson and A. Pinkus.  Positive and minimal projections in function spaces.  J. Approx. Theory, 37:182–195, 1983.
  * [21] D. Kauzlarić, J. T. Meier, P. Español, A. Greiner, and S. Succi.  Markovian equations of motion for non-Markovian coarse-graining and properties for graphene blobs.  New J. Phys., 15(12):125015, 2013.
  * [22] A. I. Khuri.  Applications of Dirac’s delta function in statistics.  Int. J. Math. Educ. Sci. Technol., 35(2):185–195, 2004.
  * [23] I. Kobyzev, S. J. D. Prince, and M. A. Brubaker.  Normalizing flows: An introduction and review of current methods.  IEEE transactions on pattern analysis and machine intelligence, 43(11):3964–3979, 2020.
  * [24] A. Lasota and M. C. Mackey.  Chaos, fractals and noise: stochastic aspects of dynamics.  Springer–Verlag, second edition, 1994.
  * [25] H. Lei, N. A. Baker, and X. Li.  Data-driven parameterization of the generalized Langevin equation.  Proc. Natl. Acad. Sci., 113(50):14183–14188, 2016.
  * [26] Q. Li, L. Chen, C. Tai, and W. E.  Maximum principle based algorithms for deep learning.  JMLR, 18:1–29, 2018.
  * [27] Q. Li, T. Lin, and Z. Shen.  Deep learning via dynamical systems: An approximation perspective.  J. Eur. Math. Soc., (published online first), 2022.
  * [28] Q. Liao and T. Poggio.  Bridging the gaps between residual learning, recurrent neural networks and visual cortex.  arXiv:1604.03640, 2016.
  * [29] K. K. Lin and F. Lu.  Data-driven model reduction, Wiener projections, and the Koopman-Mori-Zwanzig formalism.  J. Comput. Phys., 424:109864, 2021.
  * [30] F. Lu, K. Lin, and A. Chorin.  Comparison of continuous and discrete-time data-based modeling for hypoelliptic systems.  Comm. Appl. Math. Comput. Sci., 11(2):187–216, 2016.
  * [31] Y. Lu, Zhong A, Q. Li, and B. Dong.  Beyond finite layer neural networks: Bridging deep architectures and numerical differential equations.  arXiv:1710.10121, 2017.
  * [32] L. Ma, X. Li, and C. Liu.  The derivation and approximation of coarse-grained dynamics from Langevin dynamics.  J. Chem. Phys., 145(20):204117, 2016.
  * [33] Lina Ma, Xiantao Li, and Chun Liu.  Coarse-graining langevin dynamics using reduced-order techniques.  Journal of Computational Physics, 380:170–190, 2019.
  * [34] E. Minguzzi.  The equality of mixed partial derivatives under weak differentiability conditions.  eal Anal. Exch., 40(1):81–98, 2014/2015.
  * [35] H. Mori.  Transport, collective motion, and Brownian motion.  Prog. Theor. Phys., 33(3):423–455, 1965.
  * [36] D. Nozaki, D. J. Mar, P. Grigg, and J. J. Collins.  Effects of colored noise on stochastic resonance in sensory neurons.  Phys. Rev. Lett., 82(11):2402–2405, 1999.
  * [37] J. Palis and W. de Melo.  Geometric theory of dynamical systems: an introduction.  Springer–Verlag, 1982.
  * [38] A. Papoulis.  Probability, random variables and stochastic processes.  McGraw-Hill, third edition, 1991.
  * [39] R. Parhi and R. D. Nowak.  Banach space representer theorems for neural networks and ridge splines.  JMLR, (22):1–40, 2021.
  * [40] J. Price and P. Stinis.  Renormalized reduced order models with memory for long time prediction.  Multiscale Model. Simul., 17(1):68–91, 2019.
  * [41] D. Rezende and S. Mohamed.  Variational inference with normalizing flows.  In International conference on machine learning, pages 1530–1538. PMLR, 2015.
  * [42] H. Risken.  The Fokker-Planck equation: methods of solution and applications.  Springer-Verlag, second edition, 1989.  Mathematics in science and engineering, vol. 60.
  * [43] A. Sherstinsky.  Fundamentals of recurrent neural network (RNN) and long short-term memory (LSTM) network.  Phys. D: Nonlinear Phenom., 404:132306, 2020.
  * [44] S. Shreve.  Stochastic calculus for finance II: continuous-time models.  Springer, first edition, 2004.
  * [45] E. G. Tabak and E. Vanden-Eijnden.  Density estimation by dual ascent of the log-likelihood.  Comm. Math. Sci., 8(1):217–233, 2010.
  * [46] N. Črnjarić-Žic, S. Maćešić, and I. Mezić.  Koopman operator spectrum for random dynamical systems.  J. Nonlinear Sci., 30:2007–2056, 2020.
  * [47] D. Venturi.  The numerical approximation of nonlinear functionals and functional differential equations.  Physics Reports, 732:1–102, 2018.
  * [48] D. Venturi, H. Cho, and G. E. Karniadakis.  The Mori-Zwanzig approach to uncertainty quantification.  In R. Ghanem, D. Higdon, and H. Owhadi, editors, Handbook of uncertainty quantification. Springer, 2016.
  * [49] D. Venturi and A. Dektor.  Spectral methods for nonlinear functionals and functional differential equations.  Res. Math. Sci., 8(27):1–39, 2021.
  * [50] D. Venturi and G. E. Karniadakis.  Convolutionless Nakajima-Zwanzig equations for stochastic analysis in nonlinear dynamical systems.  Proc. R. Soc. A, 470(2166):1–20, 2014.
  * [51] T. Yu, Y. Yang, D. Li, T. Hospedales, and T. Xiang.  Simple and effective stochastic neural networks.  In Proc. Innov. Appl. Artif. Intell. Conf., volume 35, pages 3252–3260, 2021.
  * [52] H. Zhang and Y. Xu.  Reproducing kernel banach spaces for machine learning.  JMLR, (10):2741–2775, 2009.
  * [53] Y. Zhu, J. M. Dominy, and D. Venturi.  On the estimation of the Mori-Zwanzig memory integral.  J. Math. Phys, 59(10):103501, 2018.
  * [54] Y. Zhu and D. Venturi.  Faber approximation of the Mori-Zwanzig equation.  J. Comp. Phys., 372:694–718, 2018.
  * [55] Y. Zhu and D. Venturi.  Generalized Langevin equations for systems with local interactions.  J. Stat. Phys., 178(5):1217–1247, 2020.
  * [56] Y. Zhu and D. Venturi.  Hypoellipticity and the Mori-Zwanzig formulation of stochastic differential equations.  J. Math. Phys., 62:1035051, 2021.
  * [57] R. Zwanzig.  Memory effects in irreversible thermodynamics.  Phys. Rev., 124(4):983, 1961.