# Limitations of neural network training due to numerical instability of
backpropagation

Clemens Karner 111University of Vienna, Faculty of Mathematics. Kolingasse
14–16, 1090 Vienna, Austria.
<EMAIL_ADDRESS>Vladimir Kazeev
††footnotemark:  222University of Vienna, Research Network Data Science.
Kolingasse 14–16, 1090 Vienna, Austria. Philipp Christian Petersen
††footnotemark:  ††footnotemark:

###### Abstract

We study the training of deep neural networks by gradient descent where
floating-point arithmetic is used to compute the gradients. In this framework
and under realistic assumptions, we demonstrate that it is _highly unlikely_
to find ReLU neural networks that maintain, in the course of training with
gradient descent, _superlinearly_ many affine pieces with respect to their
number of layers. In virtually all approximation theoretical arguments which
yield high order polynomial rates of approximation, sequences of ReLU neural
networks with _exponentially_ many affine pieces compared to their numbers of
layers are used. As a consequence, we conclude that approximating sequences of
ReLU neural networks resulting from gradient descent in practice differ
substantially from theoretically constructed sequences. The assumptions and
the theoretical results are compared to a numerical study, which yields
concurring results.

Keywords: : Round-off errors, deep neural networks, gradient descent,
numerical stability.

Mathematics Subject Classification: 35A18, 65T60, 68T10.

##  1 Introduction

Deep learning is a machine learning technique based on artificial neural
networks with a large number of layers, which are trained by gradient-based
methods. This technique has been tremendously successful in a wide range of
applications [26, 24, 43, 41]. Of particular interest for applied
mathematicians are recent developments in which deep neural networks are
applied to tasks of numerical analysis such as the numerical solution of
inverse problems [1, 34, 27, 20, 38] or of (parametric) partial differential
equations [7, 12, 39, 9, 40, 25, 29, 3].

The appeal of deep neural networks for these applications is due to their
exceptional efficiency in representing functions from several approximation
classes that underlie well-established numerical methods. In terms of
approximation accuracy with respect to the number of approximation parameters,
deep neural networks have been theoretically proven to achieve approximation
rates that are at least as good as those of finite elements [15, 35, 30],
local Taylor polynomials or splines [46, 11], wavelets [42] and, more
generally, affine systems [5].

In the sequel, we consider neural networks with the rectified-linear-unit
(ReLU) activation function, which is standard in most applications. In this
case, the neural-network approximations are piecewise-affine functions. We
point out that all state-of-the-art results on the rates of approximation with
deep ReLU neural networks that achieve higher order polynomial approximation
rates are based on explicit constructions with the number of affine pieces
growing exponentially with respect to the number of layers; see, e.g., [46,
45].

In this work, we argue that this central building block, functions with
exponentially many affine pieces, _cannot be learned with the state-of-the-art
techniques_. We show theoretically that training in floating-point arithmetic
is hindered by a bound on the number of affine pieces created during an
iterative learning process. This bound is polynomial with respect to the
number of approximation parameters and also with respect to the number of
iterations. Notably, this non-computability of functions with exponentially
many affine pieces is not derived from an abstract result of computability
theory, as the related results established in [4, 6]; instead, we identify a
concrete reason why gradient descent based on backpropagation cannot produce
these functions in floating-point arithmetic. The effect of numerical
instability is also demonstrated in numerical examples. We will give a
detailed description of our findings in Subsection 1.3.

This work is strongly influenced by [14], which shows that functions with
exponentially many affine pieces with respect to the number of layers
typically do not appear in randomly initialised neural networks. The ideas
presented in that work form the basis of our analysis.

Before presenting our main results, we recapitulate the notion of floating-
point arithmetic and furnish an example illustrating how the finite precision
of floating-point arithmetic can undermine learning due to the phenomenon
known in numerical analysis as _catastrophic cancellation_.

###  1.1 Floating-point arithmetic

Computations are almost exclusively performed in floating-point arithmetic,
which typically consists in restricting the arithmetic of real numbers to a
discrete set of the form

$\displaystyle\mathbb{M}=\big{\\{}\pm b^{e}\sum_{k=0}^{p}b^{-k}c_{k}\colon
c_{0}=1,\,c_{1},\ldots,c_{s}\in\\{0,\ldots,b-1\\},\,e\in\\{e_{\mathrm{min}},\ldots,e_{\mathrm{max}}\\}\big{\\}},$
(1.1)

extended with a few special elements, such as zero, infinity and NaN elements.
Here, $b\in\mathbb{N}$ is the radix of the floating-point arithmetic
(generally $b>1$; typically $b=2$), $p\in\mathbb{N}$ is the precision and
$e_{\mathrm{min}},e_{\mathrm{max}}\in\mathbb{Z}$ are the minimum and maximum
exponents. Naturally, $\mathbb{M}$ is not closed under basic arithmetic
operations such as addition and multiplication, which necessitates rounding.
The round-to-nearest addition and multiplication are defined so that

$|a+_{\mathbb{M}}b-(a+b)|=\min_{s\in\mathbb{M}}|s-(a+b)|\,,\qquad|a\times_{\mathbb{M}}b-a\cdot
b|=\min_{s\in\mathbb{M}}|s-ab|$ (1.2)

for all $a,b\in\mathbb{M}$. The radix $b$, the precision $p$ and the above
rounding together define the so-called _machine epsilon_
$\epsilon=\frac{1}{2}b^{1-p}$, which is an accurate upper bound for the
relative accuracy of the floating-point addition and multiplication:

$|a+_{\mathbb{M}}b-(a+b)|\leq\epsilon\,|a+b|\quad\text{and}\quad|a\times_{\mathbb{M}}b-a\cdot
b|\leq\epsilon\,|a\cdot b|$ (1.3)

for all $a,b\in\mathbb{M}$. For the double- and single-precision IEEE 754
standards, we have $\epsilon=2^{-53}\approx 1.11\cdot 10^{-16}$ and
$\epsilon=2^{-24}\approx 5.96\cdot 10^{-8}$ respectively. Note that
computation in floating point arithmetic amounts to introducing a _relative
error_ in each computation.

###  1.2 Instability of deep neural networks due to catastrophic cancellation

The bounds (1.3) allow for floating-point operations to be modelled as
perturbations of the respective infinite-precision operations of the real
arithmetic.

Consider a real feed-forward ReLU neural network with $L\in\mathbb{N}$ layers
with $d=N_{0}\in\mathbb{N}$ real inputs and $N_{j}\in\mathbb{N}$ neurons in
each layer $j\in\mathbb{N}$. The evaluation
$\phi\colon\mathbb{R}^{d}\to\mathbb{R}^{N_{L}}$ of such a neural network at
any $x^{(0)}\in\mathbb{R}^{d}$ consists of iteratively applying the
transformations

$\displaystyle x^{(j)}=\varrho\,(A_{j}\,x^{(j-1)}+b_{j})$ (1.4)

sequentially for $j=1,\ldots,L$, which produces the corresponding output
$\phi(x^{(0)})=x^{(L)}\in\mathbb{R}^{N_{L}}$. Here,
$\varrho\colon\mathbb{R}\to\mathbb{R}$ given by $\varrho(x)=\max\\{x,0\\}$ for
all $x\in\mathbb{R}$ is the ReLU activation function, which is applied to the
elements of $\mathbb{R}^{N_{j}}$ with each $j\in\\{1,\ldots,L\\}$
componentwise, whereas $A_{j}\in\mathbb{R}^{N_{j-1}\times N_{j}}$ and
$b_{j}\in\mathbb{R}^{\mathbb{N}_{j}}$ with $j\in\\{1,\ldots,L\\}$ are the
weight matrices and bias vectors.

For simplicity, let us consider the case of $b_{j}=0$ for all
$j\in\\{1,\ldots,L\\}$ and assume that $x^{(0)}$ and $A_{1},\ldots,A_{L}$ are
such that all entries of $x^{(j)}$ are nonnegative for each
$j\in\\{1,\ldots,L\\}$. Then

$\phi(x)=A_{L}\cdots A_{1}x\quad\text{for each}\quad x\in\mathbb{R}^{d}\,,$
(1.5)

i.e., the evaluation of the network reduced to the multiplication of the input
by a matrix product consisting of $L$ factors.

The instability of such products for large $L$ was demonstrated in the context
of tensor decompositions, in [2, Example 3]: the first component of the tensor
considered therein yields the matrix product (1.5) with $N_{0}=N_{L}=1$,
$N_{1}=\cdots=N_{L-1}=2$ and with weight matrices

$\displaystyle A_{1}=\begin{pmatrix}(1+\varepsilon_{1})\lambda\\\
\lambda\end{pmatrix},\qquad
A_{L}=\begin{pmatrix}(1+\varepsilon_{L})(1+\lambda^{-L})\lambda&-\lambda\end{pmatrix}\quad\text{and}\quad
A_{j}=\begin{pmatrix}(1+\varepsilon_{j})\,\lambda&0\\\ 0&\lambda\end{pmatrix}$
(1.6)

for each $j=2,\dots,L-1$, where $\lambda>1$ is a fixed parameter and
$\varepsilon_{1},\ldots,\varepsilon_{L}\geq 0$ are perturbation parameters.
For all $j\in\\{1,\ldots,L\\}$, the case of
$\varepsilon_{1}=\cdots=\varepsilon_{j-1}=\varepsilon_{j+1}=\cdots=\varepsilon_{L}=0$
and $\varepsilon_{j}\geq 0$ gives
$x^{(L)}(\varepsilon)=\phi(x^{(0)})=1+\varepsilon_{j}\,(1+\lambda^{L})x^{(0)}$
for any $x^{(0)}\in\mathbb{R}^{d}$. Considering $\varepsilon$ as a
perturbation parameter, which may be at the level of the machine epsilon
$\epsilon$, we observe that a relative perturbation of magnitude $\varepsilon$
in the $j$th weight matrix leads to the error

$x^{(L)}(\varepsilon)-x^{(L)}(0)=\varepsilon\,(1+\lambda^{L})\,x^{(L)}(0)\,,$

which implies the total loss of accuracy for any $\varepsilon>0$ whenever
$\lambda^{L}\gtrsim\varepsilon$. This effect is an immediate consequence of
the subtraction encoded by $A_{L}$, in which the two terms may be perturbed
individually (just as the first is perturbed by a multiplicative factor of
$1+\varepsilon$ and the second is not perturbed at all in the above example).
This is an archetypal example of the so-called _catastrophic cancellation_ ;
see, e.g., [17, Section 1.7].

This example, based on inaccurate matrix multiplication, can be generalised
into the following statement for deep neural networks. The precise statement
is given as Proposition A.3 in the Appendix.

######  Proposition 1.1.

Let $L\in\mathbb{N}$ and $N\in\mathbb{N}$. Let $\lambda=10^{r}\geq 1$ for
$r\in\mathbb{Z}$ and let $\mathbb{M}$ be as in (1.1) $b=10$, $p=16$ and hence
$\epsilon=5\cdot 10^{-16}$. If
$(L-3)\log_{10}(N-1)+(L-1)\log_{10}(\lambda-2\epsilon)\geq 16$, then there is
a neural network $\Phi_{\lambda,N,L}$ with $L$ layers, $N$ neurons per layer,
and all weights bounded in absolute value below by 1 and above by $\lambda$
such that

$|\Phi_{\lambda,N,L}(x)-\overline{\Phi_{\lambda,N,L}(x)}|=|\Phi_{\lambda,N,L}(x)|=|x|$

for all $x\geq 0$, where $\overline{\Phi_{\lambda,N,L}(x)}$ is the neural
network $\Phi_{\lambda,N,L}$ evaluated in floating-point arithmetic with the
machine epsilon $\epsilon$.

In the finite-precision setting, we observe that the evaluation of very small
neural networks can already lead to a relative error of $1$ compared to the
neural network with arbitrary precision. Admissible parameters such that the
conclusion of Proposition 1.1 holds are, for example, $N=65$, $\lambda=10$,
and $L=8$.

Proposition 1.1 considers only the forward application of a neural network.
However, the construction of this neural network is based on the example (1.6)
and is such that the neural network is entirely linear. By the backpropagation
rule, recalled in Definition 2.3, the reversal of the neural network
constructed in Proposition 1.1 yields a neural network such that _the accurate
computation of the gradient of the weights with respect to the output of the
neural network is infeasible in floating-point arithmetic_.

###  1.3 Contribution

In this work, we analyse the effect of inaccurate computations stemming from
floating-point arithmetic in a more moderate regime than that of Proposition
1.1. In contrast to the error amplification, that was the basis of Example
1.6, we assume the errors in the following to stem mostly from _accumulation_
of many inaccurate operations. We assume that, in the computation of the
gradient of a neural network, in each layer, a small and controlled relative
error will appear. Under this assumption, we observe that neural networks
trained with gradient descent will, with high probability, _not_ exhibit
exponentially many affine pieces with respect to the number of layers.

We describe a framework for learning in floating-point arithmetic in Section
2. There, we define a gradient descent procedure where the gradient is
computed with noisy updates.

To facilitate our main results, we make two assumptions on the gradient
descent procedure (Assumptions 2.6 and 2.7), which can be intuitively
formulated as follows.

  * (A)

We assume that the average of the reciprocals of the non-zero bias updates in
an iteration is bounded by a polynomial in the maximum number $N$ of neurons
in a single layer. Moreover, every neuron with zero output (_dead neuron_) is
assumed to remain so after the iteration. Both assumptions together
essentially require the gradient of every bias to be either zero or not too
small for most neurons. We verify this assumption numerically in Section 4.

  * (B)

_For all layers $j\in\mathbb{N}$, the derivative of each coordinate of the
preactivations, i.e., $A_{j}\,x^{(j-1)}+b_{j}$, with respect to the input of
the neural network is bounded by a uniform constant._

Under these assumptions, we prove Theorem 3.4, which can be intuitively
phrased as follows: the expected number of affine pieces that the neural
network has after a single perturbed gradient descent update has a lower bound
that is polynomial (in practice quadratic) in the number of neurons, linear in
the number of layers and inversely proportional to the level of perturbation.

The result is based on the following insight of [14]. Consider as a
prototypical neuron a function of the form $x\mapsto\varrho(h(x)-b)$, where
$h$ is a piecewise affine function of bounded variation. Then, for
$\varrho(h(\cdot)-b)$ to have more affine pieces than $h$ it is necessary that
$h(x)=b$ for some $x$ in the domain of $h$. Intuitively, for many such $x$ to
exist, $h$ has to cross the $b$-level often. If $b$ is a random variable,
then, for the event $h(x)=b$ to happen with high probability, $h$ needs to
repeatedly cross a substantial part of the region where $b$ has most of its
mass. This, however, can only happen to some degree since the variation of $h$
was assumed to be bounded. The resulting estimate of expected newly generated
affine pieces in one neuron under random biases is formalised in Lemma 3.1.
Since the floating-point operations are assumed to introduce noise into the
bias updates in a way quantified in Definition 2.3, and since Assumption (B)
implies bounded variation of the outputs of each neuron, we deduce that with
high probability the number of affine pieces in a neural network is bounded.

In a gradient descent iteration, the noise in each gradient update step also
depends on the current empirical risk. Notably, if the neural network is
initialised with zero-training error, no update happens. However, under
standard assumptions on the dynamics of the loss, we can still prove that the
number of affine pieces during the iteration does not increase quickly. This
is the content of Theorem 3.6, which we state below in a simplified and
informal version:

######  Theorem 1.2 (Informal version of Theorem 3.6).

During gradient descent, where the gradient is computed by backpropagation in
floating-point arithmetic, it holds for each iteration with probability $1/2$
that the number of affine pieces of a neural network along a given line grows
at most polynomially (almost linearly) with respect to the number of
iterations, polynomial in the number of neurons, and linear in the number of
layers.

Assumptions (A) and (B) will be discussed in detail after their formal
statement in Section 2. Since Assumption (A) relates to the distribution of
updates, we numerically study if it is satisfied in practice in Section 4.

Finally, in Section 5, we numerically analyse the influence of the magnitude
of round-off errors. We observe that, while lower numerical accuracy
influences the generation of affine pieces negatively, as claimed by our main
theorem, the effect is not as pronounced as the direct application of the
theorem may suggest. The reason is that, for the approximations involved, the
number of affine pieces is already significantly lower than our upper bounds.
Hence we conclude that numerical stability is one reason that prevents the
learning of functions with exponentially many affine pieces, but there may be
many more such reasons that prevent the number of affine pieces from reaching
the theoretically established threshold. If the number of affine pieces is
large for the initial neural network, it reduces in the course of training; in
accordance with our theoretical results, that reduction occurs faster for
lower computations accuracies.

###  1.4 Related work

This work represents a counter-point to many approximation theoretical
approaches used in the analysis of deep learning. Here, we study to what
extent piecewise-affine functions with a number of affine pieces that is
exponential with respect to the number of layers can be learned by deep neural
networks using gradient descent algorithms. Functions with high numbers of
affine pieces are necessary in the approximation of the square function [46],
which forms the basis for many approximation results in the literature, e.g.,
[47, 35, 30, 37, 40, 11].

The primary motivation for the approach that we follow is taken from [14].
There it is shown that randomly initialised neural networks typically have few
affine pieces. The main idea of that work is to show that to generate many
affine pieces, the bias vectors in neural networks need to be very accurately
chosen. This is unlikely to be achieved if the biases are initialised
randomly. This very idea is also the basis behind our analysis. We cannot,
however, simply apply the main results of [14] since analysing the whole
iterative procedure of gradient descent requires keeping track of the
interdependence of the random variables modelling the weights and biases of
the network and estimating their variance throughout the process.

Floating-point arithmetic and its effect on iterated computations have been
studied in the literature before. The example given by (1.5)–(1.6) is derived
from the study of the stability of long tensor factorisations in [2]. The
explicit constructions of low-rank tensor approximations presented in [23, 21,
31, 22] all, depending on the specific implementation, may suffer from the
instability similar to that demonstrated in (1.5)–(1.6).

In [28], an empirical study is presented that finds that the accuracy affects
the overall performance of a neural-network classifier. In this context, also
an effect of the number of layers is measured, showing that deeper neural
networks are more sensitive to low accuracies. Another analysis is given in
[44], where it is shown that numerical inaccuracies introduce stochasticity
into the trajectories of trained neural networks, which is similar to noise in
stochastic gradient descent. In the context of fixed-point arithmetic, [10]
showed that in many classification problems, one can restrict the computation
accuracy significantly without harming the classification accuracy. There is
also a vast body of literature studying how neural networks can be efficiently
quantised to speed up computations or minimise storage requirements [48, 13,
19, 18]. Typically it is found that moderate quantisation or specifically
designed quantisation does not harm classification accuracies. To the best of
our knowledge, the effect on the generation of many affine pieces has not been
studied.

There have also been more abstract studies showing that on digital hardware or
in finite precision computations, there are learning scenarios that cannot be
solved to any reasonable accuracy; see [6, 4].

##  2 Framework

In this section, we introduce the framework for this article. We start by
formally introducing neural networks in Definition 2.1. Thereafter, we fix the
notion of _number of affine pieces_ of a neural network in Definition 2.2. To
study the effect of numerical accuracy on the convergence and performance of
gradient descent based learning, we define the perturbed gradient descent
update in Definition 2.3 and the full perturbed gradient descent iteration in
Definition 2.5.

To facilitate the results in the sequel, we make Assumptions 2.6 and 2.7 on
the gradient descent dynamics. Assumption 2.6 requires the absolute value of
the derivative of the objective function with respect to each of the biases to
either vanish or be bounded from below. In the event that the derivative
vanishes, we assume that the associated neuron is dead, i.e., its output is
constant on the whole domain. Assumption 2.7 stipulates that the output of
each neuron should have a bounded derivative with respect to the input
variable. We will discuss the sensibility of these assumptions at the end of
this section.

We start by defining a neural network. Here we focus only on the most widely
used activation function, the ReLU, $\varrho(x)\coloneqq\max\\{0,x\\}$, for
$x\in\mathbb{R}$.

######  Definition 2.1 ([37, 35]).

Let $d,L\in\mathbb{N}$. A _neural network (NN) with input dimension $d$ and
$L$ layers_ is a sequence of matrix-vector tuples

$\Phi\coloneqq\big{(}(A_{1},b_{1}),(A_{2},b_{2}),\dots,(A_{L},b_{L})\big{)},$

where $N_{0}\coloneqq d$ and $N_{1},\dots,N_{L}\in\mathbb{N}$, and where
$A_{j}\in\mathbb{R}^{N_{j}\times N_{j-1}}$ and $b_{j}\in\mathbb{R}^{N_{j}}$
for $j=1,...,L$. The number $N_{L}$ is referred to as the _output dimension_.

For $j^{\prime}\in\\{1,\dots,L-1\\}$, the _subnetwork from layer $1$ to layer
$j^{\prime}$ of $\Phi$_ is the sequence of matrix-vector tuples

$\Phi_{j^{\prime}}\coloneqq\big{(}(A_{1},b_{1}),(A_{2},b_{2}),\dots,(A_{j^{\prime}},b_{j^{\prime}})\big{)}.$

For a NN $\Phi$ and a domain $\Omega\subset\mathbb{R}^{d}$, we define the
associated _realisation of the NN $\Phi$_ as

$\mathrm{R}(\Phi):\Omega\to\mathbb{R}^{N_{L}}:x\mapsto
x^{(L)}:=\mathrm{R}(\Phi)(x),$

where the output $x^{(L)}\in\mathbb{R}^{N_{L}}$ results from

$\begin{split}x^{(0)}&\coloneqq x,\\\
x^{(j)}&\coloneqq\varrho(A_{j}\,x^{(j-1)}+b_{j})\quad\text{ for
}j=1,\dots,L-1,\\\ x^{(L)}&\coloneqq A_{L}\,x^{(L-1)}+b_{L}.\end{split}$ (2.1)

Here $\varrho$ is understood to act component-wise on vector-valued inputs,
i.e., for $y=(y^{1},\dots,y^{m})\in\mathbb{R}^{m}$,
$\varrho(y):=(\varrho(y^{1}),\dots,\varrho(y^{m}))$. We call $N(\Phi)\coloneqq
d+\sum_{j=1}^{L}N_{j}$ the _number of neurons of the NN_ $\Phi$,
$L(\Phi)\coloneqq L$ the _number of layers_ or _depth_. In addition, we refer
to $(d,N_{1},\dots,N_{L})$ as the _architecture of $\Phi$_.

Since the ReLU is piecewise linear, it is not hard to see that the realisation
of a NN is always a piecewise affine function. A widely-used tool to study
deep NNs is by counting the number of the affine pieces. To formalise this, we
present a definition of the number of affine pieces of a piecewise affine
function.

######  Definition 2.2.

Let $d\in\mathbb{N}$. For a piecewise affine function $f:\Omega\to\mathbb{R}$
with $\Omega\subset\mathbb{R}^{d}$ and for any $\Omega^{\prime}\subset\Omega$,
we define the _number of affine pieces of_ $f$ on $\Omega^{\prime}$ as the
smallest integer $P$ such that there exist open sets
$(U_{i})_{i=1}^{P}\subset\mathbb{R}^{d}$ satisfying

  * •

$\bigcup_{i=1}^{P}\overline{U_{i}}=\mathbb{R}^{d}$,

  * •

$f_{|\Omega^{\prime}\cap U_{i}}$ is affine for all $i\in\\{1,\ldots,P\\}$.

We denote the number of affine pieces of $f$ by
$\mathrm{pieces}(f,\Omega^{\prime})$.

Next, we introduce a model for the training process of deep NNs. It is
customary that NNs are trained by minimising an empirical risk function using
gradient descent333In practice, stochastic gradient descent is often used, in
which the gradient is only computed over random batches of training data
instead of the complete dataset. This introduces additional randomness into
the gradient dynamics and, in principle, should allow for our main results to
be strengthened. However, since this randomness does not stem from the
numerical accuracy, which is our main focus in the present work, we confine it
to the case of standard gradient descent.  over the weights. Assume we are
given a loss function
$\ell\colon\mathbb{R}^{q}\times\mathbb{R}^{q}\to\mathbb{R}^{+}$, for
$q\in\mathbb{N}$, which could, for example, be the _square loss_
$\ell(y,y^{\prime})=\|y-y^{\prime}\|^{2}$. Then, for $M\in\mathbb{N}$, we
define
$\widehat{\mathcal{R}}\colon(\mathbb{R}^{q})^{M}\times(\mathbb{R}^{q})^{M}\to\mathbb{R}$
by

$\widehat{\mathcal{R}}(Y,\hat{Y})=\frac{1}{M}\sum_{i=1}^{M}\ell(Y_{i},\hat{Y}_{i})$

for all
$Y=(Y_{i})_{i=1}^{M}\subset\mathbb{R}^{q},\hat{Y}=(\widehat{Y}_{i})_{i=1}^{M}\subset\mathbb{R}^{q}$.

For NNs $\Phi$ with a fixed architecture, input dimension $d$, and output
dimension $N_{L}=1$, we formalise the procedure to minimise
$\widehat{\mathcal{R}}(Y,\mathrm{R}(\Phi)(X))$ for given training data
$X=({x^{(L-1)}}i)_{i=1}^{M}\subset\Omega,Y=(Y_{i})_{i=1}^{M}\subset\mathbb{R}$,
where $\mathrm{R}(\Phi)(X)=(\mathrm{R}(\Phi)(X_{i}))_{i=1}^{M}$.

######  Definition 2.3 (Gradient descent update of a neural network).

Let $d\in\mathbb{N}$, $\Omega\subset\mathbb{R}^{d}$, and let $\Phi$ be a NN of
depth $L\in\mathbb{N}$, input dimension $d$, and output dimension $1$. Let
$A_{j},b_{j},N_{j}$ be as in Definition 2.1. Further, let
$\varepsilon=(\varepsilon_{j})_{j=1}^{L}>0$ be the _sequence of effective
relative perturbations_ and $\lambda>0$ be the _step size._ Moreover, let
$M\in\mathbb{N}$ and let
$(x_{i})_{i=1}^{M}\subset\Omega,(y_{i})_{i=1}^{M}\subset\mathbb{R}$ be the
training samples.

Let $j\in\\{1,\dots,L\\}$. The _exact gradient descent update of the biases in
the $j$-th layer_ is given by $u_{j}^{b}$, which is defined as

$\displaystyle u_{j}^{b}$
$\displaystyle\coloneqq\frac{1}{M}\sum_{i=1}^{M}u_{j,i}^{b},\text{ where}$
(2.2) $\displaystyle u_{j,i}^{b}$ $\displaystyle\coloneqq
I_{j}(x_{i})A_{j+1}^{T}I_{j+1}(x_{i})A_{j+2}^{T}\cdots
I_{L-1}(x_{i})A_{L}^{T}\ell^{\prime}(y_{i},\mathrm{R}(\Phi)(x_{i}))\,\text{
for }i\in\\{1,\ldots,M\\},$ (2.3)

where $I_{j}(x)\in\\{0,1\\}^{N_{j}}$ with $(I_{j}(x))_{k}=1$ if and only if
$\mathrm{R}(\Phi_{j})(x)\geq 0$. The _exact gradient descent update of the
weights in the $j$-th layer_ is defined as

$\displaystyle U_{j}^{w}$
$\displaystyle\coloneqq\frac{1}{M}\sum_{i=1}^{M}u_{j,i}^{b}\cdot\varrho({\mathrm{R}(\Phi_{j-1})({x_{i}})})^{T}.$

Let $\Theta_{j}^{b},\Theta_{j}^{w}\in\mathbb{R}^{N_{j}\times N_{j-1}}$ be two
independent random variables of diagonal matrices containing i.i.d entries
uniformly distributed on $[-0.5,0.5]$. The _perturbed bias and weight updates_
are defined as

$\displaystyle\hat{u}_{j}^{b}$
$\displaystyle\coloneqq(I+\varepsilon_{j}\Theta_{j}^{b})u_{j}^{b},$ (2.4)
$\displaystyle\hat{U}^{w}_{j}$
$\displaystyle\coloneqq(I+\varepsilon_{j}\Theta_{j}^{w})U_{j}^{w}.$

We define the _updated exact weight and bias matrices_ as

$\displaystyle\mathbb{A}_{j}\coloneqq A_{j}-\lambda
U_{j}^{w},\quad\mathbb{b}_{j}\coloneqq b_{j}-\lambda u_{j}^{b}$

and the _updated perturbed weight and bias matrices_ as

$\displaystyle\widehat{\mathbb{A}}_{j}\coloneqq
A_{j}-\lambda\hat{U}_{j}^{w},\quad\widehat{\mathbb{b}}_{j}\coloneqq
b_{j}-\lambda\hat{u}_{j}^{b}.$ (2.5)

The _update of $\Phi$ with sequence of effective relative perturbations
$\varepsilon$_ is the random variable

$\displaystyle\widehat{\Phi}^{\varepsilon}\coloneqq\big{(}(\widehat{\mathbb{A}}_{1},\widehat{\mathbb{b}}_{1}),(\widehat{\mathbb{A}}_{2},\widehat{\mathbb{b}}_{2}),\dots,(\widehat{\mathbb{A}}_{L},\widehat{\mathbb{b}}_{L})\big{)}.$

Moreover, for a domain $\Omega$ and for $j\in\\{1,\dots,L-1\\}$, we define the
_perturbed preactivations_ as
$\hat{\eta}_{j}\coloneqq\mathrm{R}(\widehat{\Phi}^{\varepsilon}_{j})\colon\Omega\to\mathbb{R}^{N_{j}}$.

######  Remark 2.4.

The perturbed updates of Definition 2.3 are a model for numerical errors
arising in floating-point arithmetic. In this model, the effective relative
perturbations $\varepsilon_{j}$ introduced in Definition 2.3 comprise all
numerical errors resulting from repeated applications of matrix
multiplications as well as from summing over potentially very large data sets.
The effective relative perturbations can therefore be assumed to be bounded
from below by the machine precision. However, as shown in Section 1.2,
unstable NNs may significantly amplify a perturbation. In addition to the
amplification, repeated matrix vector multiplications as well as the summation
over the data sets, will result in the accumulation of errors. We therefore
expect the perturbations at lower layers (corresponding to small $j$ in the
above notations) to be larger and those at final layers (corresponding to
large $j$) to be closer to the machine precision.

The full iteration of gradient descent with perturbed updates corresponds to
repeatedly applying Definition 2.3. We present a formal definition below.

######  Definition 2.5.

Let $d\in\mathbb{N},\Omega\subset\mathbb{R}^{d}$, let $\Phi$ be a NN of depth
$L\in\mathbb{N}$, input dimension $d$, and output dimension $1$. Further, let
$M\in\mathbb{N}$, and let
$X=(x_{i})_{i=1}^{M}\subset\Omega,Y=(y_{i})_{i=1}^{M}\subset\mathbb{R}$ be the
training samples. Let $\varepsilon$ be the sequence of effective relative
perturbations and $(\lambda_{n})_{n\in\mathbb{N}}\subset\mathbb{R}^{+}$ be the
sequence of step sizes.

We define a series of NNs by
$(\widehat{\Phi}^{\varepsilon,n})_{n\in\mathbb{N}}$, where for
$n\in\mathbb{N}$, $\widehat{\Phi}^{\varepsilon}_{n+1}$ is the update of
$\widehat{\Phi}^{\varepsilon,n}$ with effective relative perturbation
$\varepsilon$ and step size $\lambda_{n}$, and
$\widehat{\Phi}^{\varepsilon}_{1}=\Phi$. We call
$(\widehat{\Phi}^{\varepsilon,n})_{n\in\mathbb{N}}$ the _training sequence
with initialisation_ $\Phi$.

Next, we present two crucial assumptions on the perturbed gradient descent
updates of Definition 2.3. First, we assume that the average of the
reciprocals of the non-zero bias updates is not too large. In addition, we
assume that if the update of the bias weight of a neuron is equal to zero,
then the associated neuron is dead.

######  Assumption 2.6.

Let $d\in\mathbb{N}$ and $\Omega\subset\mathbb{R}^{d}$ be a domain. Let
$\nu,c_{0}>0$. Maintaining the notation of Definition 2.3, we assume that
$\sum_{j=1}^{L-1}\sum_{k=1}^{N_{j}}|(u_{j}^{b})_{k}|^{\dagger}\leq
N^{\nu}c_{0}$, where $x^{\dagger}=x^{-1}$ for $x>0$ and $0^{\dagger}=0$ and
$N=d+\sum_{j=1}^{L}N_{j}$.

In addition, we assume for all $j\in\\{1,\dots,L-1\\}$ and
$k\in\\{1,\dots,N_{j}\\}$ that $(u_{j}^{b})_{k}=0$ implies
$\hat{\eta}_{j}(x)\leq 0$ for all $x\in\Omega$.

The second assumption that we make is that the derivative of each of the
preactivated neurons is bounded by one.

######  Assumption 2.7.

Let $d\in\mathbb{N}$ and $\Omega\subset\mathbb{R}^{d}$ be a domain.
Maintaining the notation of Definition 2.3, we assume that for all
$j\in\\{1,\ldots,L-1\\}$ and $k\in\\{1,\ldots,N_{j}\\}$, we have that
$|(\nabla_{x}(\eta_{j})_{k})|\leq 1$ or all $x\in\Omega$.

We end this section by discussing the sensibility of the Assumptions 2.6 and
2.7 above. Assumption 2.6 requires the average of the reciprocals of the
updates of the bias vectors to bounded by $c_{0}N^{\nu-1}$. We numerically
check this assumption in Section 4 and find that for $\nu=2$ the first part of
the assumption is satisfied with probability at least $1/2$ for $c_{0}<0.1$.
We also show theoretically that under mild assumptions on the distribution of
$(u_{j}^{b})_{k}$, we can expect Assumption 2.6 to hold for moderate $\nu$.
Moreover, when sufficiently many training samples are used, any neuron such
that the derivative with respect to its bias value is zero at all points in
the training set is likely to be dead. Since this implies that the associated
bias will not be changed in the gradient descent step and since the bias value
is the most important (but not only) parameter to determine whether a neuron
lives, it is likely that the neuron remains dead after one gradient descent
step.

On the other hand, Assumption 2.7 requires the output of each neuron to have a
gradient of length not exceeding one. It will be clear from the proofs of our
main results that the bound of 1 can be replaced by any other positive number.
Note that the boundedness of the gradient is reasonable to assume in general
since a) this is a direct consequence of most initialisation schemes [14], b)
it is a sufficient condition to prevent the amplification effects that were
the basis of the examples found in Subsection 1.2 of the introduction, c)
quite often an additional regularisation, such as weight decay [8, Section
7.1], is used, which promotes smaller weights and implies bounds on the
derivatives of the outputs of each neuron.

##  3 The number of affine pieces generated in perturbed gradient descent

In this section, we demonstrate that, in the framework of Section 2 and under
Assumptions 2.6 and 2.7, a neural network trained via gradient descent with an
appropriately chosen step size will not admit a high number of affine pieces.
We demonstrate in Theorem 3.4 that the realisation of a NN after one step of
gradient descent, in expectation, admit a number of affine pieces that scales
polynomially with the number of neurons. Moreover, we show in Theorem 3.6,
that the number of affine pieces is polynomial in the number of gradient
descent steps and in the number of neurons of the neural network. The
polynomial dependence on the number of neurons in both results depends on
Assumption 2.6 and this polynomial is numerically identified to be quadratic
in Section 4.

Before we can state the results, we need some auxiliary results which will be
collected in Subsection 3.1. In Subsection 3.2, we combine these results to
obtain an upper bound on the expected number of affine pieces in a full neural
network. Finally, we present a high-probability upper bound on the number of
affine pieces during the full iteration of gradient descent in Subsection 3.3.

###  3.1 Expected number of generated affine pieces in one neuron

Intuitively speaking, a neuron of the form $x\mapsto\varrho(h(x)+b)$, where
$h$ is a piecewise affine map mapping from $\mathbb{R}^{d}$ to $\mathbb{R}$,
and where $b\in\mathbb{R}$, generates an affine piece, if $h(x)+b=0$ for a
point $x$ at which $h$ is affine and not constant. In all other cases,
$h(\cdot)+b$ is either constant or the range of $h(\cdot)+b$ locally lies in
one of the two linear regions of $\varrho$, which implies that the regularity
of $\varrho(h(x)+b)$ is the same as that of $h(x)$.

If $b$ in the argument above is chosen randomly, we expect that we can
quantify the probability of generating a given number of affine pieces. The
following lemma is a first step in this direction. To not disturb the flow of
reading, the proof of this auxiliary result was deferred to Appendix A.2.1.

######  Lemma 3.1.

Let $c>0$ and $h$ be a piecewise affine function on $[0,c]$ with
$P\in\mathbb{N}$ affine pieces. Let $t\in\mathbb{N}$ and $A$ be a Lebesgue
measurable set and assume that for every $y\in A$ it holds that

$\\#\\{x\in[0,1]\colon h(x)=y\\}\geq t.$

Then, $c\|h^{\prime}\|_{\infty}\geq\|h^{\prime}\|_{1}\geq\lambda(A)t$, where
$\lambda$ is the Lebesgue measure.

######  Remark 3.2.

Lemma 3.1 can be reformulated as saying that

$\mathbb{P}(\\#\\{x\in[0,c]\colon h(x)=U\\}\geq t)\leq c\delta/t,$

where $U$ is a uniformly distributed random variable on an interval of length
$\delta$ and $|h^{\prime}|\leq 1$. An interpretation of this is that the
probability of having many intersections of a uniform random variable with a
function of bounded gradient is bounded by the inverse width of the uniform
distribution and the number of intersections.

######  Remark.

Note that the estimate in Lemma 3.1 is independent from the number of affine
pieces $P$.

In the framework of Definition 2.3, we have that the updates at each gradient
descent step are randomly perturbed. It follows from Lemma 3.1, that a random
bias vector is unlikely to generate a high number of affine pieces. As a
consequence, we obtain the following result, which is proved in Appendix
A.2.2.

######  Proposition 3.3.

Let $L,d\in\mathbb{N}$, $N_{1},\dots N_{L-1}\in\mathbb{N}$, $N_{L}=1$, and let
$\Phi$ be a NN with architecture $(d,N_{1},\dots,N_{L})$. Let
$\widehat{\Phi}^{\varepsilon}$ be the NN after one backpropagation step with
sequence of effective relative perturbation
$\varepsilon=(\varepsilon_{1},\ldots,\varepsilon_{L})\in(0,\infty)^{L}$
satisfying Assumption 2.7 and step size $\lambda>0$. Let $\kappa$ be a line in
$[0,1]^{d}$. For every $j\in\\{1,\ldots,L-1\\}$ and
$k\in\\{1,\ldots,N_{j}\\}$, let $\omega_{j,k}\subset\kappa$ be a set that
contains all points along which
$\mathrm{R}(\widehat{\Phi}^{\varepsilon}_{j})_{k}$ restricted to $\kappa$ is
not locally affine. Furthermore, let $\hat{\omega}_{j,k}\subset\kappa$ be the
smallest set that contains all points along which
$\varrho(\mathrm{R}(\widehat{\Phi}^{\varepsilon}_{j})_{k})$ restricted to
$\kappa$ is not locally affine.

Then, we have for all $q\in\mathbb{N}$ that

$\displaystyle\mathbb{P}(\\#\left(\hat{\omega}_{j,k}\setminus\omega_{j,k}\right)\geq
q)\leq
2\mathcal{L}(\kappa)\cdot(\lambda\varepsilon_{j}q|(u_{j}^{b})_{k}|)^{-1}\,,$
(3.1)

where $\mathcal{L}(\kappa)$ is the length of $\kappa$. Moreover,

$\displaystyle\mathbb{E}(\\#\left(\hat{\omega}_{j,k}\setminus\omega_{j,k}\right))\leq\frac{2\mathcal{L}(\kappa)}{\lambda\,\varepsilon_{j}}\,|(u_{j}^{b})_{k}|^{-1}\,\ln(N^{j+1})\,.$
(3.2)

###  3.2 Expected number of affine pieces of neural networks after one
gradient descent step

Next, we bound from above the expected number of affine pieces of a full NN
after a single gradient descent update according to Definition 2.3. In
Proposition 3.3, we bounded from above the expected number of added affine
pieces by one neuron. The expected number of affine pieces added in the full
NN is found by adding all the added affine pieces corresponding to individual
neurons. The following theorem then follows by the linearity of the expected
value. We present a proof in Appendix A.3.

######  Theorem 3.4.

Let $L,d\in\mathbb{N}$, and $N_{1},\dots N_{L-1}\in\mathbb{N}$, $N_{L}=1$, and
let $\Phi$ be a NN with architecture $(d,N_{1},\dots,N_{L})$. Further, let
$\widehat{\Phi}^{\varepsilon}$ be a NN after one backpropagation step with
sequence of effective relative perturbation $\varepsilon\in(0,\infty)^{L}$
satisfying Assumption 2.6 with $c_{0},\nu>0$ and Assumption 2.7 and with step
size $\lambda>0$. Then, for every line $\kappa\subset[0,1]^{d}$ of length
$\mathcal{L}(\kappa)$, we have for every $j^{\prime}=1,\dots,L$

$\displaystyle\mathbb{E}(\mathrm{pieces}(\mathrm{R}(\widehat{\Phi}^{\varepsilon}),\kappa))\leq\left(1+\frac{2c_{0}}{\lambda}\,\mathcal{L}(\kappa)\,j^{\prime}\hat{\varepsilon}_{j^{\prime}}^{-1}N^{\nu}\ln(N)\right)\cdot(2N)^{L-j^{\prime}},$
(3.3)

where
$\hat{\varepsilon}_{j^{\prime}}\coloneqq\min_{j<j^{\prime}}\epsilon_{j}$.

######  Remark 3.5.

To interpret (3.3), we note that depending on the choice of $j^{\prime}$ the
upper bound varies between fully exponential in $L$, for $j^{\prime}=1$, and
linear in $L$, for $j^{\prime}=L$. Moreover, if $L-j^{\prime}$ is fixed, then
the whole estimate is a polynomial upper bound on the number of affine pieces
in terms of the number of neurons $N$.

On a theoretical level, this demonstrates that the asymptotic scaling of the
number of affine pieces of a NN after one backpropagation step does, in
expectation, not scale exponentially with the number of layers.

For practical purposes, it should be mentioned that this upper bound is void
if $\epsilon_{j}$ is too small for some $j<j^{\prime}$. However, we expect
that, for NNs with many neurons per layer, already for reasonably small
$j^{\prime}$, performing $j^{\prime}$ matrix-vector multiplications as in
(2.2) as well as summations over large data sets should result in the
accumulation of errors such that the effective perturbations $\epsilon_{j}$
with $j<j^{\prime}$ are sufficiently large.

Theorem 3.4 demonstrates that the number of affine pieces that are expected
after one gradient descent step is polynomial in the number of layers. This is
problematic for constructions of NNs that create an exponential number of
affine pieces with respect to the number of layers: _when the effective
perturbations are not too small_ , it is unlikely that a NN with a number of
affine pieces exponentially large with respect to the number of layers is
found after a single gradient descent step.

###  3.3 Generation of affine pieces during full iteration

We now take a closer look at the implications of Theorem 3.4 when optimising a
NN using gradient descent. As a consequence of Theorem 3.4 and Markov’s
inequality for a given line $\kappa$ it holds that

$\displaystyle\mathbb{P}\left(\mathrm{pieces}(\mathrm{R}(\widehat{\Phi}^{\varepsilon}),\kappa)\geq
2\min_{j^{\prime}\leq
L}\big{(}1+\frac{2c_{0}}{\lambda}\,\mathcal{L}(\kappa)\,j^{\prime}\hat{\varepsilon}_{j^{\prime}}^{-1}N^{\nu}\ln(N)\big{)}\cdot(2N)^{L-j^{\prime}}\right)\leq\frac{1}{2}.$
(3.4)

We continue by studying the impact of Equation 3.4 when training a NN with the
gradient descent algorithm and obtain the following main result.

######  Theorem 3.6.

Let $d\in\mathbb{N}$ and $\Omega\subset[0,1]^{d}$. Let $\Phi$ be a NN of depth
$L\in\mathbb{N}$ and input dimension $d$ and output dimension $1$,
$M\in\mathbb{N}$ and
$X=(x_{i})_{i=1}^{M}\subset\Omega,Y=(y_{i})_{i=1}^{M}\subset\mathbb{R}$ be the
training samples. Let $\varepsilon>0$ be the effective relative perturbation
and $(\lambda_{n})_{n\in\mathbb{N}}$ be the sequence of step sizes. Let
$(\widehat{\Phi}^{\varepsilon,n})_{n\in\mathbb{N}}$ be the associated
_training sequence with initialisation_ $\Phi$. Assume that
$\kappa\subset[0,1]^{d}$ is a line.

Then, for each $n\in\mathbb{N}$ satisfying Assumption 2.6 with
$c_{0}=c_{0}^{(n)}$ and Assumption 2.7 in the update step of
$\widehat{\Phi}^{\varepsilon,n}$, we have that

$\displaystyle\mathbb{P}\left(\mathrm{pieces}(\mathrm{R}(\widehat{\Phi}^{\varepsilon,n}),\kappa)\geq
2\min_{j^{\prime}\leq
L}\Big{(}1+\frac{2c_{0}^{(n)}}{\lambda_{n}}\,\mathcal{L}(\kappa)\,j^{\prime}\,\hat{\varepsilon}_{j^{\prime}}^{-1}\,N^{\nu}\ln
N\Big{)}\cdot(2N)^{L-j^{\prime}}\right)\leq\frac{1}{2}\,.$

######  Remark 3.7.

As presented in [33, 32], gradient descent as used in most machine learning
applications typically achieves a convergence order of $1/\sqrt{n}$ with a
corresponding step size of $1/\sqrt{n}$. Therefore, for a converging training
sequence $(\widehat{\Phi}^{\varepsilon,n})_{n\in\mathbb{N}}$, we can assume
that

$\displaystyle\widehat{\mathcal{R}}(R(\widehat{\Phi}^{\varepsilon,n}(X),Y))\gtrsim\frac{1}{\sqrt{n}}\text{
and }\widehat{\mathcal{R}}(R(\widehat{\Phi}^{\varepsilon,n}(X),Y))\to 0\text{
for }n\to\infty.$ (3.5)

Considering the square loss, the scaling of (3.5) shows that it is reasonable
to assume that

$\displaystyle\lambda_{n}\,\bigg{|}\frac{1}{M}\sum_{i=1}^{M}\ell^{\prime}(y_{i},\mathrm{R}(\widehat{\Phi}^{\varepsilon}_{n-1}))\bigg{|}\gtrsim
n^{-5/4}\,.$ (3.6)

Indeed, assuming the reverse inequality of (3.6) shows by a direct computation
that the weights of the NNs
$(\widehat{\Phi}^{\varepsilon,n})_{n\in\mathbb{N}}$ would converge with a rate
faster than $n^{-1/4}$. This, by the local Lipschitz property of the
realisation function [36, Proposition 4.1] yields a pointwise convergence of
$(\mathrm{R}(\widehat{\Phi}^{\varepsilon,n}))_{n\in\mathbb{N}}$ to a limit
with a rate faster than $n^{-1/4}$. Due to the quadratic dependence of the
square loss on the pointwise error, this would violate (3.5).

As a consequence of (3.6) as well as the linearity of (2.2) and (2.3), it is
sensible to claim that $c_{0}^{(n)}$ of Assumption 2.6 should scale as
$(n^{-5/4}/\lambda_{n})^{-1}=n^{3/4}$ with respect to $n\in\mathbb{N}$. For
other loss functions than the square loss, similar arguments can be made, that
may yield different exponents depending on the modulus of continuity of the
loss.

Based on (2.3), we, therefore, conclude that an admissible scaling law for
Assumption 2.6 for $n\in\mathbb{N}$ is

$\displaystyle\lambda_{n}^{-1}c_{0}^{(n)}\lesssim n^{5/4}.$

We observe that, if $\lambda_{n}c_{0}^{(n)}$ scales polynomially with respect
to $n\in\mathbb{N}$, then Theorem 3.6 states that NNs obtained by training
with gradient descent, with probability $1/2$, will not have a number of
affine pieces exponentially large with respect to the number of layers, unless
an exponential number of iterations is performed. Remark 3.7 shows that
$c_{0}^{(n)}$ can be assumed to grow not faster than $n^{5/4}$ with respect to
$n\in\mathbb{N}$.

Since, by construction, the NNs
$(\widehat{\Phi}^{\varepsilon,n})_{n\in\mathbb{N}}$ are independent random
variables, it holds that _the event that the sequence
$(\mathrm{R}(\widehat{\Phi}^{\varepsilon,n}))_{n\in\mathbb{N}}$ maintains many
affine pieces for a long series of gradient descent steps is highly unlikely_.

##  4 Experimental study of the assumptions

$N^{-2}\,\sum_{j=1}^{L-1}\sum_{k=1}^{N_{\ell}}|(u_{j}^{b})_{k}|^{\dagger}$

Figure 1: Cumulative distribution function for results of the experiment to
estimate the value of (4.1) in practice. In about half of the cases the value
(4.1) is bounded by $0.1$.

Figure 2: Relative frequencies of the bias updates $(u_{j}^{b})_{k}$ computed
over 10’000 neural network instantiations with $N=200$ neurons per layer and,
from top to bottom, $L=4,5,6$ layers. We only include updates with absolute
value less than one since we are interested in the asymptotic behaviour at
zero. The constants $C_{1},C_{2},C_{3}$ in the plots are $1/160,1/270,1/160$,
respectively. The right panel shows the same distribution in a log-log plot
and indicates that the density of the bias updates has a mild algebraic pole,
of order $1/2$, at the origin since the slope of the distribution appears to
be not steeper than $-1/2$.

To verify the sensibility of Assumption 2.6, we numerically check if the
gradient descent update of the biases of a randomly initialised NN is on
average appropriately far away from zero. For $L=3,4,5,6,7$, we construct
random NNs with $L$ layers and input and output dimensions equal to $1$ and
with various numbers of neurons between $100$ and $400$ per layer. We
initialise the NNs according to the He initialiser [16]. Next, we choose $500$
random training points $(x_{i})_{i=1}^{500}$ uniformly in $[0,2\pi]$ and
associate as labels $(\sin(x_{i}))_{i=1}^{500}$. We perform one step of
gradient descent and collect the size of the updates of the bias vectors.
Finally, we compute for each NN the value

$\displaystyle\frac{1}{N^{2}}\sum_{j=1}^{L-1}\sum_{k=1}^{N_{\ell}}|(u_{j}^{b})_{k}|^{\dagger},$
(4.1)

where $N\in\mathbb{N}$ is the total number of neurons of the NN and
$(u_{j}^{b})_{k}$ are the updates of the bias vectors as in Assumption 2.6.
The term (4.1) corresponds to $\nu=2$ in Assumption 2.6.

In Figure 1, we depict for each number of layers, the resulting values of
(4.1). More precisely, for number of layers $L$ and for each number
$N_{1}\in\\{100,110,\dots,390,400\\}$, we initialise $500$ neural networks
with architecture $(1,N_{1},\dots,N_{1},1)$ and compute the associated values
of (4.1). We average these values over the 500 runs. Then, we depict in Figure
1, the maximum value of these averages for all values $N_{1}$. We observe
that, in at least half of the cases the values (4.1), are bounded by
$c_{0}<0.1$.

Theoretically, we can also justify Assumption 2.6, if we assume that all bias
updates are independent identically distributed random variables with
compactly supported densities $\sigma$ with
$\sigma(x)\in\mathcal{O}(x^{-(1-\alpha)})$ for $x\to 0$, where
$\alpha\in(0,1)$.

In this case, we have that for a $q>1/\alpha$ with $q\in\mathbb{N}$,

$\displaystyle\left(\frac{1}{N^{q}}\sum_{j=1}^{L-1}\sum_{k=1}^{N_{\ell}}|(u_{j}^{b})_{k}|^{\dagger}\right)^{1/q}\leq\frac{1}{N}\sum_{j=1}^{L-1}\sum_{k=1}^{N_{\ell}}\left(|(u_{j}^{b})_{k}|^{\dagger}\right)^{1/q}$

by the binomial theorem. Moreover, for large $N\in\mathbb{N}$, we will have by
the law of large numbers that with high probability

$\displaystyle\frac{1}{N}\sum_{j=1}^{L-1}\sum_{k=1}^{N_{\ell}}\left(|(u_{j}^{b})_{k}|^{\dagger}\right)^{1/q}\approx\int_{(0,\infty)}x^{-1/q}\sigma(x)dx\lesssim\int_{(0,c)}x^{-1/q+\alpha-1}dx<\infty,$

where $c$ is such that $\operatorname{supp\,}\sigma\subset[0,c]$ and we used
that $-1/q+\alpha>0$.

We can study the distribution of $(u_{j}^{b})_{k}$ numerically to see if it,
indeed, has a mild pole at zero. We use the same setup as in the experiment of
Figure 1 and observe in Figure 2 that the distribution $\sigma$ of
$(u_{j}^{b})_{k}$, with some $\delta>0$, appears to satisfy $\sigma(x)\simeq
x^{-1/2+\delta}$ as $x\to 0$. By the previous argument, this implies that
Assumption 2.6 should hold with $\nu=2$. This estimate coincides with the
findings of the experiment that was shown in Figure 1.

##  5 Experimental study of the main theorem

In this section, we study the number of pieces of NN realisations in practical
scenarios as well as the effect of the accuracy in floating-point
computations. Computing the number of pieces of a NN by identifying regions
with gradients differing by a threshhold is extremely sensitive to the choice
of the threshhold, and therefore does not yield reliable results. Instead, we
will count so-called activation regions, which were considered as a proxy for
affine pieces in [14]. Indeed, the number of activation regions is an upper
bound on the number of affine pieces. We refer to [14, Definition 1] for the
definition of activation regions and [14, Lemma 3] for the connection between
activation regions and affine pieces.

In our first experiment, depicted in Figures 3 and 4, we train NNs with
various numbers of layers and neurons to approximate the functions $x\mapsto
1-x^{2}/2$ and $x\mapsto\cos(x)$, respectively. For that, the mean square loss
on $500$ randomly chosen training points is minimised via gradient descent. As
a step size, we use $\lambda_{n}=0.02/(1+\sqrt{n}/8)$ at iteration $n$. This
step size rule was empirically found to yield the fastest convergence. The NNs
are randomly initialised according to the He initialisation [16]. In the
computation of the gradient as well as when evaluating the NNs, we relatively
perturb the output of each matrix-vector multiplication with relative noise
the amplitude of which will be specified below.

Figure 3: Approximation of the function $x\mapsto 1-x^{2}/2$ on $[0,1]$ by
neural networks with $50$ neurons per layer and $5$ layers (top row) and $50$
neurons per layer and $5$ layers (bottom row). Each matrix-vector
multiplication during training is perturbed by relative noise of size
$\\{10^{-5},5\cdot 10^{-5},10^{-4},5\cdot 10^{-4},10^{-3}\\}$. All results are
averaged over five runs.

Figure 4: Approximation of the function $x\mapsto\cos(x)$ on $[0,1]$ by neural
networks with $50$ neurons per layer and $5$ layers (top row) and $50$ neurons
per layer and $8$ layers (bottom row). Each matrix-vector multiplication
during training is perturbed by relative noise of size $\\{10^{-5},5\cdot
10^{-5},10^{-4},5\cdot 10^{-4},10^{-3}\\}$. All results are averaged over five
runs.

In Figures 3 and 4, we observe that the number of activation regions _always
stays far below the upper bound of Theorem 3.6_. However, there is no clear
effect of the accuracy in the computations. Note that, in the regime far away
from the upper bound of Theorem 3.6, it is conceivable that noise in the
updates of the biases can even create new affine pieces, since, in this case,
as seen in Remark 3.2, the generation of affine pieces is not necessarily a
low probability event. Hence, it should not be surprising that no clear effect
of the level of numerical accuracy is visible in this regime. Moreover, we
expect that in addition to the accuracy of the computations, there are further
reasons which prevent the creation of many pieces in practice. To make the
effect of the accuracy in floating-point arithmetic visible, we need to enter
a regime, where the underlying NNs already have a lot of affine pieces and
then observe the effect of the accuracy on the number of pieces. We propose
the following example: As famously shown in [46, Proposition 2], the function
$x\mapsto x^{2}$ can be very efficiently approximated by ReLU NNs. The NN with
$L$ layers and four neurons per layer given by

$\displaystyle b_{\ell}$ $\displaystyle\coloneqq
4^{1-\ell}\left(0,0,-1,-2\right)^{T}\text{ for }\ell=1,\dots,L-1,\qquad
b_{L}\coloneqq 0,$ (5.1) $\displaystyle A_{1}$
$\displaystyle\coloneqq\left(1,2,2,2\right)^{T},\qquad
A_{\ell}\coloneqq\left(\begin{array}[]{cccc}1&-1&2&-1\\\ 0&-1/2&1&-1/2\\\
0&-1/2&1&-1/2\\\ 0&-1/2&1&-1/2\end{array}\right)\text{ for }\ell=2,\dots,L-1,$
(5.6) $\displaystyle A_{L}$ $\displaystyle\coloneqq(1,-1,2,-1),$

satisfies $|\mathrm{R}(\Phi)(x)-x^{2}|\leq 4^{-L}$ for all $x\in[0,1]$ and has
$2^{L-1}$ affine pieces.

We initialise a NN in the way above, and train it with step size
$\lambda_{n}=0.1/(1+\sqrt{n})$, for $n=1,\dots,200$, on training data which is
comprised of 5’000 uniformly randomly sampled elements of $[0,1]$ with labels
described by the function $x\mapsto cx^{2}$, where $c=0.99999$. Note that,
simply changing $A_{L}$ to $\widetilde{A}_{L}=c\cdot(1,-1,2,-1)$ would already
achieve zero training error. Hence, the target function is very close to the
initialisation. In this example, we _clearly observe the effect of the
numerical accuracy_ , with low accuracy computations maintaining activation
regions much better and in some cases even increasing the number of activation
regions further. The results are depicted in Figure 5.

Figure 5: Training of neural networks initialised with many pieces according
to (5.1). In the top row 12 layers are used, in the bottom row 14 layers are
used. Each matrix-vector multiplication during training is perturbed by
relative noise of size $\\{5\cdot 10^{-5},10^{-4},5\cdot
10^{-4},10^{-3},5\cdot 10^{-3}\\}$. All results are averaged over five runs.

We observe in Figure 5 that if higher accuracy computations are used, then
fewer activation regions are lost during the iteration compared to iterations
with low accuracies. In all iterations, the mean squared error does _not_
significantly change during the iteration.

## References

  * [1] S. Arridge, P. Maass, O. Öktem, and C.-B. Schönlieb.  Solving inverse problems using data-driven models.  Acta Numerica, 28:1–174, 2019.
  * [2] M. Bachmayr and V. Kazeev.  Stability and preconditioning of elliptic PDEs with low-rank multilevel structure.  Foundations of Computational Mathematics, 20:1175–1236, 2020.
  * [3] K. Bhattacharya, B. Hosseini, N. B. Kovachki, and A. M. Stuart.  Model reduction and neural networks for parametric PDEs.  arXiv preprint arXiv:2005.03180, 2020.
  * [4] H. Boche, A. Fono, and G. Kutyniok.  Limitations of deep learning for inverse problems on digital hardware.  arXiv preprint arXiv:2202.13490, 2022.
  * [5] H. Bölcskei, P. Grohs, G. Kutyniok, and P. C. Petersen.  Optimal approximation with sparsely connected deep neural networks.  SIAM Journal on Mathematics of Data Science, 1:8–45, 2019.
  * [6] M. J. Colbrook, V. Antun, and A. C. Hansen.  The difficulty of computing stable and accurate neural networks: On the barriers of deep learning and Smale’s 18th problem.  Proceedings of the National Academy of Sciences, 119(12):e2107151119, 2022.
  * [7] W. E and B. Yu.  The deep Ritz method: a deep learning-based numerical algorithm for solving variational problems.  Communications in Mathematics and Statistics, 6(1):1–12, 2018.
  * [8] I. Goodfellow, Y. Bengio, and A. Courville.  Deep Learning.  MIT Press, 2016.
  * [9] P. Grohs, F. Hornung, A. Jentzen, and P. Von Wurstemberger.  A proof that artificial neural networks overcome the curse of dimensionality in the numerical approximation of Black-Scholes partial differential equations.  Memoirs of the American Mathematical Society, 2020.
  * [10] S. Gupta, A. Agrawal, K. Gopalakrishnan, and P. Narayanan.  Deep learning with limited numerical precision.  In International conference on machine learning, pages 1737–1746. PMLR, 2015.
  * [11] I. Gühring, G. Kutyniok, and P. C. Petersen.  Error bounds for approximations with deep ReLU neural networks in ${W}^{s,p}$\- norms.  arXiv preprint arXiv:1902.07896, 2019.
  * [12] J. Han, A. Jentzen, and W. E.  Solving high-dimensional partial differential equations using deep learning.  Proceedings of the National Academy of Sciences, 115(34):8505–8510, 2018.
  * [13] S. Han, H. Mao, and W. J. Dally.  Deep compression: Compressing deep neural networks with pruning, trained quantization and huffman coding.  arXiv preprint arXiv:1510.00149, 2015.
  * [14] B. Hanin and D. Rolnick.  Complexity of linear regions in deep networks.  In International Conference on Machine Learning, pages 2596–2604. PMLR, 2019.
  * [15] J. He, L. Li, J. Xu, and C. Zheng.  ReLU deep Neural Networks and Linear Finite Elements.  arXiv preprint arXiv:1807.03973, 2018.
  * [16] K. He, X. Zhang, S. Ren, and J. Sun.  Delving deep into rectifiers: Surpassing human-level performance on imagenet classification.  In Proceedings of the IEEE international conference on computer vision, pages 1026–1034, 2015.
  * [17] N. J. Higham.  Accuracy and Stability of Numerical Algorithms.  Society for Industrial and Applied Mathematics, second edition, 2002.
  * [18] I. Hubara, M. Courbariaux, D. Soudry, R. El-Yaniv, and Y. Bengio.  Quantized neural networks: Training neural networks with low precision weights and activations.  The Journal of Machine Learning Research, 18(1):6869–6898, 2017\.
  * [19] B. Jacob, S. Kligys, B. Chen, M. Zhu, M. Tang, A. Howard, H. Adam, and D. Kalenichenko.  Quantization and training of neural networks for efficient integer-arithmetic-only inference.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 2704–2713, 2018.
  * [20] K. H. Jin, M. T. McCann, E. Froustey, and M. Unser.  Deep convolutional neural network for inverse problems in imaging.  IEEE Transactions on Image Processing, 26(9):4509–4522, 2017.
  * [21] V. Kazeev, I. Oseledets, M. Rakhuba, and C. Schwab.  QTT-finite-element approximation for multiscale problems I: model problems in one dimension.  Advances in Computational Mathematics, 43(2):411–442, 2017.
  * [22] V. Kazeev, I. Oseledets, M. V. Rakhuba, and C. Schwab.  Quantized tensor FEM for multiscale problems: Diffusion problems in two and three dimensions.  Multiscale Modeling & Simulation, 20(3):893–935, 2022.
  * [23] V. Kazeev and C. Schwab.  Quantized tensor-structured finite elements for second-order elliptic PDEs in two dimensions.  Numerische Mathematik, 138:133–190, January 2018.
  * [24] A. Krizhevsky, I. Sutskever, and G. Hinton.  Imagenet classification with deep convolutional neural networks.  In Adv. Neural Inf. Process. Syst. 25, pages 1097–1105. Curran Associates, Inc., 2012.
  * [25] G. Kutyniok, P. Petersen, M. Raslan, and R. Schneider.  A theoretical analysis of deep neural networks and parametric pdes.  Constructive Approximation, 55(1):73–125, 2022.
  * [26] Y. LeCun, Y. Bengio, and G. Hinton.  Deep learning.  Nature, 521(7553):436–444, 2015.
  * [27] H. Li, J. Schwab, S. Antholzer, and M. Haltmeier.  Nett: Solving inverse problems with deep neural networks.  Inverse Problems, 36(6):065005, 2020.
  * [28] Z. Li, Y. Ma, C. Vajiac, and Y. Zhang.  Exploration of numerical precision in deep neural networks.  arXiv preprint arXiv:1805.01078, 2018.
  * [29] L. Lu, P. Jin, G. Pang, Z. Zhang, and G. E. Karniadakis.  Learning nonlinear operators via deeponet based on the universal approximation theorem of operators.  Nature Machine Intelligence, 3(3):218–229, 2021.
  * [30] C. Marcati, J. A. Opschoor, P. C. Petersen, and C. Schwab.  Exponential ReLU neural network approximation rates for point and edge singularities.  arXiv preprint arXiv:2010.12217, 2020.
  * [31] C. Marcati, M. Rakhuba, and C. Schwab.  Tensor rank bounds for point singularities in $\mathbb{R}^{3}$.  Advances in Computational Mathematics, 48(3):17–57, 2022.
  * [32] A. Nemirovski, A. Juditsky, G. Lan, and A. Shapiro.  Robust stochastic approximation approach to stochastic programming.  SIAM Journal on Optimization, 19(4):1574–1609, 2009.
  * [33] A. S. Nemirovsky and D. B. Yudin.  Problem Complexity and Method Efficiency in Optimization.  Wiley-Interscience Series in Discrete Mathematics. Wiley, 1983.
  * [34] G. Ongie, A. Jalal, C. A. Metzler, R. G. Baraniuk, A. G. Dimakis, and R. Willett.  Deep learning techniques for inverse problems in imaging.  IEEE Journal on Selected Areas in Information Theory, 1(1):39–56, 2020.
  * [35] J. Opschoor, P. Petersen, and C. Schwab.  Deep ReLU networks and high-order finite element methods.  Analysis and Applications, 18, 12 2019.
  * [36] P. Petersen, M. Raslan, and F. Voigtlaender.  Topological properties of the set of functions generated by neural networks of fixed size.  Foundations of computational mathematics, 21(2):375–444, 2021.
  * [37] P. Petersen and F. Voigtlaender.  Optimal approximation of piecewise smooth functions using deep ReLU neural networks.  Neural Networks, 108:296–330, 2018.
  * [38] A. F. L. Pineda and P. C. Petersen.  Deep neural networks can stably solve high-dimensional, noisy, non-linear inverse problems.  arXiv preprint arXiv:2206.00934, 2022.
  * [39] M. Raissi, P. Perdikaris, and G. E. Karniadakis.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  Journal of Computational Physics, 378:686–707, 2019.
  * [40] C. Schwab and J. Zech.  Deep learning in high dimension: Neural network expression rates for generalized polynomial chaos expansions in UQ.  Analysis and Applications, 17(01):19–55, 2019.
  * [41] A. W. Senior, R. Evans, J. Jumper, J. Kirkpatrick, L. Sifre, T. Green, C. Qin, A. Žídek, A. W. Nelson, and A. Bridgland.  Improved protein structure prediction using potentials from deep learning.  Nature, 577(7792):706–710, 2020.
  * [42] U. Shaham, A. Cloninger, and R. R. Coifman.  Provable approximation properties for deep neural networks.  Applied and Computational Harmonic Analysis, 44(3):537–557, 2018\.
  * [43] D. Silver, A. Huang, C. J. Maddison, A. Guez, L. Sifre, G. Van Den Driessche, J. Schrittwieser, I. Antonoglou, V. Panneershelvam, and M. Lanctot.  Mastering the game of go with deep neural networks and tree search.  Nature, 529(7587):484–489, 2016.
  * [44] Y. Sun, D. Lao, G. Sundaramoorthi, and A. Yezzi.  Surprising instabilities in training deep networks and a theoretical analysis.  arXiv preprint arXiv:2206.02001, 2022.
  * [45] M. Telgarsky.  Representation benefits of deep feedforward networks.  arXiv preprint arXiv:1509.08101, 2015.
  * [46] D. Yarotsky.  Error bounds for approximations with deep ReLU networks.  Neural Networks, 94:103–114, 2017.
  * [47] D. Yarotsky.  Optimal approximation of continuous functions by very deep ReLU networks.  In Conference on learning theory, pages 639–649. PMLR, 2018.
  * [48] S.-C. Zhou, Y.-Z. Wang, H. Wen, Q.-Y. He, and Y.-H. Zou.  Balanced quantization: An effective and efficient approach to quantized neural networks.  Journal of Computer Science and Technology, 32(4):667–682, 2017\.

##  Appendix A Appendix

###  A.1 Example of unstable neural network

Here, we introduce the notion of a _finite precision realisation_ of a given
neural network $\Phi$ and give an example of how numerical instability can
affect such NNs. We assume the floating-point addition and multiplication to
be defined according to (1.2), where the the minima are well defined since
$\mathbb{M}$ has no accumulation points. Further, we consider the floating-
point realisation of matrix-vector multiplication that is given by

$(A\times_{\mathbb{M}}b)_{j}\coloneqq(A_{j,1}\times_{\mathbb{M}}b_{1})+_{\mathbb{M}}(A_{j,2}\times_{\mathbb{M}}b_{2})+_{\mathbb{M}}\dots+_{\mathbb{M}}(A_{j,n}\times_{\mathbb{M}}b_{n}),$

for any matrix $A\in\mathbb{R}^{m\times n}$ and a vector $b\in\mathbb{R}^{n}$
with $m,n\in\mathbb{N}$, where $j=1,\dots m$ and the summation is performed
from left to right. Next, we define the finite precision realisation of a
neural network.

######  Definition A.1.

Let $\mathbb{M}\subset\mathbb{R}$ be a set without accumulation points. For a
neural network $\Phi$ with input dimension $d\in\mathbb{N}$, $L\in\mathbb{N}$
layers, output dimension $N_{L}\in\mathbb{N}$ and an activation function
$\varrho_{\mathbb{M}}\colon\mathbb{M}\to\mathbb{M}$, we define the _finite
precision realisation_ by

$\displaystyle\mathrm{R}_{\mathbb{M}}(\Phi)\colon{\mathbb{R}^{d}}\to\mathbb{M}^{N_{L}},x\mapsto
x^{(L)},$

where $x^{(L)}$ is defined as

$\begin{split}x^{(0)}&\coloneqq x,\\\
x^{(j)}&\coloneqq\varrho_{\mathbb{M}}(A_{\ell}\times_{\mathbb{M}}x^{(j-1)}+_{\mathbb{M}}b_{\ell})\quad\text{
for }j=1,\dots,L-1,\\\ x^{(L)}&\coloneqq
A_{L}\times_{\mathbb{M}}x^{(L-1)}+_{\mathbb{M}}b_{L}.\end{split}$ (A.1)

######  Remark A.2.

Observe that in this section we only require the operations of the NNs to be
elements of $\mathbb{M}$ but not the weights. This greatly simplifies the
mathematical argument while still allowing NNs to exhibit instabilities.

Having defined finite precision operations we continue by presenting a neural
network $\Phi$ such that for every non-negative input $x$ it holds that

$|\mathrm{R}(\Phi)(x)-\mathrm{R}_{\mathbb{M}}(\Phi)(x)|=|\mathrm{R}(\Phi)(x)|=|x|.$

This implies that the relative error of passing to the finite precision
network is 1.

To illustrate numerical instability, we introduce the following neural network
$\Phi$: Let $\lambda\geq 1$, $N,L\in\mathbb{N}$ and

$\displaystyle A_{1}=\left(\begin{array}[]{c}1\\\ \lambda\\\ \vdots\\\
\lambda\end{array}\right),\qquad
A_{\ell}=\left(\begin{array}[]{cccc}1&0&\dots&0\\\ 0&\lambda&\dots&\lambda\\\
0&\vdots&\ddots&\vdots\\\
0&\lambda&\dots&\lambda\end{array}\right),\qquad\text{ for }\ell=2,\dots L-3$

and

$\displaystyle A_{L-2}=\left(\begin{array}[]{cccc}1&0&\dots&0\\\
0&\lambda&\dots&\lambda\end{array}\right),\qquad
A_{L-1}=\left(\begin{array}[]{cc}1&\lambda\\\ 0&\lambda\\\
\end{array}\right),\qquad
A_{L}=\left(\begin{array}[]{cc}1&-1\end{array}\right).$

Furthermore, all biases satisfy $b_{\ell}=0$ for $\ell=1,\dots,L$. We call the
neural network defined above $\Phi_{\lambda,N,L}$, an illustration is provided
below in Figure 6.

Figure 6: The unstable neural network $\Phi_{\lambda,N,L}$. As defined
earlier, the weights between the green dots are always $\lambda$.

######  Proposition A.3.

Let $\lambda=10^{r}\geq 1$ for $r\in\mathbb{Z}$ and let $\mathbb{M}$ be as in
(1.1), with $b_{\min}\leq 0,b_{\max}>16+r$, $b=10$, $p=16$ and hence
$\epsilon=5\cdot 10^{-16}$. Furthermore, let $\varrho(x)=\max\\{0,x\\}$ and
$\varrho_{\mathbb{M}}(z)\coloneqq\varrho(z)$ for all $z\in\mathbb{M}$.
Additionally, let $L\in\mathbb{N}$ and $N\in\mathbb{N}$ be such that
$(L-3)\log_{10}(N-1)+(L-1)\log_{10}(\lambda-2\epsilon)\geq 16$, then the
neural network $\Phi_{\lambda,N,L}$ as defined above satisfies

$|\mathrm{R}(\Phi_{\lambda,N,L})(x)-\mathrm{R}_{\mathbb{M}}(\Phi_{\lambda,N,L})(x)|=|\mathrm{R}(\Phi_{\lambda,N,L})(x)|=|x|$

for all $x\in\mathbb{M}$ such that $x\geq 0$.

###### Proof.

It follows from the definition of the NN $\Phi_{\lambda,N,L}$ that
$\mathrm{R}(\Phi)(x)=x$. Thus, to prove the theorem we have to show that

$\displaystyle\mathrm{R}_{\mathbb{M}}(\Phi_{\lambda,N,L})(x)=0$ (A.2)

for all nonnegative $x$. Let, for $x\geq 0$, $x^{(0)},\dots,x^{(L)}$ be as in
(A.1). The construction of $\Phi_{\lambda,N,L}$ implies that $x^{(L)}=0$ if
$(x^{(L-1)})_{1}=(x^{(L-1)})_{2}$. Additionally,
$(x^{(L-1)})_{1}=(x^{(L-1)})_{2}$ follows if
$(x^{(L-2)})_{1}<\lambda(x^{(L-2)})_{2}\epsilon/2$.

As a consequence, we have to find a lower bound of $(x^{(L-2)})_{2}$. To begin
with, we observe that for $x\in\mathbb{M}$ the multiplication with $\lambda$
is exact and thus $(x^{(2)})_{k}=\lambda x$ for all $k=2,\dots,N$. Also, we
observe that for $j=2,\dots,L-3$ it holds for $k\in\\{2,\dots,n\\}$ that

$\displaystyle(N-1)\lambda\max_{i\in\\{2,\ldots,N-1\\}}(x^{(j-1)})_{i}+\epsilon(N-1)\lambda\max_{i\in\\{2,\ldots,N-1\\}}(x^{(j-1)})_{i}$
$\displaystyle\qquad\geq(A_{j}x^{(j-1)})_{k}$
$\displaystyle\qquad\geq(N-1)\lambda\min_{i\in\\{2,\ldots,N-1\\}}(x^{(j-1)})_{i}-\epsilon(N-1)\lambda\max_{i\in\\{2,\ldots,N-1\\}}(x^{(j-1)})_{i}.$

This inequality results from the assumption that in the summation of each row
of $A_{j}$ we create $(N-2)$ times an error at most of size
$\epsilon\cdot\max_{i\in\\{2,\ldots,N-1\\}}(x^{(j-1)})_{i}$. This leads to

$(x^{(L-3)})_{k}\geq\lambda(N-1)^{L-4}(\lambda-2\epsilon)^{L-4}x$

and

$(x^{(L-2)})_{2}\geq\lambda(N-1)^{L-3}(\lambda-2\epsilon)^{L-3}x$

for all $k\in\\{2,\dots,N\\}$. Due to the construction of our NN, we also have
that $(x^{(L-2)})_{1}=x$.

Thus, the inequality
$(x^{(L-2)})_{1}<\lambda(x^{(L-2)})_{2}\,\frac{\epsilon}{2}$ holds if

$\frac{1}{2}\lambda^{2}(N-1)^{L-3}(\lambda-2\epsilon)^{L-3}>\epsilon^{-1}.$

Note that

$\displaystyle\frac{1}{2}\lambda^{2}(N-1)^{L-3}(\lambda-2\epsilon)^{L-3}>\frac{1}{2}(N-1)^{L-3}(\lambda-2\epsilon)^{L-1}$
(A.3)

and hence, it suffices to show that the right-hand side of (A.3) is bounded
below by $\epsilon^{-1}=210^{15}$. Taking logarithms, implies that this
inequality is implied by

$(L-3)\log_{10}(N-1)+(L-1)\log_{10}(\lambda-2\epsilon)-\log_{10}(2)>\log_{10}(2)+15.$

Observing that $2\log_{10}(2)<1$ completes the proof. ∎

###  A.2 Auxiliary results

####  A.2.1 Proof of Lemma 3.1

###### Proof.

We denote the affine pieces of $h$ on $[0,c]$ by $(I_{p})_{p=1}^{P}$. Here we
choose the pieces as semi-open intervals $I_{p}=(r_{p},s_{p}]$ so that they
are disjoint. Moreover, we choose $x_{p}\in I_{p}$ for $p=1,\dots,P$. We claim
that

$\displaystyle\sum_{p=1}^{P}\lambda(h(I_{p}))\geq t\lambda(A).$ (A.4)

To prove this, we consider the auxiliary function

$\displaystyle\tilde{h}\colon[0,1]$ $\displaystyle\to
A\times\left(\\{1,\dots,P\\}\cup\\{\infty\\}\right),$
$\displaystyle\tilde{h}(x)$
$\displaystyle\coloneqq\left(h(x),\\#\\{z\in[0,x]\colon h(z)=h(x)\\}\right).$

Note that for $x\in[0,1]$ either $\\#\\{z\in[0,x]\colon h(z)=h(x)\\}\leq P$ or
$\\#\\{z\in[0,x]\colon h(z)=h(x)\\}=\infty$, and hence $\tilde{h}$ is well
defined. We have by assumption that

$\tilde{h}([0,1])\supset A\times\\{1,\dots,t\\}.$

If we equip $A\times(\\{1,\dots,P\\}\cup\\{\infty\\})$ with the measure
$\lambda\times\nu$, where $\lambda$ is the Lebesgue and $\nu$ is the counting
measure, then we conclude by the monotonicity of measures that

$\left(\lambda\times\nu\right)\left(\tilde{h}([0,1])\right)\geq\lambda(A)t.$

It is easy to see that

$\displaystyle\lambda(h(I_{p}))\leq(\lambda\times\nu)\left(\tilde{h}(I_{p})\right)\leq\sup_{\begin{subarray}{c}A_{1},\dots,A_{P+1}\subset
I_{p},\\\ \text{ measurable}\end{subarray}}\sum_{q=1}^{P+1}\lambda(h(A_{q})),$
(A.5)

If $h^{\prime}(x_{p})\neq 0$ on $I_{p}$, then the right-hand side of (A.5)
equals $\lambda(h(I_{p}))$. If $h^{\prime}(x_{p})=0$, then the right-hand side
of (A.5) is $0$ and by the non-negativity of the Lebesgue measure the left
hand-side is $0$ as well. We conclude that

$\displaystyle\lambda(h(I_{p}))=(\lambda\times\nu)(\tilde{h}(I_{p})).$ (A.6)

By the subadditivity of measures and (A.6), we have that

$\lambda(A)t\leq\sum_{p=1}^{P}(\lambda\times\nu)\left(\tilde{h}(I_{p})\right)=\sum_{p=1}^{P}\lambda(h(I_{p})),$

which yields (A.4). Finally, by standard properties of the Lebesgue measure
under Lipschitz maps, we have that

$\sum_{p=1}^{P}\lambda(h(I_{p}))\leq\sum_{p=1}^{P}h^{\prime}(x_{p})\lambda(I_{p})=\|h^{\prime}\|_{1}\leq
c\|h^{\prime}\|_{\infty},$

where we used that $I_{p}$ are disjoint. ∎

####  A.2.2 Proof of Proposition 3.3

###### Proof.

Part 1: (Proof of (3.1)) Let $\kappa\subset[0,1]^{d}$ be a straight line.
Recalling equation (2.5), we observe that

$\displaystyle\widehat{\mathbb{A}}_{j}\mathbbm{1}_{\hat{\eta}_{j-1}\geq
0}\,\hat{\eta}_{j-1}(x)+\widehat{\mathbb{b}}_{j}$
$\displaystyle=\widehat{\mathbb{A}}_{j}\mathbbm{1}_{\hat{\eta}_{j-1}\geq
0}\,\hat{\eta}_{j-1}(x)+(b_{j}-\lambda\hat{u}_{j}^{b})\,.$

Inserting $\hat{u}_{j}^{b}=(I+\varepsilon_{j}\Theta_{j}^{b})u_{j}^{b}$ from
equation (2.4) leads to

$\displaystyle\widehat{\mathbb{A}}_{j}\mathbbm{1}_{\hat{\eta}_{j-1}\geq
0}\,\hat{\eta}_{j-1}(x)+\widehat{\mathbb{b}}_{j}$
$\displaystyle=\widehat{\mathbb{A}}_{j}\mathbbm{1}_{\hat{\eta}_{j-1}\geq
0}\,\hat{\eta}_{j-1}(x)+(\mathbb{b}_{j}-\lambda\varepsilon_{j}\Theta_{j}u_{j}^{b})$
$\displaystyle\eqqcolon\widehat{\mathbb{A}}_{j}\mathbbm{1}_{\hat{\eta}_{j-1}\geq
0}\,\hat{\eta}_{j-1}(x)+\mathbb{b}_{j}+\theta_{j}.$

Let us introduce for $k\in\\{1,\ldots,N_{j}\\}$ the random variable

$Q_{\widehat{\mathbb{A}}_{j},\hat{\eta}_{j-1},\theta_{j},k}:=\\#\\{x\in\kappa\colon\big{(}\widehat{\mathbb{A}}_{j}\mathbbm{1}_{\hat{\eta}_{j-1}\geq
0}\,\hat{\eta}_{j-1}(x)+\mathbb{b}_{j}+\theta_{j}\big{)}_{k}=0\\}$

and observe that

$\\#\left(\hat{\omega}_{j,k}\setminus\omega_{j,k}\right)\leq
Q_{\widehat{\mathbb{A}}_{j},\hat{\eta}_{j-1},\theta_{j},k}$

since the number of not locally affine points added after applying the
activation function is less or equal than the number of $x\in\kappa$
satisfying $(\widehat{\mathbb{A}}_{j}\mathbbm{1}_{\hat{\eta}_{j-1}\geq
0}\hat{\eta}_{j-1}(x)+\mathbb{b}_{j})_{k}+(\theta_{j})_{k}=0$. Per Definition
2.3, we have that $\widehat{\mathbb{A}}_{j}$, $\hat{\eta}_{j-1}$, and
$\theta_{j}$ are independent random variables which allows the following
computation for $t\in\mathbb{N}$

$\displaystyle\mathbb{P}(\\#\left(\hat{\omega}_{j,k}\setminus\omega_{j,k}\right)\geq
t)$
$\displaystyle\leq\mathbb{P}_{\widehat{\mathbb{A}}_{j},\theta_{j},\hat{\eta}_{j-1}}(Q_{\widehat{\mathbb{A}}_{j},\theta_{j},k}\geq
t)$
$\displaystyle=\mathbb{E}_{\widehat{\mathbb{A}}_{j},\hat{\eta}_{j-1}}\mathbb{E}_{\theta_{j}}\mathbbm{1}_{Q_{\widehat{\mathbb{A}}_{j},\hat{\eta}_{j-1},\theta_{j},k}\geq
t}$
$\displaystyle=\mathbb{E}_{\widehat{\mathbb{A}}_{j},\hat{\eta}_{j-1}}\mathbb{P}\left(Q_{\widehat{\mathbb{A}}_{j},\hat{\eta}_{j-1},\theta_{j},k}\geq
t\right).$ (A.7)

By Assumption 2.7 we have that for all $x\in\Omega$

$\displaystyle 1$ $\displaystyle\geq|\nabla_{x}(\eta_{j})_{k}|$
$\displaystyle=\left|\nabla_{x}\left(\widehat{\mathbb{A}}_{j}\mathbbm{1}_{\hat{\eta}_{j-1}\geq
0}\hat{\eta}_{j-1}(x)+\mathbb{b}_{j}\right)_{k}\right|$
$\displaystyle=\left|\nabla_{x}\left(\widehat{\mathbb{A}}_{j}\mathbbm{1}_{\hat{\eta}_{j-1}\geq
0}\hat{\eta}_{j-1}(x)\right)_{k}\right|.$ (A.8)

We pick $v,w\in\mathbb{R}^{d}$ with $\|v\|=1$ such that $\kappa=\\{tv+w\colon
t\in[0,\mathcal{L}(\kappa)]\\}$ where $\mathcal{L}(\kappa)$ is the length of
$\kappa$. Next, we define

$\displaystyle h\colon[0,\mathcal{L}(\kappa)]$ $\displaystyle\to[0,1]^{d}$
$\displaystyle t$
$\displaystyle\mapsto\left(\widehat{\mathbb{A}}_{j}\mathbbm{1}_{\hat{\eta}_{j-1}\geq
0}\,\hat{\eta}_{j-1}(tv+w)\right)_{k}\,.$

Then, by the chain rule and (A.8), we have that $\|h^{\prime}\|_{\infty}\leq
1$. Moreover,

$\displaystyle
Q_{\widehat{\mathbb{A}}_{j},\hat{\eta}_{j-1},\theta_{j},k}=\\#\\{x\in[0,\mathcal{L}(\kappa)]:h(x)+(\theta_{j})_{k}=0\\}.$

Since by construction $(\theta_{j})_{k}$ is uniformly distributed on the
interval
$[-\lambda\varepsilon_{j}|(u_{j}^{b})_{k}|/2,\lambda\varepsilon_{j}|(u_{j}^{b})_{k}|/2]$,
we conclude from Remark 3.2 that

$\displaystyle\mathbb{P}\left(Q_{\widehat{\mathbb{A}}_{j},\hat{\eta}_{j-1},\theta_{j},k}\geq
t\right)\leq\frac{2\mathcal{L}(\kappa)}{\lambda}\cdot\left(\varepsilon_{j}\,q\,|(u_{j}^{b})_{k}|\right)^{-1}$

applied to (A.7) yields (3.1).

Part 2: (Proof of (3.2)) We recall the following upper bound on the total
number of pieces of the realisation of an arbitrary NN:

######  Theorem A.4 ([45]).

Let $L\in\mathbb{N}$ and $\sigma$ be piecewise affine with $p$ pieces. Then,
for every NN $\Phi$ with $d=1,N_{L}=1$ and $N_{1},\ldots,N_{L-1}\leq N$, we
have that $\mathrm{R}(\Phi)$ has at most $(pN)^{(L-1)}$ affine pieces.

Invoking Theorem A.4 results in
$\max_{k}\\{\\#\left(\hat{\omega}_{j,k}\setminus\omega_{j,k}\right)\\}\leq
N^{j-1}$ for $j=1,\dots,L$. We observe that

$\displaystyle\mathbb{E}(\\#(\hat{\omega}_{j,k}\setminus\omega_{j,k}))$
$\displaystyle=\sum_{q=1}^{N^{j}}\bigg{[}\mathbb{P}(\\#\left(\hat{\omega}_{j,k}\setminus\omega_{j,k}\right)\geq
q)-\mathbb{P}(\\#\left(\hat{\omega}_{j,k}\setminus\omega_{j,k}\right)\geq
q+1)\bigg{]}q$
$\displaystyle=\sum_{q=1}^{N^{j}}\mathbb{P}(\\#\left(\hat{\omega}_{j,k}\setminus\omega_{j,k}\right)\geq
q)q-\sum_{q=2}^{N^{j}+1}\mathbb{P}(\\#\left(\hat{\omega}_{j,k}\setminus\omega_{j,k}\right)\geq
q)(q-1)$
$\displaystyle=\sum_{q=1}^{N^{j}}\mathbb{P}(\\#\left(\hat{\omega}_{j,k}\setminus\omega_{j,k}\right)\geq
q)(q-(q-1))\,.$

Next, we invoke (3.1) and obtain that

$\displaystyle\mathbb{E}(\\#(\hat{\omega}_{j,k}\setminus\omega_{j,k}))$
$\displaystyle\leq\sum_{q=1}^{N^{j}}\frac{2\mathcal{L}(\kappa)}{\lambda\varepsilon_{j}|(u_{j}^{b})_{k}|q}(q-(q-1))$
$\displaystyle=\frac{2\mathcal{L}(\kappa)}{\lambda\varepsilon_{j}|(u_{j}^{b})_{k}|}\sum_{q=1}^{N^{j}}\frac{1}{q}$
$\displaystyle\leq\frac{2\mathcal{L}(\kappa)}{\lambda\varepsilon_{j}|(u_{j}^{b})_{k}|}\left(1+\int_{1}^{N^{j}}\frac{1}{x}dx\right)$
$\displaystyle\leq\frac{2\mathcal{L}(\kappa)}{\lambda\varepsilon_{j}|(u_{j}^{b})_{k}|}(1+\ln(N^{j}))\leq\frac{\mathcal{L}(\kappa)}{2\lambda\varepsilon_{j}|(u_{j}^{b})_{k}|}(\ln(eN^{j}))\leq\frac{2\mathcal{L}(\kappa)}{\lambda\varepsilon_{j}|(u_{j}^{b})_{k}|}(\ln(N^{j+1}))\,,$

where we used that $e\leq 3\leq N$. This completes the proof of (3.2). ∎

###  A.3 Proof of Theorem 3.4

###### Proof.

Part 1: We start by proving the result for $j^{\prime}=L$. Using the notation
of Proposition 3.3, we start by assuming that

$\displaystyle\mathrm{pieces}(\mathrm{R}(\widehat{\Phi}^{\varepsilon}),\kappa))\leq
1+\sum_{j\in[L-1]}\sum_{k\in
N_{j}}\\#\left(\hat{\omega}_{j,k}\setminus\omega_{j,k}\right),$ (A.9)

if $\omega_{j,k}$ is chosen as $\omega_{j,k}=\bigcup_{k\in
N_{j-1}}\hat{\omega}_{j-1,k}$ for $j>1$ and $\omega_{1,k}=\varnothing$.

We will prove (A.9) at the end this step of this proof to not distract from
the main argument. By (A.9) and the linearity of the expected value, we
obtain,

$\mathbb{E}(\mathrm{pieces}(\mathrm{R}(\widehat{\Phi}^{\varepsilon}),\kappa)))\leq
1+\sum_{j\in[L-1]}\sum_{k\in
N_{j}}\mathbb{E}(\\#\left(\hat{\omega}_{j,k}\setminus\omega_{j,k}\right))\,.$

By Assumption 2.6, it holds that if $(u_{j}^{b})_{k}=0$, then
$\mathrm{R}((\widehat{\Phi}^{\varepsilon})_{j})_{k}=0$ and therefore
$\\#\left(\hat{\omega}_{j,k}\setminus\omega_{j,k}\right)=0$. Hence, invoking
Proposition 3.3, we conclude, using Assumption 2.6, that

$\displaystyle\mathbb{E}(\mathrm{pieces}(\mathrm{R}(\widehat{\Phi}^{\varepsilon}),\kappa)))$
$\displaystyle\leq
1+2\mathcal{L}(\kappa)\lambda^{-1}\sum_{j=1}^{L-1}\sum_{k\in
N_{j}}\frac{1}{\varepsilon_{j}}|(u_{j}^{b})_{k}|^{\dagger}\log(N^{j+1})$
$\displaystyle\leq
1+2\mathcal{L}(\kappa)\lambda^{-1}\sum_{j=1}^{L-1}\sum_{k\in
N_{j}}\frac{j+1}{\varepsilon_{j}}|(u_{j}^{b})_{k}|^{\dagger}\log(N)$
$\displaystyle\leq
1+2\mathcal{L}(\kappa)\lambda^{-1}\frac{L}{\hat{\epsilon}_{L}}\sum_{j=1}^{L-1}\sum_{k\in
N_{j}}|(u_{j}^{b})_{k}|^{\dagger}\log(N)$ $\displaystyle\leq
1+2\mathcal{L}(\kappa)\lambda^{-1}c_{0}\frac{L}{\hat{\epsilon}_{L}}N^{\nu}\log(N)\,.$

We complete the proof of the first part by showing (A.9). Note that since

$\displaystyle\mathrm{R}(\widehat{\Phi}^{\varepsilon})_{k^{\prime}}=\sum_{k\in
N_{L-1}}(\widehat{\mathbb{A}}_{L})_{k^{\prime}k}\,\varrho(\mathrm{R}((\widehat{\Phi}^{\varepsilon})_{L-1})_{k})+(\widehat{\mathbb{b}}_{L})_{k^{\prime}}\,,$
(A.10)

we have that

$\displaystyle\mathrm{pieces}(\mathrm{R}(\widehat{\Phi}^{\varepsilon}),\kappa))\leq
1+\\#\bigcup_{k\in N_{L-1}}\hat{\omega}_{L-1,k}\,,$

where we use that a piecewise affine function on a line has exactly one more
affine piece than the smallest number of distinct points where the function is
not affine.

We will show by induction over $L$ that

$\displaystyle\\#\bigcup_{k\in
N_{L-1}}\hat{\omega}_{L-1,k}\leq\sum_{j\in[L-1]}\sum_{k\in
N_{j}}\\#\left(\hat{\omega}_{j,k}\setminus\omega_{j,k}\right)$ (A.11)

which then yields the claim. For $L=2$, the claim follows directly by the
choice $\omega_{1,k}=\varnothing$.

Now let $L>2$, then

$\displaystyle\bigcup_{k\in N_{L-1}}\hat{\omega}_{L-1,k}$
$\displaystyle=\bigcup_{k\in
N_{L-1}}\left(\hat{\omega}_{L-1,k}\setminus\omega_{L-1,k}\right)\cup\omega_{L-1,k}$
$\displaystyle=\omega_{L-1,1}\cup\bigcup_{k\in
N_{L-1}}\left(\hat{\omega}_{L-1,k}\setminus\omega_{L-1,k}\right),$ (A.12)

since $\omega_{L-1,1}=\omega_{L-1,k}$ for all $k\in N_{L-1}$ by definition. We
have that

$\displaystyle\mathrm{R}((\widehat{\Phi}^{\varepsilon})_{L-1})_{k^{\prime}}=\sum_{k\in
N_{L-2}}(\widehat{\mathbb{A}}_{L-1})_{k^{\prime}k}\,\varrho(\mathrm{R}(\Phi_{L-2})_{k})+(\widehat{\mathbb{b}}_{L-1})_{k^{\prime}}\,,$

and hence

$\displaystyle\omega_{L-1,1}\subset\bigcup_{k\in
N_{L-2}}\hat{\omega}_{L-2,k}\,.$ (A.13)

Using (A.13) and applying the induction hypothesis (A.11) to (A.12), we obtain
that

$\displaystyle\\#\bigcup_{k\in N_{L-1}}\hat{\omega}_{L-1,k}$
$\displaystyle\leq\\#\omega_{L-1,1}+\\#\bigcup_{k\in
N_{L-1}}\left(\hat{\omega}_{L-1,k}\setminus\omega_{L-1,k}\right)$
$\displaystyle\leq\\#\bigcup_{k\in
N_{L-2}}\hat{\omega}_{L-2,k}+\\#\bigcup_{k\in
N_{L-1}}\left(\hat{\omega}_{L-1,k}\setminus\omega_{L-1,k}\right)$
$\displaystyle\leq\sum_{j\in[L-2]}\sum_{k\in
N_{j}}\\#\left(\hat{\omega}_{j,k}\setminus\omega_{j,k}\right)+\sum_{k\in
N_{L-1}}\\#\left(\hat{\omega}_{L-1,k}\setminus\omega_{L-1,k}\right)\,,$

which yields (A.11).

Part 2: Let $j^{\prime}>0$. We have that

$\displaystyle\mathrm{R}(\widehat{\Phi}^{\varepsilon})=\mathrm{R}(\widehat{\Phi}^{\varepsilon,c}_{j^{\prime}})\circ\varrho\circ\mathrm{R}(\widehat{\Phi}^{\varepsilon}_{j^{\prime}}),$

where

$\displaystyle\widehat{\Phi}^{\varepsilon,c}_{j^{\prime}}=\big{(}(\widehat{\mathbb{A}}_{j^{\prime}+1},\widehat{\mathbb{b}}_{j^{\prime}+1}),\dots,(\widehat{\mathbb{A}}_{L},\widehat{\mathbb{b}}_{L})\big{)}.$
(A.14)

For $q,p\in\mathbb{N}$ and two functions $f:\kappa\to\mathbb{R}^{q}$ and
$g:\mathbb{R}^{q}\to\mathbb{R}^{p}$ such that $f$ has $r\in\mathbb{N}$ affine
pieces on $\kappa$, and $g$ has at most $s\in\mathbb{N}$ many pieces along all
possible lines it holds that $g\circ f\colon\kappa\to\mathbb{R}^{q}$ has at
most $rs$ affine pieces, since on the image of each affine piece of $f$ no
more than $s$ affine pieces can be generated by $g$. Note that
$\varrho\colon\mathbb{R}^{N_{j^{\prime}}}\to\mathbb{R}^{N_{j^{\prime}}}$ has
at most $2N$, pieces along each line, since every coordinate of a line can
change sign at most once. Invoking again Proposition A.4 as well as Part 1 of
this proof and the monotonicity of the expected value therefore yields that

$\displaystyle\mathbb{E}(\mathrm{pieces}(\mathrm{R}(\widehat{\Phi}^{\varepsilon}),\kappa))$
$\displaystyle\leq(2N)^{L-j^{\prime}-1}\cdot
2N_{j^{\prime}}\cdot\left(1+2\mathcal{L}(\kappa)\lambda^{-1}c_{0}\frac{j^{\prime}}{\hat{\epsilon}_{j^{\prime}}}N^{\nu}\log(N)\right)$
$\displaystyle\leq(2N)^{L-j^{\prime}}\cdot\left(1+2\mathcal{L}(\kappa)\lambda^{-1}c_{0}\frac{j^{\prime}}{\hat{\epsilon}_{j^{\prime}}}N^{\nu}\log(N)\right),$

which completes the proof. ∎