# Activation function design for deep networks: linearity and effective
initialisation

$\dagger{\ddagger}$ M. Murray, $\dagger$ V. Abrol, $\dagger{\ddagger}$ J.
Tanner
$\dagger$ Mathematical Institute, University of Oxford, UK
${\ddagger}$ The Alan Turing Institute, London, UK
<EMAIL_ADDRESS>

###### Abstract

The activation function deployed in a deep neural network has great influence
on the performance of the network at initialisation, which in turn has
implications for training. In this paper we study how to avoid two problems at
initialisation identified in prior works: rapid convergence of pairwise input
correlations, and vanishing and exploding gradients. We prove that both these
problems can be avoided by choosing an activation function possessing a
sufficiently large linear region around the origin, relative to the bias
variance $\sigma_{b}^{2}$ of the network’s random initialisation. We
demonstrate empirically that using such activation functions leads to tangible
benefits in practice, both in terms test and training accuracy as well as
training time. Furthermore, we observe that the shape of the nonlinear
activation outside the linear region appears to have a relatively limited
impact on training. Finally, our results also allow us to train networks in a
new hyperparameter regime, with a much larger bias variance than has
previously been possible.

_Keywords_ Activation Function Design, Deep Learning, Initialisation, Random
Networks.

##  1 Introduction

The preactivations $\mathbf{h}^{(l)}\in\mathbb{R}^{N_{l}}$ and activations
$\mathbf{z}^{(l)}(\mathbf{x})\in\mathbb{R}^{N_{l}}$ of an input vector
$\mathbf{x}\in\mathbb{R}^{N_{0}}$ at each layer $l\in[L]$ of a fully-
connected, deep, feedforward neural network, with weight matrices
$\mathbf{W}^{(l)}\in\mathbb{R}^{N_{l}\times N_{l-1}}$ and bias vectors
$\mathbf{b}^{(l)}\in\mathbb{R}^{N_{l}}$, are computed by the following pair of
recurrence relations,

$\mathbf{z}^{(l)}(\mathbf{x})=\mathbf{W}^{(l)}\mathbf{h}^{(l-1)}+\mathbf{b}^{(l-1)},\;\mathbf{h}^{(l)}(\mathbf{x})=\phi(\mathbf{z}^{(l)}(\mathbf{x})),$
(1)

where $\phi:\mathbb{R\rightarrow R}$ is the pointwise nonlinearity, referred
to as the activation function, deployed at each layer. The impact of the
choice of activation function $\phi$ on a wide range of a properties,
including for instance the trainability and generality of the network, is an
active area of research [1, 2, 3, 4, 5, 6]. In this paper we study
specifically how the choice of activation function impacts the performance of
the network at initialisation, and the subsequent implications of this on
training. We build on prior work [7, 8, 9, 10] concerning randomly initialised
deep neural networks (DNNs), focusing on two objectives: first, ensuring
deeper information propagation by slowing the rate of convergence of pairwise
input correlations, and second, avoiding vanishing and exploding gradients
[11] by at least approximately achieving dynamical isometry. In [7] estimates
for how a) the length of the preactivations of an input and b) the correlation
between the preactivations of two distinct inputs evolve with depth were
derived. In both [7, 8] the authors investigated how the choice of the
activation function $\phi$ and the variance of the random weights
$\sigma_{w}^{2}$ and biases $\sigma_{b}^{2}$ effect the lengths of and
correlations between preactivations. A key finding of [8] was the
identification of a condition on $(\sigma_{w}^{2},\sigma_{b}^{2})$, termed
initialisation on the edge of chaos (EOC), which is necessary in order to
avoid asymptotic exponentially fast convergence of both the lengths and
correlations. Furthermore, [10] identified for a broad range of activation
function that by letting $\sigma_{b}^{2}\rightarrow 0$ then the rate of
convergence in correlation can be reduced away from the fixed point, i.e.,
non-asymptotically. We review this line of research in detail in Section 2.1.
A related line of work [9, 12, 13] proposed that the problem of vanishing and
exploding gradients can be mitigated by ensuring that the spectrum of the
input-output Jacobian of the network is concentrated around one, a property
termed dynamical isometry. In addition, it was shown that initialisation on
the EOC implies that the mean of the spectrum of the input-output Jacobian is
one. Therefore, dynamical isometry can at least approximately be achieved by
making, through careful selection of $(\phi,\sigma_{w}^{2},\sigma_{b}^{2})$,
the variance $\sigma_{JJ^{T}}^{2}$ of the spectrum close to zero. More
specifically, the authors highlighted how this can be achieved for particular
activation functions by using an orthogonal initialisation scheme for the
weights and again by letting $\sigma_{b}^{2}\rightarrow 0$. A number of
algorithmic and architectural methods have been proposed and deployed in
practice to overcome the problem of vanishing and exploding gradients across a
variety of contexts, see for instance [1, 5, 14, 15, 16, 17, 18, 19]. For
brevity we do not review these here, but remark that the line of work
described above can be viewed as complimentary to these methods, focusing
specifically on the network at initialisation.

Letting $\sigma_{b}^{2}\rightarrow 0$ is unsatisfactory for a number of
reasons, first, $\sigma_{b}^{2}\rightarrow 0$ typically results in the length
of the preactivations converging towards zero with depth. Second, and as
demonstrated by our experiments in Section 4, more affine initialisations with
a small but not overly small $\sigma_{b}^{2}$ appear to provide the best
training outcomes in practice. Third and finally, this strategy does not
leverage our freedom in designing the activation function deployed in order to
achieve these two goals. In this paper we introduce and study a novel set of
scaled, bounded activation functions which are linear in some neighbourhood of
the origin. For functions in this set we derive upper bounds controlling both
the rate of convergence of correlations between preactivations, as well as the
variance of the spectrum of the input-output Jacobian. Notably these bounds
depend on the ratio $\sigma_{b}^{2}/a^{2}$, where $a$ denotes the size of the
linear region of the activation function, and approach zero as
$\sigma_{b}^{2}/a^{2}\rightarrow 0$. The theory presented in this paper
therefore provides a rigorous explanation for the observation that networks
equipped with an activation function which approximates the identity at the
origin often train well (see for example [20]). Furthermore, our experiments
demonstrate the benefits of deploying these scaled, bounded activations in
practice, in terms of train and test accuracy, training time and allowing us
to train networks with a much larger bias variance than was previously
possible.

The structure of this paper is as follows, in Section 2 we review the
foundational results of prior works and summarise the contributions of this
paper, in particular Theorem 2.1. In Section 3 we derive Theorem 2.1, then in
Section 4 we investigate the implications of our theory for training deep
networks in practice.

##  2 Principles for initialising deep networks

As referenced to in Section 1, the prior works on which our results are based
can be divided into two distinct themes, one concerned with the dynamics of
the preactivation correlations in the forward pass [7, 8, 10], and the other
dynamical isometry [9, 12, 13]. We review these two themes in detail in
Sections 2.1 and 2.2 respectively, then in Section 2.3 we introduce and
present the key contribution of this paper, Theorem 2.1. Note that in what
follows $Z,Z_{1},Z_{2}\sim\mathcal{N}(0,1)$ are considered to be independent
and identically distributed standard Gaussian random variables. Furthermore,
the standard Gaussian measures in one and two dimensions are be denoted by
$\gamma$ and $\gamma^{(2)}$ respectively.

###  2.1 Dynamics of the preactivation correlations in the forward pass

In [21] it was established that the outputs of certain random, feedforward,
single hidden layer neural networks converge in distribution to centred
Gaussian processes as the width of the hidden layer goes to infinity. More
recently in [22, 23], this result was extended to random multilayer neural
networks, with the preactivations at each layer being found to converge in
distribution to centred Gaussian processes in the large width limit. The
behaviour of a centred Gaussian process is fully described by its covariance
matrix. In particular, for networks whose forward pass is described by (1),
with the weights and biases at each layer $l\in[L]$ being mutually
independent, centred Gaussian random variables with variances
$\sigma_{w}^{2}/N_{l-1}$ and $\sigma_{b}^{2}$ respectively, and assuming that
the same activation function $\phi$ is deployed at each neuron, then the
kernel used to compute the entries of the covariance matrix of the Gaussian
process at the $l$th layer is defined by the following recurrence relation,

$\kappa^{(l)}(\mathbf{x}^{(\alpha)},\mathbf{x}^{(\beta)})=\sigma_{w}^{2}\mathbb{E}[\phi(z_{i}^{(l)}(\mathbf{x}^{(\alpha)}))\phi(z_{i}^{(l)}(\mathbf{x}^{(\beta)}))]+\sigma_{b}^{2}.$
(2)

Note here that $i\in[N_{l}]$ can be any entry of the preactivation. This
recurrence relation has been highlighted a number of times in a variety of
contexts, the one most relevant for this paper being a mean field
approximation used to study signal propagation in the forward pass [7, 8, 10].
In this work we do not focus on the correspondence between wide, random neural
networks and Gaussian processes as in [22, 23]. Rather, as in [7, 8, 10], we
adopt the Gaussian process or mean field approximation in order to better
understand and analyse key statistics of the network at initialisation. In
particular, and as highlighted in [7], in the infinite width limit one can
interpret the variance of the preactivations of an input at a given layer,
corresponding to (2) with $\alpha\neq\beta$, as the expected euclidean length
of said input at said layer. Adopting the notation of [7], then the sequence
of variances, or expected lengths, $(q_{\alpha}^{(l)})_{l=1}^{L}$, associated
with an arbitrary input $\mathbf{x}^{(\alpha)}$, are generated by the
following recurrence relation,

$\displaystyle
q_{\alpha}^{(l)}=V_{\phi}(q_{\alpha}^{(l-1)}):=\sigma_{w}^{2}\int_{\mathbb{R}}\phi\left(\sqrt{q_{\alpha}^{(l-1)}}z\right)^{2}d\gamma(z)+\sigma_{b}^{2}=\sigma_{w}^{2}\mathbb{E}\left[\phi\left(\sqrt{q_{\alpha}^{(l-1)}}Z\right)^{2}\right]+\sigma_{b}^{2}$
(3)

for $l\in[L]L$ with
$q_{\alpha}^{(0)}\vcentcolon=||\mathbf{x}^{(\alpha)}||_{2}^{2}$. As
demonstrated in Figure 1(a), analysis of the variance function
$V_{\phi}:\mathbb{R}_{\geq 0}\rightarrow\mathbb{R}_{\geq 0}$ provides insight
into the impact of $(\phi,\sigma_{w}^{2},\sigma_{b}^{2})$ on the dynamics of
the expected length of the preactivations with depth. Indeed, unless
$(\phi,\sigma_{w}^{2},\sigma_{b}^{2})$ is chosen appropriately then
$q_{\alpha}^{(l)}$ can either rapidly converge towards zero or diverge. This
problem of vanishing or exploding activation lengths can be mitigated by
choosing $(\phi,\sigma_{w}^{2},\sigma_{b}^{2})$ in order that there exists a
stable fixed point $q^{*}=V_{\phi}(q^{*})$ for which
$q_{\alpha}^{(l)}\rightarrow q^{*}$ for all inputs $\mathbf{x}^{(\alpha)}$,
see for example Tanh with $\sigma_{b}^{2}=0.1$ in Figure 1(a). We remark that
this problem was also studied in [24] in the context of finite width deep ReLU
networks.

(a)

(b)

Figure 1: (a) Variance functions for ReLU, Tanh and ELU with
$\sigma_{b}^{2}\in\\{0,0.1\\}$. With $\sigma_{b}^{2}=0.1$ the variance
function associated with ReLU has no fixed point and
$q^{(l)}\rightarrow\infty$ from any initial variance
$q^{(0)}\in\mathbb{R}_{\geq 0}$. With $\sigma_{b}^{2}=0$ the variance function
associated with ELU has a unique and marginally stable point $q^{*}=0$,
therefore $q^{(l)}\rightarrow 0$ from any initial variance
$q^{(0)}\in\mathbb{R}_{\geq 0}$. Finally, with $\sigma_{b}^{2}=0.1$, the
variance function associated with Tanh has a unique and stable fixed point
$q^{*}>0$, therefore $q^{(l)}\rightarrow q^{*}$ from any initial variance
$q^{(0)}\in\mathbb{R}_{\geq 0}$. (b) Correlation functions associated with
Relu, Tanh and ELU with $\sigma_{b}^{2}\in\\{0,0.1\\}$. For ReLU, and as
illustrated by (a), for there to exist a fixed point of the variance function
it is necessary that $\sigma_{b}^{2}=0$. Note that without a fixed point
$q^{*}$ it is not possible to define a correlation function which is fixed
with depth. Relative to the other activation functions considered, the
correlation function of ReLU is further from the identity function, and as a
result correlations converge faster with the depth. This can also be observed
in Figure 10 in Appendix B.6. Note that in both (a) and (b) $\sigma_{w}^{2}$
is chosen so that $\chi_{1}=1$, see (7).

In [7] the covariance and correlation between preactivations of different
inputs was also analysed and a recurrence relation for the covariance derived,
which is equivalent to (2) with $\alpha\neq\beta$. In the infinite width limit
the covariance can be interpreted as the inner product between the
preactivations of two inputs, and the correlation as the angle. The evolution
of these quantities with depth has important implications for the performance
of the network at initialisation as well as subsequent training. Denoting the
covariance between the preactivations of $\mathbf{x}^{(\alpha)}$ and
$\mathbf{x}^{(\beta)}$ at the $l$th layer as $q_{\alpha\beta}^{(l)}$, with
respective variances $q_{\alpha}^{(l)}$ and $q_{\beta}^{(l)}$, then it was
shown in [7] that

$\displaystyle
q_{\alpha\beta}^{(l)}=\sigma_{w}^{2}\int_{\mathbb{R}^{2}}\phi(u_{1})\phi(u_{2})d\gamma^{(2)}(z_{1},z_{2})+\sigma_{b}^{2},$
(4) $\displaystyle
u_{1}:=\sqrt{q_{\alpha}^{(l-1)}}z_{1},\;u_{2}:=\sqrt{q_{\beta}^{(l-1)}}(\rho_{\alpha\beta}^{(l-1)}z_{1}+\sqrt{1-(\rho_{\alpha\beta}^{(l-1)})^{2}}z_{2}),$

where
$\rho_{\alpha\beta}^{(l-1)}=q_{\alpha\beta}^{(l-1)}/(q_{\alpha}^{(l-1)}q_{\beta}^{(l-1)})$
is the correlation at the previous layer. If the variance function defined in
(3) has a fixed point $q^{*}>0$ then an additional key benefit beyond that
already discussed is that the analysis of the evolution of the covariance and
correlation with depth can be simplified. In particular, assuming that the
inputs are normalised so that $q_{\alpha}^{(0)}=q_{\beta}^{(0)}=q^{*}>0$, then
$q_{\alpha}^{(l)}=q_{\beta}^{(l)}=q^{*}$ for all $l\in[L]$. In this setting
the sequence of correlations $(\rho_{\alpha\beta}^{(l)})_{l=1}^{L}$ are
generated by the following recurrence relation,

$\displaystyle\rho_{\alpha\beta}^{(l)}=R_{\phi,q^{*}}(\rho_{\alpha\beta}^{(l-1)}):=\frac{q_{\alpha\beta}^{(l)}}{q^{*}}=\frac{\sigma_{w}^{2}}{q^{*}}\int_{\mathbb{R}}\phi(u_{1})\phi(u_{2})d\gamma^{(2)}(z_{1},z_{2})+\frac{\sigma_{b}^{2}}{q^{*}}=\frac{\sigma_{w}^{2}}{q^{*}}\mathbb{E}[\phi(U_{1})\phi(U_{2})]+\frac{\sigma_{b}^{2}}{q^{*}}.$
(5)

Here we refer to $R_{\phi,q^{*}}:[-1,1]\rightarrow[-1,1]$ as the correlation
function, and $U_{1}:=\sqrt{q^{*}}Z_{1}$, $U_{2}:=\sqrt{q^{*}}(\rho
Z_{1}+\sqrt{1-\rho^{2}}Z_{2})$, where $\rho$ is the input argument
$R_{\phi,q^{*}}(\rho)$, are dependent Gaussian random variables with
$U_{1},U_{2}\sim\mathcal{N}(0,q^{*})$.

Analysing the univariate recurrence relation
$\rho^{(l)}=R_{\phi,q^{*}}(\rho^{(l-1)})$ allows for the identification of
both depth limits, beyond which information cannot propagate, as well as
issues around stability, in a manner analogous to that of the dynamical
systems perspective given in [25]. In order to appreciate these points,
observe, as in [7], that by construction as long as the correlation function
is well defined on $[-1,1]$ then one is always a fixed point,

$\displaystyle
R_{\phi,q^{*}}(1)=\frac{\sigma_{w}^{2}}{q^{*}}\mathbb{E}[\phi(U_{1})^{2}]+\frac{\sigma_{b}^{2}}{q^{*}}=\frac{V(q^{*})}{q^{*}}=1.$

The general shape of the correlation function, as well as the stability of the
fixed point at one and the existence of other fixed points, is determined by
the choice of $(\phi,\sigma_{w}^{2},\sigma_{b}^{2})$. This can be observed in
Figure 1(b), in which the correlation functions of three commonly used
activation functions are plotted side by side for comparison. Figure 2(a)
shows the correlation function for three different values of $\sigma_{w}^{2}$
with $\sigma_{b}^{2}=0.1$ and $\phi(z)=\text{htanh}(z)$ fixed. For
$\sigma_{w}^{2}=0.63$ and $\sigma_{w}^{2}=1.26$ it is clear from Figure 2(a)
that $\rho^{*}=1$ is the only positive fixed point. Furthermore, for these two
choices of $\sigma_{w}^{2}$ it is also clear, by inspection, that $\rho^{*}=1$
is a stable and marginally stable fixed point respectively, and that from any
initial correlation the sequence of correlations $\rho^{(l)}\rightarrow 1$
with depth. As a result, in the infinite width limit any pair of inputs, no
matter how large the angle between them is, are mapped to the same point by
the network asymptotically with depth. However, for $\sigma_{w}^{2}=0.63$ a
new, stable fixed point $\rho^{*}<1$ is introduced. In this case, then for any
initial correlation the sequence $\rho^{(l)}$ converges to the limit
$\rho^{*}<1$. Therefore, in the infinite width limit any pair of inputs, no
matter how small the initial angle between them is, are mapped to different
points by the network asymptotically with depth. Figure 2(a) therefore
illustrates the existence of two very distinct regimes for a given choice of
activation function, one ordered, in which all inputs are mapped
asymptotically to the same output, and the other chaotic, in which even
arbitrarily small perturbations of an input lead to a different output. Figure
2(b) illustrates the limiting $\rho^{*}$ throughout the
$(\sigma_{w},\sigma_{b})$ plane for $\phi=\text{htanh}(\cdot)$.

(a)

(b)

(c)

(d)

Figure 2: (a) Correlation map for htanh activation function, plotted on the
interval $[0,1]$ to better display fixed point behaviour (there are no fixed
points in the interval $[-1,0]$). When $\chi\leq 1$ then one is the unique
fixed point of the correlation function. Furthermore, for $\chi<1$ then one is
a stable fixed point, and for $\chi=1$ then one is a marginally stable fixed
point. When $\chi>1$ then although one is still a fixed point it is unstable,
and a new stable fixed point less than one is introduced. (b) Stable fixed
point $\rho^{*}$ as a function of $(\sigma_{w}^{2},\sigma_{b}^{2})$ for the
tanh activation function. (c) The EOC critical line $\chi=1$ for different
activation functions. (d) $\sigma_{w}^{2}$ as a function of $q^{*}$, as
$\sigma_{b}^{2}\rightarrow 0$ then $q^{*}\rightarrow 0$.

In order to understand and separate these two regions [7, 8] studied the slope
of $R_{\phi,q^{*}}(\rho)$ at $\rho=1$: for a sufficiently smooth activation
function $\phi$ the derivative is given by

$R_{\phi,q^{*}}^{\prime}(\rho)=\sigma_{w}^{2}\mathbb{E}[\phi^{\prime}(U_{1})\phi^{\prime}(U_{2})],$
(6)

the slope at $\rho=1$ is therefore

$\chi_{1}:=R_{\phi,q^{*}}^{\prime}(1)=\sigma_{w}^{2}\mathbb{E}[\phi^{\prime}(U_{1})^{2}].$
(7)

The set of pairs $(\sigma_{w}^{2},\sigma_{b}^{2})\in\mathbb{R}_{\geq
0}\times\mathbb{R}_{\geq 0}$ such that $\chi_{1}=1$ is colloquially referred
to in the literature as the edge of chaos (EOC) and in Figure 2(c) this set,
or curve, is plotted for three different activation functions. In [8] it was
proved that if $(\sigma_{w}^{2},\sigma_{b}^{2})$ does not lie on the EOC then
the sequences generated by the variance and correlation functions, (3) and (5)
respectively, approach their respective fixed points asymptotically
exponentially fast. As a result, the EOC curve can be interpreted as a
transition boundary between order, corresponding to $\chi_{1}<1$, in which
pairwise correlations converge asymptotically exponentially fast to $1$ with
depth, and chaos, corresponding to $\chi_{1}>1$, in which even pairwise input
correlations that are arbitrarily close to $1$ diverge asymptotically
exponentially fast with depth. The dynamics of the preactivation correlations
can also be understood and interpreted in terms of a network’s sensitivity to
perturbations of the input, with networks initialised in the chaotic regime
carrying the risk of being overly sensitive, and those initialised in the
ordered regime being insensitive. These observations are also of practical
relevance for training, indeed experimental evidence [8] illustrates that
initialising closer to the EOC consistently results in reduced training times.

To recap the discussion presented so far, rapid convergence in correlations to
a fixed point appears to result in poor training outcomes, with the network
being either highly sensitive or insensitive to perturbations of the input. To
avoid asymptotically exponentially fast convergence it is necessary to choose
$(\sigma_{w}^{2},\sigma_{b}^{2})$ so that $\chi_{1}=1$. The specific
asymptotic rate however depends on the activation function deployed. In [10],
and under the assumption that initialisation is on the EOC, the asymptotic
convergence rate of the sequence of correlations for different families of
activation functions was analysed. A key finding of this work is that while
ReLU like activation functions have asymptotic convergence
$|1-\rho^{(l)}|=\mathcal{O}(l^{-2})$, a broad class of smooth activation
functions, including, tanh, elu and swish, have convergence
$|1-\rho^{(l)}|=\mathcal{O}(l^{-1})$. However, it seems highly reasonable that
avoiding fast non-asymptotic, i.e., away for $\rho^{*}=1$, convergence of
correlations is also potentially of value. Figure 1(b) illustrates, assuming
initialisation on the EOC, that certain combinations of
$(\phi,\sigma_{b}^{2})$ can result in a correlation function which is closer
to the identity function. This in turn implies a slower convergence of the
correlation sequence throughout the layers, rather than just at the depth
limit. This idea was promoted in [10, Appendix B.2], where it was suggested
that one should choose a small $\sigma_{b}^{2}$ and an activation function
which is approximately linear, for example, one which is the weighted sum of a
linear and nonlinear activation function.

###  2.2 Avoiding vanishing and exploding gradients via dynamical isometry

A key ingredient behind the success of neural networks is that their modular,
compositional structure enables easy computation of the gradients of the
network output with respect to the model parameters. The computation of these
gradients can efficiently be performed using the backpropagation algorithm
[26]. Backpropagation, also referred to as the backward pass of the network,
generates a sequence of vectors $(\delta^{(l)}(\mathbf{x}))_{l=1}^{L}$ at each
layer with respect to an input $\mathbf{x}\in\mathbb{R}^{N}$, which measure
the error associated with minimising the loss $\mathcal{L}(\theta,x)$. Here
$\theta$ denotes the trainable network parameters,
$\theta\vcentcolon=\bigcup_{l=1}^{L}\\{\mathbf{W}^{(l)}\\}\cup\\{\mathbf{b}^{(l)}\\}$.
For a suitably smooth loss function $\mathcal{L}(\theta,x)$, the
backpropagation error vector $\delta^{(l)}(\mathbf{x})$, and the gradient of
$\mathcal{L}(\theta,x)$ with respect to biases $b_{j}^{(l)}$ and weights
$w^{(l)}_{j,i}$, at each layer $l\in[L]$ can be computed using the following
recurrence equations,

$\displaystyle\delta^{(L)}(\mathbf{x})$
$\displaystyle=\mathbf{D}^{(L)}(\mathbf{x})\nabla_{\mathbf{h}^{(L)}}\mathcal{L}(\theta,\mathbf{x}),$
(8) $\displaystyle\delta^{(l)}(\mathbf{x})$
$\displaystyle=\left(\mathbf{D}^{(l)}(\mathbf{x})\mathbf{W}^{(l+1)}\right)^{T}\delta^{(l+1)}(\mathbf{x}),$
$\displaystyle\frac{\partial\mathcal{L}(\theta,\mathbf{x})}{\partial
b_{j}^{(l)}}$ $\displaystyle=\delta_{j}^{(l)}(\mathbf{x}),$
$\displaystyle\frac{\partial\mathcal{L}(\theta,\mathbf{x})}{\partial
w_{j,i}^{(l)}}$
$\displaystyle=h_{i}^{(l-1)}(\mathbf{x})\delta_{j}^{(l)}(\mathbf{x}).$

Here, for each $l\in[L]$, $\mathbf{D}^{(l)}(\mathbf{x})$ is a diagonal matrix
with
$D_{ii}^{(l)}(\mathbf{x})\vcentcolon=\phi^{{}^{\prime}}(z_{i}^{(l)}(\mathbf{x}))$
for all $i\in[N_{l}]$. Consequently, the weights and biases of the network can
be iteratively updated via gradient descent type algorithms, or some variant
thereof, in order that the forward pass improves its performance, as measured
typically by averaging $\mathcal{L}(\theta,\mathbf{x})$ over batches of the
training data, at executing the task in hand. We emphasise that in practice,
instead of updating the parameters with respect to a single data point,
gradients are calculated and averaged over batches of the training data. It is
important therefore not to confuse the input-output Jacobian of the network,
which we will discuss presently, with the Jacobian of the loss surface over
the training data, or batches of the training data.

We now turn our attention to reviewing the findings of the literature
concerning the propagation of $\delta^{(l)}$ in the backwards pass. Analogous
to the recurrence relation for the variance function (3) in the forward pass,
in [8] a mean field approximation for $\delta^{(l)}$ was developed.
Specifically, the recurrence relation for the variance,
$\tilde{q}^{l}\vcentcolon=\mathbb{E}[(\delta_{i}^{l})^{2}]$, of entries of the
error vectors $\delta_{i}^{l}$ was shown to be

$\tilde{q}^{l}=\frac{N^{l+1}}{N^{l}}\tilde{q}^{l+1}\chi_{1}$ (9)

with $\chi_{1}$ defined as in (7). Although certain assumptions needed for
this approximation are undesirable, notably the weights used during forward
propagation are drawn independently from the weights used in backpropagation,
empirical evidence indicates that it is nonetheless a useful model [8]. In
particular, unlike the analysis of the gradients used to update the network
parameters. The reappearance of $\chi_{1}$ in (9) shows the benefit of
selecting $(\phi,\sigma_{w}^{2},\sigma_{b}^{2})$ such that $\chi=1$, as this
ensures, at least in expectation, that the magnitudes of the error vectors are
stable during the backward pass, i.e., they neither converge to zero or
diverge.

Unfortunately, initialising on the EOC suffices only to ensure that vanishing
and exploding gradients are avoided on average, and in practice does not
guarantee good training performance. In particular, if, as it propagates
backwards, an error vector is projected onto a smaller and smaller subspace
then the directions in which the parameters in the lower layers can be updated
becomes limited, thereby resulting in model degeneracy [12, 13]. In order to
improve training, in [12, 13] the authors proposed that the product of
Jacobians at each layer should act as an isometry on as large a subspace as
possible, implying that the singular values of the input-output Jacobian at
any layer should concentrate around one. As in [9, 13], we denote the input-
output Jacobian of the network as
$\mathbf{J}:\mathbb{R}^{N_{0}}\rightarrow\in\mathbb{R}^{N_{0}\times N_{L}}$,

$\mathbf{J}(\mathbf{x})=\prod_{l=1}^{L}\mathbf{D}^{(l)}(\mathbf{x})\mathbf{W}^{(l)}.$
(10)

Comparing (10) with (8), then the Jacobian of $\mathbf{x}$ can be viewed as
the transpose of the product of linear backpropagation operators used to
compute the error vectors associated with $\mathbf{x}$ at each layer. In [9,
13] the authors considered two initialisation schemes, in both the biases
$\mathbf{b}^{(l)}_{i}$ at each layer are mutually independent and identically
distributed Gaussians with mean $0$ and variance $\sigma_{b}^{2}$. The weight
matrices at each layer however either have mutually independent identically
distributed Gaussian entries with mean 0 and variance
$\sigma_{w}^{2}/N_{l-1}$, or are drawn from a uniform distribution over scaled
orthogonal matrices such that
$(\mathbf{W}^{(l)})^{T}\mathbf{W}^{(l)}=\sigma_{w}^{2}\mathbf{I}$. To analyse
the Jacobian of the network at an arbitrary point $\mathbf{x}$, the authors
considered again the large width setting, in which, as per the discussion in
Section 2.1, the preactivations at layer $l$ layer can be modelled as mutually
independent, identically distributed, centred Gaussian random variables with
variance $q^{(l)}$ 111Note that (3) was derived in the context of Gaussian
initialisation, for now we also assume it holds true also for the Orthogonal
case.. For both of the initialisation schemes described, the limiting spectral
density of $\mathbf{J}(\mathbf{x})$, in terms of its moment generating
transform

$M_{JJ^{T}}(z)=\sum_{k=1}^{\infty}\frac{m_{k}}{z^{k}},$ (11)

was computed and analysed using tools from free probability [9, 13, 27]. We
refer the reader to Appendix 6 of [9] for further details. In the context of
ensuring dynamical isometry, then of particular interest are the first and
second moments, which can be expressed in terms of the moments $\mu_{k}$ of
$\mathbf{D}^{2}$ and $s_{k}$ of $(\mathbf{W}^{T}\mathbf{W})$. Here we drop the
dependence on the layer index $l$ by assuming that the condition
$q^{(0)}=q^{*}$ holds true for all inputs $\mathbf{x}$. Adopting again the
notation used in [9],

$\displaystyle\mu_{k}$
$\displaystyle:=\int\phi^{\prime}(\sqrt{q^{*}z})^{2k}d\gamma(z)=\mathbb{E}[\phi^{\prime}(\sqrt{q^{*}})Z)^{2k}],$
(12) $\displaystyle m_{1}$
$\displaystyle:=(\sigma_{w}^{2}\mu_{1})^{L}=(\chi_{1})^{L},$ $\displaystyle
m_{2}$
$\displaystyle:=(\chi_{1})^{2L}L\Big{(}\frac{\mu_{1}}{\mu_{2}}+\frac{1}{L}-1-s_{1}\Big{)}.$

It is evident from (12) that the mean squared singular value $m_{1}$ of the
Jacobian either exponentially explodes or vanishes unless the network is
initialised on the EOC. However, although initialisation on the EOC makes
$m_{1}$ independent of depth, the variance

$\sigma_{{JJ}^{T}}^{2}=m_{2}-m_{1}^{2}=L\Big{(}\frac{\mu_{2}}{\mu_{1}^{2}}-1-s_{1}\Big{)}$
(13)

grows linearly with depth. As a result, under the limiting width assumption
then for an arbitrary, normalised input vector initialisation on the EOC may
still result in an increasingly ill-conditioned Jacobian with depth $L$.
Indeed, we require that $m_{1}=1$ and $\sigma_{{JJ}^{T}}^{2}\approx 0$ so as
to at least approximately achieve dynamical isometry. Equation (13) shows that
using a Gaussian initialisation scheme, corresponding to $s_{1}=-1$, results
in a linear growth in $\sigma_{JJ^{T}}^{2}$ with depth regardless of the
activation function used. As a result, deep, dense feed-forward networks
cannot achieve dynamical isometry using Gaussian initialisation. In contrast,
orthogonal initialisation, corresponding to $s_{1}=0$, can in theory be used
to ensure that $\sigma_{JJ^{T}}^{2}$ is arbitrarily close to one by
controlling the moment ratio $\mu_{2}/\mu_{1}^{2}$. However, the ability to
control the moment ratio depends on the activation function deployed: in the
case of ReLU for instance, whose EOC is a singleton, the moment ratio is a
constant and hence once again dynamical isometry cannot be achieved. As
highlighted in [9], for certain activation functions, e.g., erf and tanh,
which have a non-singleton EOC as illustrated by Figure 2(c), then the moment
ratio can be reduced by shrinking $q^{*}$, which in turn can be achieved by
reducing $\sigma_{b}^{2}$, see Figure 2 (d). Finally we remark that similar
analyses were conducted in [28] in the context of finite width deep ReLU
networks.

###  2.3 Contribution: activations which approximately preserve correlations
and achieve dynamical isometry

To recap, in Section 2.1 we discussed how, in order to achieve deep
information propagation and avoid the network being either highly sensitive or
insensitive to perturbations of the input, it is important to avoid a rapid
rate of convergence of the correlation. To this end it is necessary to choose
$(\phi,\sigma_{w}^{2},\sigma_{b}^{2})$ so that $\chi_{1}=1$ and the associated
correlation function is close to identity. In Section 2.2 we discussed the
problem of vanishing and exploding gradients. Prior works suggest that this
problem can be avoided if the singular values of the input-output Jacobian at
each layer concentrate around 1. The condition $\chi_{1}=1$ is equivalent to
ensuring that the mean of the input-output Jacobian’s spectrum at any layer is
one. Therefore, if the variance $\sigma_{JJ^{T}}^{2}$ of the spectrum, defined
in (13), is close to zero, then the spectrum of the input-output Jacobian is
guaranteed to concentrate around one. Furthermore, with orthogonal
initialisation on the EOC then so long as $\mu_{2}/\mu_{1}^{2}\approx 1$ then
$\sigma_{JJ^{T}}^{2}\approx 0$.

In this paper we present principles for choosing the activation function
$\phi$, see Definition 2.1, and an additional linear scaling parameter which
simultaneously enables uniform convergence of the correlation function
$R_{\phi,q^{*}}(\rho)$ to the identity map, and convergence of the variance
$\sigma_{JJ^{T}}^{2}$ of the input-output Jacobian’s spectrum to zero, without
requiring $\sigma_{b}^{2}$ to shrink towards zero. For further details we
refer the reader to Theorem 2.1. From prior work it was unclear how one could
achieve this. For example, initialising on the EOC, using orthogonal
initialisation and deploying the modulus activation function achieves perfect
dynamical isometry, but results in rapid convergence of correlations in the
forward pass, leading to poor training outcomes. Treated separately, a key
theme that emerges in the prior works towards achieving both goals is that of
reducing $\sigma_{b}^{2}$ in order to make $q^{*}$ small. As mentioned briefly
in section 2.1, in [10] the authors highlight that, for a smooth class of
activation functions, taking this action can reduce the gap between the
correlation map and the identity. The authors also provide a rule for
selecting $\sigma_{b}^{2}$ just small enough so that pairwise correlations
avoid being within some $\epsilon>0$ of one at the output layer. Furthermore,
experiments conducted in [10] indicate that selecting the variance
hyperparameters in this manner accelerates training in practice. However, the
authors also highlight that as the depth of the network increases this rule
still requires $\sigma_{b}^{2}\rightarrow 0$. Likewise, and as mentioned in
section 2.2, in [9] the authors show, for the htanh and erf activation
functions, that reducing $q^{*}$ by shrinking $\sigma_{b}^{2}$ results in the
moment ratio converging towards one. However, this solution is not entirely
satisfactory, first if the expected length of the activations goes to zero
then asymptotically the network will on average annihilate inputs and fail to
pass information to the output. Second, experimental evidence indicates that a
small, but not overly small $\sigma_{b}^{2}$, gives the best results in
practice, which we demonstrate in Section 4. Third, this solution relies only
on the choice of $\sigma_{b}^{2}$ and ignores the role and design of $\phi$.
Our work seeks to address this issue by achieving both goals through the
design of $\phi$, even when $\sigma_{b}^{2}$ and $q^{*}$ are fixed away from
zero.

In light of the above, and assuming $\sigma_{b}^{2}$ is fixed and that
$\sigma_{w}^{2}$ is chosen in order that $\chi_{1}=1$, or goal is to design or
choose $\phi$ in order that both
$\max_{\rho\in[-1,1]}|R_{\phi,q^{*}}(\rho)-\rho|\approx 0$ and
$|\mu_{2}/\mu_{1}^{2}-1|\approx 0$. As highlighted in [10], linear activation
functions are advantageous from the perspective of slowing the convergence of
the pairwise correlations in the forward pass. Likewise, [27, 12] highlight
that linear networks with orthogonal initialisation achieve perfect dynamical
isometry. However, linear activations are not a suitable option as a linear
network can only represent linear functions. To preserve the rich
approximation capabilities of nonlinear networks we analyse a set of
activation functions characterised by being odd, bounded, Lipschitz continuous
and linear around the origin. We refer to activation functions of this type,
defined in Definition 2.1, as scaled-bounded activations.

######  Definition 2.1 (scaled-bounded activations).

We refer to the set of activation functions
$\phi:\mathbb{R}\rightarrow\mathbb{R}$ which satisfy the following properties
as scaled-bounded activations.

  1. 1.

Continuous.

  2. 2.

Odd, meaning that $\phi(z)=-\phi(-z)$ for all $z\in\mathbb{R}$.

  3. 3.

Linear around the origin and bounded: in particular there exists
$a,k\in\mathbb{R}_{>0}$ such that $\phi(z)=kz$ for all $z\in[-a,a]$ and
$\phi(z)\leq ak$ for all $z\in\mathbb{R}$.

  4. 4.

Twice differentiable at all points $z\in\mathbb{R}\backslash\mathcal{D}$,
where $\mathcal{D}\subset\mathbb{R}$ is a finite set. Furthermore
$|\phi^{\prime}(z)|\leq k$ for all $z\in\mathbb{R}\backslash\mathcal{D}$.

As illustrated in Figure 3, the structure of scaled-bounded activation
functions outside of their linear region $[-a,a]$ can vary substantially. We
note that for any scaled-bounded activation $\phi$ there exists
$a,k\in\mathbb{R}_{>0}$ such that $|\phi(z)|\leq|\text{Shtanh}_{a,k}(z)|$,
where

$\text{Shtanh}_{a,k}(z)=\begin{cases}zk,\;\;|z|<a\\\ ak,\;\;|z|\geq
a.\end{cases}$ (14)

Here Shtanh stands for scaled-bounded hard tanh, and we further note that
$\text{Shtanh}_{a,k}(\cdot)$ is a scaled-bounded activation for any
$a,k\in\mathbb{R}_{>0}$. Our main contribution is Theorem 2.1, which asserts,
for scaled-bounded activations, that the correlation function $R_{\phi,q^{*}}$
can be made arbitrarily close to the identity on $[0,1]$, and the moment ratio
$\mu_{2}/\mu_{1}^{2}$ arbitrarily close to one, by choosing
$\sigma_{b}^{2}/a^{2}$ to be small.

Figure 3: examples of scaled-bounded activation functions with $k=1$ and
$a\in\\{1,2,4\\}$. The prefix S refers to the scaling of the linear region via
the parameter $a$.

######  Theorem 2.1.

Let $\phi$ be a scaled-bounded activation, see Definition 2.1,
$\sigma_{b}^{2}>0$ and suppose that

$\chi_{1}\vcentcolon=\sigma_{w}^{2}\mathbb{E}[\phi^{\prime}(\sqrt{q^{*}}Z)^{2}]=1,$

where $q^{*}>0$ is a fixed point of the associated variance function
$V_{\phi}$. In addition, assume that all inputs $\mathbf{x}$ are normalised so
that $||\mathbf{x}||_{2}^{2}=q^{*}$. With
$y\vcentcolon=\frac{\sigma_{b}^{2}}{a^{2}}$ and $\Lambda$ defined as in Lemma
3.4, then

$\max_{\rho\in[0,1]}|R_{\phi}(\rho)-\rho|<\left(\frac{8}{\pi}\right)^{1/3}y^{1/3}$
(15)

and

$\left|\frac{\mu_{2}}{\mu_{1}^{2}}-1\right|\leq\text{erf}\left(\frac{\Lambda(y)}{\sqrt{2}}\right)^{-2}-1.$
(16)

We emphasise that as $y\vcentcolon=\sigma_{b}^{2}/a^{2}\rightarrow 0$ then
both

$\max_{\rho\in[0,1]}|R_{\phi,q^{*}}(\rho)-\rho|,\left|\mu_{2}/\mu_{1}^{2}-1\right|\rightarrow
0.$

We first remark that the uniform bound on the correlation provided by Theorem
2.1 is only for nonnegative correlations. However, Figure 6 indicates that
letting $\sigma_{b}^{2}/a^{2}\rightarrow 0$ results
$\max_{\rho\in[-1,1]}|R_{\phi,q^{*}}(\rho)-\rho|\rightarrow 0$. The key
takeaway of Theorem 2.1 for practitioners is that an improved initialisation
can be achieved by using an activation function which has a sufficiently large
linear region spanning either side of the origin. The size of this linear
region, $a$, needs to be selected on the basis of both the bias variance
hyperparameter $\sigma_{b}^{2}$ and the depth of the network. The deeper the
network or the larger $\sigma_{b}^{2}$ is, the larger $a$ needs to be in order
to avoid both convergence of the correlations and achieve approximate
dynamical isometry. In practice, as per Figure 5, modest increases in $a$
typically suffice to capture many of potential benefits at initialisation. We
emphasise that there is a tension between the linear region being large enough
to achieve a good initialisation, while being small enough so that the
expressivity of the network is not reduced. In particular, if the network is
initialised so that it acts as linear transform on the input data, and if the
optimiser becomes trapped in a local minimum close to the initialisation
point, then the network may continue to act as a linear transform throughout
training.

##  3 Analysis of scaled-bounded activation functions

###  3.1 Derivation of Theorem 2.1

We start our analysis by revisiting what it means to initialise on the EOC:
the condition in (7) presupposes the existence of a fixed point $q^{*}$ of
$V_{\phi}$. To resolve this matter we introduce and study the fixed points of
a related function, $W$. Before presenting this analysis we introduce the
following slight abuse of notation,

$\phi^{\prime}(z):=\frac{1}{2}\left(\frac{d\phi}{dz}(z^{-})+\frac{d\phi}{dz}(z^{+})\right)$
(17)

where $\frac{d\phi}{dz}(z^{-})$ and $\frac{d\phi}{dz}(z^{+})$ are the left and
right derivatives of $\phi$ at $z\in\mathbb{R}$ respectively. This is
convenient for what follows as it allows us to define a notion of a derivative
of a scaled-bounded activation $\phi$ over the whole of the domain: we remark,
by Definition 2.1, that the true derivative of $\phi$, and
$\phi^{\prime}(\cdot)$, as defined above, are equal almost everywhere.

######  Lemma 3.1.

Let $\phi$ be a scaled-bounded activation, see Definition 2.1, and
$\sigma_{b}^{2}>0$. Define

$W_{\phi}(q):=\frac{\mathbb{E}[\phi(\sqrt{q}Z)^{2}]}{\mathbb{E}[\phi^{\prime}(\sqrt{q}Z)^{2}]}+\sigma_{b}^{2}$
(18)

for all $q\in\mathbb{R}_{\geq 0}$. Then $W_{\phi}:\mathbb{R}_{\geq
0}\rightarrow\mathbb{R}_{\geq 0}$ and $W$ has a fixed point $q^{*}>0$.

A proof of Lemma 3.1 is provided in Appendix B.1. As a consequence of Lemma
3.1, we are able to make the following claims concerning the existence of
fixed points of the variance function $V_{\phi}$ for any scaled-bounded
activation $\phi$.

######  Corollary 3.1.1.

Let $\phi$ be a scaled-bounded activation, see Definition 2.1,
$\sigma_{b}^{2}>0$ and suppose

$\chi_{1}\vcentcolon=\sigma_{w}^{2}\mathbb{E}[\phi^{\prime}(\sqrt{q^{*}}Z)^{2}]=1,$
(19)

where $q^{*}>0$ is a fixed point of $W_{\phi}$, defined in (18). Then $q^{*}$
is a fixed point of the associated variance function $V_{\phi}$, defined in
(3).

A proof of corollary 3.1.1 is included in Appendix B.2 for completeness. The
key takeaway of Corollary 3.1.1 is that for any scaled-bounded activation
there exists a fixed point $q^{*}$ of $V_{\phi}$ satisfying $q^{*}>0$. In
fact, as $V_{\phi}(0)=\sigma_{b}^{2}>0$ then any fixed point of $V_{\phi}$ is
greater than 0. We emphasise that such analysis is necessary as in general
care is required when making assumptions concerning the existence of fixed
points. For example, the variance function for ReLU like functions do not have
any fixed points unless $\sigma_{b}^{2}=0$.

To simplify the analysis we will now proceed under the assumption that the
input data is normalised to have euclidean length $q^{*}$ and that
$q^{(l)}=q^{*}$ for all $l\in[L]$. If one where interested in initialising in
this manner in practice, it might also be of interest to explore guarantees
concerning the ease with which a stable fixed point of $V_{\phi}$ can be
computed: in particular if the fixed point is not attractive then numerical
precision issues might mean that the sequence $q^{(l)}$ diverges from $q^{*}$.
We do not pursue this line of inquiry and instead leave it as potential future
work. Now that we have identified the existence of fixed points of the
variance functions of scaled-bounded activations, it is possible to study the
associated correlation functions.

######  Lemma 3.2.

Let $\phi$ be a scaled-bounded activation, see Definition 2.1,
$\sigma_{b}^{2}>0$ and suppose

$\chi_{1}\vcentcolon=\sigma_{w}^{2}\mathbb{E}[\phi^{\prime}(\sqrt{q^{*}}Z)^{2}]=1,$
(20)

where $q^{*}>0$ is a fixed point of $W_{\phi}$. In addition, assume that all
inputs $\mathbf{x}$ are normalised so that $||\mathbf{x}||_{2}^{2}=q^{*}$.
Then the associated correlation map $R_{\phi,q*}$, defined in (5), is fixed at
each layer $l\in[L]$, satisfies $R_{\phi,q*}:[-1,1]\rightarrow[-1,1]$ and is
differentiable with

$R_{\phi,q^{*}}^{\prime}(\rho)=\sigma_{w}^{2}\mathbb{E}[\phi^{\prime}(U_{1})\phi^{\prime}(U_{2})]$
(21)

for all input correlations $\rho\in(-1,1)$.

A proof of Lemma 3.2 is given in Appendix B.3. This expression for the
correlation is equivalent to that given in [7], however we emphasise, due to
the fact that $\phi$ is not necessarily continuously differentiable, that
significantly more care is required to derive it. Indeed, if $\phi$ were not
odd then this equivalence would not hold and additional terms would appear,
potentially complicating the analysis downstream. We are now ready to present
a key lemma, which provides a tight uniform bound on the interval $[0,1]$
between the correlation functions associated with scaled-bounded activations
and the identity function.

######  Lemma 3.3.

Under the same conditions and assumptions as in Lemma 3.2, it holds that

$\max_{\rho\in[0,1]}|R_{\phi,q^{*}}(\rho)-\rho|=\frac{\sigma_{b}^{2}}{q^{*}}.$
(22)

A proof of Lemma 3.3 is provided in Appendix B.4. The final lemma we present
before proving Theorem 2.1 concerns the relationship between $a$, the size of
the linear region, and $q^{*}$, the fixed point of the variance function. This
lemma provides lower and upper bounds on the ratio $\frac{a}{\sqrt{q^{*}}}$
for all scaled-bounded activations as a function of $a$ and $\sigma_{b}^{2}$.
The fact that this is possible indicates that the particular shape of the
tails of a scaled-bounded activation do not play a key role in achieving a
good initialisation. In Figure 4 we plot these bounds, as well as example
$a/\sqrt{q^{*}}$ ratios, for a number of activation functions as a function of
$a$.

Figure 4: Bounds on $\frac{a}{\sqrt{q^{*}}}$ given in (23) vs.
$\frac{a}{\sqrt{q^{*}}}$, computed numerically for the same scaled-bounded
activations visualised in Figure 3.

######  Lemma 3.4.

Under the same conditions and assumptions as in Lemma 3.2, and defining
$y:=\frac{\sigma_{b}^{2}}{a^{2}}$, then

$\Lambda(y)<\frac{a}{\sqrt{q^{*}}}<\left(\frac{8}{\pi}\right)^{1/6}y^{-1/3},$
(23)

where $\Lambda(y)$ is defined as

$\left(W_{0}\left(\frac{2}{\pi}\left(\sqrt{\frac{8}{\pi}}\exp\left(-\frac{\left(\sqrt{W_{0}\left(\frac{2}{\pi}y^{-2}\right)}\right)^{2}}{2}\right)\left(\frac{1}{\left(\sqrt{W_{0}\left(\frac{2}{\pi}y^{-2}\right)}\right)}+\left(\sqrt{W_{0}\left(\frac{2}{\pi}y^{-2}\right)}\right)\right)\right)^{-2}\right)\right)^{1/2}$

and $W_{0}$ denotes the principal branch of the Lambert W function.

A proof of Lemma 3.4 is given in Appendix B.5. We note that while the upper
bound in Lemma 3.4 is easy to interpret, the lower bound is not immediately
interpretable. However, this lower bound still allows us to compute a
numerical lower bound for any scaled-bounded activation as per Figure 4. In
terms of asymptotic behaviour it is easy to check that as $y\rightarrow 0$
then $y^{-1/3},\Lambda(y)\rightarrow\infty$, and as $y\rightarrow\infty$ then
$y^{-1/3},\Lambda(y)\rightarrow 0$. We also observe from Figure 4 that, at
least empirically, the lower bound seems to be far tighter: we leave it as
potential future work to investigate deriving a tighter upper bound. We are
now ready to prove Theorem 2.1.

###### Proof.

To derive the inequality concerning the correlation function given (15) we use
the upper bound

$\frac{a}{\sqrt{q^{*}}}<\left(\frac{8}{\pi}\right)^{1/6}\frac{a^{2/3}}{\sigma_{b}^{2/3}}$

derived in Lemma 3.4. Squaring both sides, dividing by $a^{2}$ and multiplying
by $\sigma_{b}^{2}$ then

$\frac{\sigma_{b}^{2}}{q^{*}}<\left(\frac{8}{\pi}\right)^{1/3}\frac{\sigma_{b}^{2/3}}{a^{2/3}}.$

To conclude we apply Lemma 3.3. We now turn our attention to proving the
inequality concerning the moment ratio provided in equation 16. Analysing the
$l$th moment,

$\displaystyle\mu_{l}$
$\displaystyle:=\mathbb{E}[\phi^{\prime}(\sqrt{q^{*}}Z)^{2l}]$
$\displaystyle=2\left(k^{2l}\int_{0}^{a/\sqrt{q^{*}}}d\gamma(z)+\int_{a/\sqrt{q^{*}}}^{\infty}\phi^{\prime}(\sqrt{q^{*}}z)^{2l}d\gamma(z)\right).$

By assumption $|\phi^{\prime}(z)|\leq k$, hence we can bound this quantity as

$k^{2l}\text{erf}\left(\frac{a}{\sqrt{2q^{*}}}\right)\leq\mu_{l}\leq k^{2l},$

where the lower bound arises simply by zeroing the second integral term. It
therefore follows that

$\displaystyle\frac{k^{4}\text{erf}\left(\frac{a}{\sqrt{2q^{*}}}\right)}{k^{4}}$
$\displaystyle\leq\frac{\mu_{2}}{\mu_{1}^{2}}$
$\displaystyle\leq\frac{k^{4}}{k^{4}\text{erf}\left(\frac{a}{\sqrt{2q^{*}}}\right)^{2}},$
$\displaystyle\text{erf}\left(\frac{a}{\sqrt{2q^{*}}}\right)$
$\displaystyle\leq\frac{\mu_{2}}{\mu_{1}^{2}}$
$\displaystyle\leq\text{erf}\left(\frac{a}{\sqrt{2q^{*}}}\right)^{-2}$

Observe that
$\mu_{2}/\mu_{1}^{2}-1\leq\text{erf}\left(\frac{a}{\sqrt{2q^{*}}}\right)^{-2}-1$
and $1-\mu_{2}/\mu_{1}^{2}\leq
1-\text{erf}\left(\frac{a}{\sqrt{2q^{*}}}\right)$. We now prove, by
contradiction, that $\text{erf}\left(\frac{a}{\sqrt{2q^{*}}}\right)^{-2}-1\geq
1-\text{erf}\left(\frac{a}{\sqrt{2q^{*}}}\right)$. Indeed, assume that
$\text{erf}\left(\frac{a}{\sqrt{2q^{*}}}\right)^{-2}-1<1-\text{erf}\left(\frac{a}{\sqrt{2q^{*}}}\right)$,
then

$\displaystyle 1-\text{erf}\left(\frac{a}{\sqrt{2q^{*}}}\right)^{2}$
$\displaystyle<\text{erf}\left(\frac{a}{\sqrt{2q^{*}}}\right)\left(\text{erf}\left(\frac{a}{\sqrt{2q^{*}}}\right)-\text{erf}\left(\frac{a}{\sqrt{2q^{*}}}\right)^{2}\right)$
$\displaystyle\leq\text{erf}\left(\frac{a}{\sqrt{2q^{*}}}\right)\left(1-\text{erf}\left(\frac{a}{\sqrt{2q^{*}}}\right)^{2}\right).$

As $\text{erf}\left(\frac{a}{\sqrt{2q^{*}}}\right)<1$ for
$\frac{a}{\sqrt{2q^{*}}}<\infty$ then this is a contradiction. Therefore

$\left|\frac{\mu_{2}}{\mu_{1}^{2}}-1\right|\leq\text{erf}\left(\frac{a}{\sqrt{2q^{*}}}\right)^{-2}-1.$

As erf is monotonically increasing it suffices to lower bound the quantity
$\frac{a}{\sqrt{q^{*}}}$. The result claimed is then recovered by applying the
lower bound on $\frac{a}{\sqrt{q^{*}}}$ derived in Lemma 3.4. ∎

Figure 5: The left and right plots display the correlation and moment ratio
bounds, given in Equations (15) and (16) respectively, vs. the equivalent
numerically computed quantities for a variety of different activation
functions $\phi\in\Omega$, as shown in Figure 3. We note that the right hand
curve for hard saw is accurate only up to an $a$ of 6 due to numerical
precision issues arising from numerical integration steps.

Figure 6: impact of scale parameter $a$ on correlation function of scaled (or
adapted) htanh (left-hand plot) and sinusoid (right-hand plot).

###  3.2 Discussion and practical takeaways

Equation (15) implies, by choosing $a$ sufficiently large relative to
$\sigma_{b}^{2}$, that problems arising as a result of limits on the depth of
information propagation as well as network sensitivity can be mitigated
without the need for $q^{*},\sigma_{b}^{2}\rightarrow 0$. Equation (16)
likewise implies, as long as orthogonal initialisation is used, that
increasing $a$ also mitigates the problem of model degeneracy, again without
the need for $q^{*},\sigma_{b}^{2}\rightarrow 0$. Our numerics support these
conclusions: in Figure 5 the bounds in (15) and (16) are plotted numerically
with $\sigma_{b}^{2}$ fixed, and converge to $0$ as $a$ increases. Likewise,
Figure 6 shows that increasing $a$ moves $R_{\phi,q^{*}}$ closer to the
identity. We remark that (3) and (5) where derived in the context of Gaussian
initialisation. As a result one might question whether (15) is immediately
applicable to orthogonally initialised networks. Indeed, although referenced
to in [9, 13], the correspondence between orthogonally initialised networks
and Gaussian processes is to our knowledge yet to be rigorously established.
However, it seems highly likely that the same correspondence holds, an
assertion supported both by empirical observation and the fact that large
random orthogonal matrices are well approximated by Gaussian matrices (see
e.g., [29]). We defer a detailed study of this correspondence to later work.
We also note that the uniform bound we provide holds for non-negative
correlations only, this is due to the fact that the upper bound on this
interval, $\sigma_{b}^{2}/q^{*}$ as per Lemma 3.3, is relatively easy to
analyse. However, we suspect for scaled-bounded activations that letting
$\sigma_{b}^{2}/a^{2}\rightarrow 0$ results in a uniform convergence of the
correlation function to the identity for $\rho\in[-1,1]$. This hypothesis is
supported by Figure 6.

We conclude this section by considering the relevance and importance of each
of the four properties of scaled-bounded activation functions in regard to
achieving a good initialisation in practice. Property 1), $\phi$ being
continuous, does not seem controversial, indeed most of the commonly used
activation functions are continuous. Furthermore, it seems reasonable that
continuous activation functions result in a loss landscape that is easier for
an optimiser to navigate, compared with that induced by non-continuous ones.
Property 2), $\phi$ being odd, is beneficial in theory as it removes
additional constant terms from (22) in Lemma 3.3. However, as per the proof of
Lemma 3.3, it can be observed that these terms will decay exponentially fast
as $a$ grows due to the fact that the locations of the points in $\mathcal{D}$
all have magnitude at least $a$. As for property 3), while the existence of a
linear region around the origin is critical to our results, we suspect that
boundedness is more an artefact of our proof than a necessity in practice.
Assuming that $\phi$ is bounded simplifies certain pieces of analysis, and
crucially allows us to formulate upper and lower bounds on $V_{\phi}(q)$,
needed for the proof of Lemma 3.4. We hypothesise that this condition could be
relaxed to $|\phi(z)|<|z|$. In regard to property 4) and contrary to the
conclusion one might be inclined to draw from [10], $\phi$ being non-
differentiable demonstrates that smoothness is not necessary for a good
initialisation. Finally, the bound on the derivative of $\phi$ was introduced
to allow us to to derive bounds on $V_{\phi}$.

##  4 Experiments

###  4.1 Experimental setup

Across all experiments we train networks with depths
$L\in\\{20,50,100,200\\}$, with each layer having a fixed width $N=400$, for
$100$ epochs on CIFAR-10. The results are averaged over 10 trials. Variance
hyperparameters $(\sigma_{w}^{2},\sigma_{b}^{2})$ are selected to lie on the
EOC, with $\sigma_{b}^{2}\in\\{1,10^{-1},10^{-2},10^{-3},10^{-4}\\}$ and
$\sigma_{w}^{2}$ computed in in order that $\chi_{1}=1$, see (7). The
parameters of the network are initialised according to the following two
schemes.

  * •

Gaussian initialisation: all parameters are drawn mutually independent of one
another. The weights in each layer are identically distributed with
$w_{i,j}^{(l)}\sim\mathcal{N}(0,\sigma_{w}^{2}/N)$. The biases are all
identically distributed with $b_{i,j}^{(l)}\sim\mathcal{N}(0,\sigma_{b}^{2})$.

  * •

Orthogonal initialisation: all weight matrices and bias parameters are drawn
mutually independent of one another. The weight matrix at each layer is drawn
according to the Haar measure, i.e., uniformly, over the orthogonal group of
$N\times N$ matrices such that
$(\mathbf{W}^{(l)})^{T}\mathbf{W}^{(l)}=\sigma_{w}^{2}\textbf{I}_{N}$. The
biases are all identically distributed with
$b_{i,j}^{(l)}\sim\mathcal{N}(0,\sigma_{b}^{2})$.

Finally, optimisation is performed using SGD with a batch size of 64. A
learning rate of 10-4 was found to be provide good results across all
experiments.

###  4.2 Experimental validation of the results of Section 3

In order to test and validate the results of Section 3, we deploy the scaled
Shtanh activation function, defined in (14), with $a\in\\{1,2,5,10\\}$ and
$k=1$, and measure the test accuracy of the networks described in Section 4.1
at various stages of training. The results of these experiments are summarised
in the heatmaps provided in Figure 7. As per the implications of Theorem 2.1,
it is clear from Figure 7 that as $\sigma_{b}^{2}$ increases a larger but not
necessarily maximal value of $a$ gives the best results. This is illustrated
by the ridges, indicating the highest test accuracy, running from top left to
bottom right in the heatmap subplots. This therefore highlights the trade-off
between activation functions which are linear enough to allow for a good
initialisation, while not being overly linear that the network loses
approximation power. We emphasise that these results are not unique to Shtanh
and that the same conclusions can be drawn for other scaled-bounded activation
functions 222Implementation and additional results can be found at
https://github.com/Cross-Caps/AFLI. As established in prior works, we also
observe an advantage in using orthogonal over Gaussian initialisation.
However, we also note that this performance gap is relatively small. We leave
it to future work to investigate the impact of optimising over, or
regularising with respect to, the set of orthogonal weight matrices so that
orthogonality is preserved, at least approximately, throughout training.

(a) Shtanh with Gaussian initialisation (b) Shtanh with orthogonal
initialisation Figure 7: test accuracy on CIFAR-10 at different stages of
training. Note that all sub-heatmap plots share the same colour scale.

We also provide some very preliminary results concerning a new training
protocol, in which a scaled-bounded activation function is deployed and the
value of $a$ is decreased slowly over the first few epochs. To motivate this,
observe in Figure 8 that, as per our theory, a large value of $a$ always gives
the best performance at initialisation. However, for reasons already
discussed, an overly large value of $a$ has a negative impact on training long
term. This strategy can therefore be interpreted as computing a good
initialisation for the final network. The purple curve in Figure 8 displays
the outcome of this strategy, which appears promising over a wide range of
depths and hyperparameter choices.

(a) Shtanh with Gaussian initialisation (b) Shtanh with orthogonal
initialisation Figure 8: test accuracy during training on CIFAR-10. Networks
are trained either with a fixed value of $a$ or with $a$ linearly decreasing
from 10$\rightarrow$2 over the first 30 epochs.

##  5 Conclusion and avenues for future work

In this paper we considered the role of the activation function in avoiding
certain problems at initialisation, namely limited information propagation
with depth, high sensitivity or insensitivity to perturbations of the input,
and vanishing and exploding gradients. The first two of these we investigated
by studying the dynamics of the preactivation correlations, and the third by
analysing the spectrum of the input-output Jacobian of the network.
Previously, these problems have been ameliorated by shrinking $\sigma_{b}^{2}$
in relation to the depth of the network. This is unsatisfactory for two
reasons: firstly it results in the expected euclidean length of the
activations shrinking towards 0 with depth, and second it places constraints
on the initialisation regime, which can result in suboptimal training
outcomes. Our theory and experiments clarify that shrinking $\sigma_{b}^{2}$
is not necessary, instead, it is possible to avoid these problems at
initialisation by ensuring that the activation function deployed has a
sufficiently large linear region around the origin. This work therefore
provides a rigorous explanation for the observation that activation functions
which approximate the identity near the origin perform well, particularly at
initialisation.

Avenues for future work include characterising more precisely how large the
linear region needs to be for a given depth, a more comprehensive
investigation as to the potential benefits of a more affine initialisation,
and an analysis of the approximation capabilities of a neural network whose
parameters are constrained to lie in some neighbourhood of a given
initialisation point. In addition, a quantitative description of how the
preactivation correlation dynamics impact training is desirable.

## Acknowledgements

This work is supported by the Alan Turing Institute under the EPSRC grant
EP/N510129/1 and the Ana Leaf Foundation.

## References

  * [1] Xavier Glorot and Yoshua Bengio.  Understanding the difficulty of training deep feedforward neural networks.  In International Conference on Artificial Intelligence and Statistics, pages 249–256, May 2010.
  * [2] Prajit Ramachandran, Barret Zoph, and Quoc V. Le.  Searching for activation functions.  In 6th International Conference on Learning Representations, ICLR 2018, Workshop Track Proceedings, 2018.
  * [3] Ludovic Trottier, Philippe Giguère, and Brahim Chaib-draa.  Parametric exponential linear unit for deep convolutional neural networks.  In 16th IEEE International Conference on Machine Learning and Applications, ICMLA 2017, pages 207–214, 2017.
  * [4] Xinyu Liu and Xiaoguang Di.  Tanhexp: A smooth activation function with high convergence speed for lightweight neural networks.  IET Computer Vision, 15(2):136–150, 2021.
  * [5] Günter Klambauer, Thomas Unterthiner, Andreas Mayr, and Sepp Hochreiter.  Self-normalizing neural networks.  In Advances in Neural Information Processing Systems 30, pages 971–980, December 2017.
  * [6] Djork-Arné Clevert, Thomas Unterthiner, and Sepp Hochreiter.  Fast and accurate deep network learning by exponential linear units (ELUs).  International Conference on Learning Representations (ICLR), May 2015.
  * [7] Ben Poole, Subhaneil Lahiri, Maithra Raghu, Jascha Sohl-Dickstein, and Surya Ganguli.  Exponential expressivity in deep neural networks through transient chaos.  In International Conference on Neural Information Processing Systems (NIPS), pages 3368–3376, 2016.
  * [8] Samuel S. Schoenholz, Justin Gilmer, Surya Ganguli, and Jascha Sohl-Dickstein.  Deep information propagation.  In International Conference on Learning Representations (ICLR), April 2017.
  * [9] Jeffrey Pennington, Samuel Schoenholz, and Surya Ganguli.  The emergence of spectral universality in deep networks.  In International Conference on Artificial Intelligence and Statistics (AISTATS), pages 1924–1932, April 2018.
  * [10] Soufiane Hayou, Arnaud Doucet, and Judith Rousseau.  On the impact of the activation function on deep neural networks training.  In International Conference on Machine Learning (ICML), pages 2672–2680, June 2019.
  * [11] Sepp Hochreiter.  The vanishing gradient problem during learning recurrent neural nets and problem solutions.  International Journal of Uncertainty, Fuzziness and Knowledge-Based Systems, 6(2):107–116, 1998.
  * [12] Andrew M. Saxe, James L. McClelland, and Surya Ganguli.  Exact solutions to the nonlinear dynamics of learning in deep linear neural networks.  In International Conference on Learning Representations(ICLR), April 2014.
  * [13] Jeffrey Pennington, Samuel S. Schoenholz, and Surya Ganguli.  Resurrecting the sigmoid in deep learning through dynamical isometry: Theory and practice.  In International Conference on Neural Information Processing Systems (NeurIPS), page 4788–4798, December 2017.
  * [14] Sepp Hochreiter and Jürgen Schmidhuber.  Long short-term memory.  Neural Computing, 9(8):1735–1780, November 1997.
  * [15] K. He, X. Zhang, S. Ren, and J. Sun.  Deep residual learning for image recognition.  In International Conference on Computer Vision and Pattern Recognition (CVPR), pages 770–778, June 2016.
  * [16] Rupesh . Srivastava, Klaus Greff, and Jürgen Schmidhuber.  Training very deep networks.  In Advances in Neural Information Processing Systems (NeurIPS), pages 2377–2385, 2015.
  * [17] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Delving deep into rectifiers: Surpassing human-level performance on imagenet classification.  In IEEE International Conference on Computer Vision (ICCV), pages 1026–1034, 2015.
  * [18] Sergey Ioffe and Christian Szegedy.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  In International Conference on International Conference on Machine Learning (ICML), pages 448–456, 2015.
  * [19] Jimmy Lei Ba, Jamie Ryan Kiros, and Geoffrey E Hinton.  Layer normalization.  arXiv preprint arXiv:1607.06450, 2016.
  * [20] Hamed Aghdam and Elnaz Heravi.  Guide to Convolutional Neural Networks: A Practical Application to Traffic-Sign Detection and Classification.  Springer, January 2017.
  * [21] Radford M. Neal.  Bayesian Learning for Neural Networks.  Springer-Verlag, Berlin, Heidelberg, 1996.
  * [22] Alexander G. de G. Matthews, Jiri Hron, Mark Rowland, Richard E. Turner, and Zoubin Ghahramani.  Gaussian process behaviour in wide deep neural networks.  In International Conference on Learning Representations, 2018.
  * [23] Jaehoon Lee, Yasaman Bahri, Roman Novak, Samuel S. Schoenholz, Jeffrey Pennington, and Jascha Sohl-Dickstein.  Deep neural networks as gaussian processes.  In International Conference on Learning Representations (ICLR), pages 1–11, April 2018.
  * [24] Boris Hanin and David Rolnick.  How to start training: The effect of initialization and architecture.  In Advances in Neural Information Processing Systems (NeurIPS), pages 571–581, December 2018.
  * [25] Eldad Haber and Lars Ruthotto.  Stable architectures for deep neural networks.  Inverse Problems, 34(1):014004, Dec 2017.
  * [26] David E. Rumelhart, Geoffrey E. Hinton, and Ronald J. Williams.  Learning representations by back-propagating errors.  Nature, 323(6088):533–536, 1986.
  * [27] Jeffrey Pennington and Yasaman Bahri.  Geometry of neural network loss surfaces via random matrix theory.  In International Conference on Machine Learning (ICML), pages 2798–2806, 2017.
  * [28] Boris Hanin.  Which neural net architectures give rise to exploding and vanishing gradients?  In Advances in Neural Information Processing Systems, volume 31, 2018.
  * [29] Elizabeth S. Meckes.  The Random Matrix Theory of the Classical Compact Groups.  Cambridge Tracts in Mathematics. Cambridge University Press, 2019.
  * [30] G. B. Folland.  Real analysis : modern techniques and their applications.  Wiley, New York, 1999.
  * [31] G. K. Karagiannidis and A. S. Lioumpas.  An improved approximation for the gaussian q-function.  IEEE Communications Letters, 11(8):644–646, 2007.
  * [32] R. M. Corless, G. H. Gonnet, D. E. G. Hare, D. J. Jeffrey, and D. E. Knuth.  On the lambertw function.  Advances in Computational Mathematics, 5(1):329–359, 1996.

##  Appendix A Supporting Lemmas

For the sake of clarity and completeness we recall here two well known Lemmas
concerning Lebesgue integrals. Both can be proved using a combination of
Lebesgue’s dominated convergence theorem and the mean value theorem (See e.g.,
Chapter 2 of [30]).

######  Lemma A.1.

Let $(X,\mathcal{F},\mu)$ be a measure space and $Y\subset\mathbb{R}$ an open
interval. Consider a function $f:X\times Y\rightarrow\mathbb{R}$ such that the
following are true.

  1. 1.

For each $y\in Y$ then the function $f_{y}:X\rightarrow\mathbb{R}$ with
$f_{y}(x)\vcentcolon=f(x,y)$ satisfies $f_{y}(x)\in L^{1}(X,\mathcal{F},\mu)$.

  2. 2.

For each $y\in Y$ then $\lim_{y^{\prime}\rightarrow y}f(x,y^{\prime})=f(x,y)$
almost everywhere in $X$.

  3. 3.

For each $y\in Y$ there exists an open interval $K$, with $y\in K$, and a
$g_{K}(x)\in L^{1}(X,\mathcal{F},\mu)$ such that

$|f(x,y^{\prime})|\leq g_{K}(x)$

for all $y^{\prime}\in K$.

Then the function $F:Y\rightarrow\mathbb{R}$ with $F(y):=\int_{X}f(x,y)d\mu$
is continuous on $Y$.

######  Lemma A.2.

Let $(X,\mathcal{F},\mu)$ be a measure space and $Y\subset\mathbb{R}$ an open
interval. Consider a function $f:X\times Y\rightarrow\mathbb{R}$ such that the
following are true.

  1. 1.

For each $y\in Y$ then the function $f_{y}:X\rightarrow\mathbb{R}$ with
$f_{y}(x)\vcentcolon=f(x,y)$ satisfies $f_{y}(x)\in L^{1}(X,\mathcal{F},\mu)$.

  2. 2.

For each $y\in Y$ then $\frac{\partial f}{\partial y}(x,y)$ exists almost
everywhere in $X$.

  3. 3.

For each $y\in Y$ there exists an open interval $K$, with $y\in K$, and a
$g_{K}(x)\in L^{1}(X,\mathcal{F},\mu)$ such that

$|\frac{\partial f}{\partial y}(x,y^{\prime})|\leq g_{K}(x)$

for all $y^{\prime}\in K$.

Then the function $F:Y\rightarrow\mathbb{R}$ with $F(y):=\int_{X}f(x,y)dx$ is
differentiable on $Y$ with

$F^{\prime}(y)=\int_{X}\frac{\partial f}{\partial y}(x,y)d\mu.$ (24)

We also now present a specific adaptation of integration by parts for
piecewise continuously differentiable functions of Gaussian random variables
with a finite number of discontinuities and bounded derivative.

######  Lemma A.3.

Suppose $f:\mathbb{R}\rightarrow\mathbb{R}$ is bounded, piecewise continuously
differentiable at all but a finite number $T\in\mathbb{N}$ of non-
differentiable points, $t_{1}<t_{2}<...<t_{T}$, and has bounded derivative.
Then

$\int_{\mathbb{R}}zf(z)d\gamma(z)=\sum_{i}^{T}\left[-\frac{\exp(-\frac{1}{2}z^{2})}{\sqrt{2\pi}}f(z)\right]_{t_{i}^{+}}^{t_{i}^{-}}+\int_{\mathbb{R}}f^{\prime}(z)d\gamma(z)$
(25)

###### Proof.

We first note that as $f$ and $f^{\prime}$ are continuous and continuous
almost everywhere but at a finite number of points respectively, then they are
clearly measureable with respect to the completion of the Borel sigma algebra
on $\mathbb{R}$. Additionally, under the assumptions that $f$ and $f^{\prime}$
are bounded, it follows that $f(z),zf(z),f^{\prime}(z)\in
L^{1}(\mathbb{R},\mathcal{B}(\mathbb{R}),\gamma)$ where $\gamma$ is the
standard one dimensional Gaussian measure. Defining, for typographical ease,
$g(z)\vcentcolon=zf(z)$, $g_{-}(z)\vcentcolon=\max\\{-zf(z),0\\}$ and
$g_{+}(z)\vcentcolon=\max\\{zf(z),0\\}$, then as $|g_{\pm}(z)|\leq|g(z)|$ it
follows that $g_{-}(z),g_{+}(z)\in
L^{1}(\mathbb{R},\mathcal{B}(\mathbb{R}),\gamma)$. As both $g_{-}$ and $g_{+}$
are nonnegative functions then $(g_{-}\textbf{1}_{[-n,n]})_{n\in\mathbb{N}}$
and $(g_{+}\textbf{1}_{[-n,n]})_{n\in\mathbb{N}}$ are sequences of non-
decreasing functions converging to $g_{-}$ and $g_{+}$ respectively.
Therefore, by monotone convergence

$\displaystyle\int_{\mathbb{R}}zf(z)d\gamma(z)$
$\displaystyle=\int_{\mathbb{R}}g_{+}(z)d\gamma(z)-\int_{\mathbb{R}}g_{-}(z)d\gamma(z)$
$\displaystyle=\lim_{n\rightarrow\infty}\int_{-n}^{n}g_{+}(z)d\gamma(z)-\lim_{n\rightarrow\infty}\int_{-n}^{n}g_{-}(z)d\gamma(z)$
$\displaystyle=\lim_{n\rightarrow\infty}\int_{-n}^{n}zf(z)d\gamma(z).$

Let $n\in\mathbb{N}$ be any integer such that $n>|t_{T}|$. For typographical
ease let $t_{0}\vcentcolon=-n$ and $t_{T+1}\vcentcolon=n$, we proceed to
analyse the integral of interest over $[-n,n]$.

$\displaystyle\int_{-n}^{n}zf(z)d\gamma(z)$
$\displaystyle=\sum_{i=0}^{T}\int_{t_{i}}^{t_{i+1}}f(z)z\frac{\exp(-\frac{1}{2}z^{2})}{\sqrt{2\pi}}dz.$

By construction $f$ is continuously differentiable on each of the above
intervals of integration. Standard integration by parts gives

$\int_{t_{i}}^{t_{i+1}}f(z)z\frac{\exp(-\frac{1}{2}z^{2})}{\sqrt{2\pi}}dz=\left[-\frac{\exp(-\frac{1}{2}z^{2})}{\sqrt{2\pi}}f(z)\right]_{t_{i}^{+}}^{t_{i+1}^{-}}+\int_{t_{i}}^{t_{i+1}}f^{\prime}(z)d\gamma(z).$

Collecting terms it follows that

$\displaystyle\int_{-n}^{n}zf(z)d\gamma(z)$
$\displaystyle=\sum_{i=0}^{T+1}\left(\left[-\frac{\exp(-\frac{1}{2}z^{2})}{\sqrt{2\pi}}f(z)\right]_{t_{i}^{+}}^{t_{i+1}^{-}}+\int_{t_{i}}^{t_{i+1}}f^{\prime}(z)d\gamma(z)\right)$
$\displaystyle=\left[-\frac{\exp(-\frac{1}{2}z^{2})}{\sqrt{2\pi}}f(z)\right]_{-n}^{n}+\sum_{i=1}^{T}\left[-\frac{\exp(-\frac{1}{2}z^{2})}{\sqrt{2\pi}}f(z)\right]_{t_{i}^{+}}^{t_{i}^{-}}+\int_{-n}^{n}f^{\prime}(z)d\gamma(z).$

Defining $f_{-}^{\prime}(z)\vcentcolon=\max\\{-f^{\prime}(z),0\\}$ and
$f_{+}^{\prime}(z)\vcentcolon=\max\\{f^{\prime}(z),0\\}$, then as
$|f_{\pm}(z)|\leq|f(z)|$ it follows that $f_{-}(z),f_{+}(z)\in
L^{1}(\mathbb{R},\mathcal{B}(\mathbb{R}),\gamma)$. As both $f_{-}$ and $f_{+}$
are nonnegative functions, $(f_{-}\textbf{1}_{[-n,n]})_{n\in\mathbb{N}}$ and
$(f_{+}\textbf{1}_{[-n,n]})_{n\in\mathbb{N}}$ are sequences of non-decreasing
functions converging to $f^{+}$ and $f^{-}$ respectively. Therefore by
monotone convergence

$\displaystyle\int_{\mathbb{R}}f^{\prime}(z)d\gamma(z)$
$\displaystyle=\int_{\mathbb{R}}f_{-}^{\prime}(z)d\gamma(z)-\int_{\mathbb{R}}f_{+}^{\prime}(z)d\gamma(z)$
$\displaystyle=\lim_{n\rightarrow\infty}\int_{-n}^{n}f_{-}^{\prime}(z)d\gamma(z)-\lim_{n\rightarrow\infty}\int_{-n}^{n}f_{+}^{\prime}(z)d\gamma(z)$
$\displaystyle=\lim_{n\rightarrow\infty}\int_{-n}^{n}f^{\prime}(z)d\gamma(z).$

As a result

$\displaystyle\int_{\mathbb{R}}zf(z)d\gamma(z)$
$\displaystyle=\lim_{n\rightarrow\infty}\int_{-n}^{n}zf(z)d\gamma(z)$
$\displaystyle=\lim_{n\rightarrow\infty}\left[-\frac{\exp(-\frac{1}{2}z^{2})}{\sqrt{2\pi}}f(z)\right]_{-n}^{n}+\sum_{i=1}^{T}\left[-\frac{\exp(-\frac{1}{2}z^{2})}{\sqrt{2\pi}}f(z)\right]_{t_{i}^{+}}^{t_{i}^{-}}+\lim_{n\rightarrow\infty}\int_{-n}^{n}f^{\prime}(z)d\gamma(z)$
$\displaystyle=\sum_{i=1}^{T}\left[-\frac{\exp(-\frac{1}{2}z^{2})}{\sqrt{2\pi}}f(z)\right]_{t_{i}^{+}}^{t_{i}^{-}}+\int_{\mathbb{R}}f^{\prime}(z)d\gamma(z)$

as claimed. ∎

##  Appendix B Section 3 proofs

###  B.1 Lemma 3.1

Lemma 3.1. Let $\phi$ be a scaled-bounded activation, see Definition 2.1, and
$\sigma_{b}^{2}>0$. Define

$W_{\phi}(q):=\frac{\mathbb{E}[\phi(\sqrt{q}Z)^{2}]}{\mathbb{E}[\phi^{\prime}(\sqrt{q}Z)^{2}]}+\sigma_{b}^{2}$

for all $q\in\mathbb{R}_{\geq 0}$. Then $W_{\phi}:\mathbb{R}_{\geq
0}\rightarrow\mathbb{R}_{\geq 0}$ and $W$ has a fixed point $q^{*}>0$.

###### Proof.

We first prove that $W_{\phi}:\mathbb{R}_{\geq 0}\rightarrow\mathbb{R}_{\geq
0}$. To bound the numerator term then for any scaled-bounded activation
function $\phi$, see Definition 2.1, it follows that there exists
$a,k\in\mathbb{R}{>0}$ such that

$0\leq\mathbb{E}[\phi(\sqrt{q}Z)^{2}]<a^{2}k^{2}<\infty.$ (26)

for all $q\in\mathbb{R}_{\geq 0}$. As
$\mathbb{E}[\phi^{\prime}(\sqrt{q}Z)^{2}]\geq 0$ it suffices to show that the
denominator is nonzero. For a given $q\in\mathbb{R}_{\geq 0}$. as
$\phi^{\prime}(\sqrt{q}z)^{2}\geq 0$ for all
$z\in\mathbb{R}\backslash\mathcal{D}$, then
$\mathbb{E}[\phi^{\prime}(\sqrt{q}Z)^{2}]=0$ implies that
$\phi^{\prime}(\sqrt{q}z)=0$ almost everywhere for $z\in\mathbb{R}$. This is a
contradiction however as there exists by construction an $a>0$ such that
$\phi(\sqrt{q}z)=kz$ for $z\in[-a/\sqrt{q},a/\sqrt{q}]$. As a result
$\mathbb{E}[\phi^{\prime}(\sqrt{q}Z)^{2}]>0$. We therefore conclude that
$W_{\phi}:\mathbb{R}_{\geq 0}\rightarrow\mathbb{R}_{\geq 0}$.

To prove that $W_{\phi}$ has a fixed point $q^{*}>0$ we need to lower and
upper bound $\mathbb{E}[\phi^{\prime}(\sqrt{q}Z)^{2}]$. A lower bound can be
derived as follows,

$\displaystyle\mathbb{E}[\phi^{\prime}(\sqrt{q}Z)^{2}]$
$\displaystyle=\int_{\mathbb{R}}\phi^{\prime}(\sqrt{q}Z)^{2}d\gamma$
$\displaystyle=2\int_{0}^{\infty}\phi^{\prime}(\sqrt{q}Z)^{2}d\gamma$
$\displaystyle=2\int_{0}^{a/\sqrt{q}}k^{2}d\gamma+2\int_{a/\sqrt{q}}^{\infty}\phi^{\prime}(\sqrt{q}Z)^{2}d\gamma$
$\displaystyle>2k\int_{0}^{a/\sqrt{q}}\phi^{\prime}(\sqrt{q}Z)^{2}d\gamma$
$\displaystyle=k^{2}\text{erf}\left(\frac{a}{\sqrt{2q}}\right).$

Here the second equality follows from the fact that integrand is an even
function and the equality on the third line by the construction of $\phi$. The
inequality on the fourth line follows from zeroing the second integral which
is positive. Furthermore $|\phi^{\prime}(z)|\leq k$ almost everywhere by
construction, therefore we conclude that

$k^{2}\text{erf}\left(\frac{a}{\sqrt{2q}}\right)\leq\mathbb{E}[\phi^{\prime}(\sqrt{q}Z)^{2}]<k^{2}<\infty$
(27)

for all $q\in\mathbb{R}_{0}$.

We now prove that $W_{\phi}$ is continuous. Due to fact that
$\mathbb{E}[\phi^{\prime}(\sqrt{q}Z)^{2}]>0$ then $W_{\phi}$ has no
singularities. It therefore suffices to prove that both
$q\mapsto\mathbb{E}[\phi(\sqrt{q}Z)^{2}]$ and
$q\mapsto\mathbb{E}[\phi^{\prime}(\sqrt{q}Z)^{2}]$ are continuous functions on
$\mathbb{R}_{\geq 0}$. In both cases we achieve this by applying Lemma A.1 in
the context of the measure space
$(\mathbb{R},\mathcal{B}(\mathbb{R}),\gamma)$, where $\gamma$ denotes the
standard one dimensional Gaussian measure and $\mathcal{B}(\mathbb{R})$ the
completion of the Borel $\sigma$ algebra on the real numbers. For
$q\mapsto\mathbb{E}[\phi(\sqrt{q}Z)^{2}]$ with $q\in\mathbb{R}_{\geq 0}$ then
condition 1 is satisfied due to (26), condition 2 follows from the continuity
of $\phi$ and condition 3 is satisfied by $g(z)\vcentcolon=a^{2}k^{2}$. For
$q\mapsto\mathbb{E}[\phi^{\prime}(\sqrt{q}Z)^{2}]$ with $q\in\mathbb{R}_{\geq
0}$ then condition 1 is satisfied due to (27), condition 2 follows from the
fact that $\phi^{\prime}$ is continuous almost everywhere in $\mathbb{R}$ and
condition 3 is satisfied for all $q\in\mathbb{R}_{\geq 0}$ by $g(z):=k^{2}$.
We conclude then that $W_{\phi}$ is continuous on $\mathbb{R}_{\geq 0}$.

To prove that $0$ is not a fixed point, observe that
$\frac{\mathbb{E}[\phi(\sqrt{q}Z)^{2}]}{\mathbb{E}[\phi^{\prime}(\sqrt{q}Z)^{2}]}\geq
0$ and therefore $W_{\phi}(q)\geq\sigma_{b}^{2}>0$ for all
$q\in\mathbb{R}_{\geq 0}$. Using (26) and (27) we may in addition derive the
following upper bound on $W_{\phi}$,

$\sigma_{b}^{2}\leq
W_{\phi}(q)<\frac{a^{2}}{\text{erf}\left(\frac{a}{\sqrt{2q}}\right)}+\sigma_{b}^{2}=\vcentcolon
U(q)$ (28)

for all $q\in\mathbb{R}_{\geq 0}$. Observe that as $W_{\phi}$ is continuous
and $W_{\phi}(0)>0$, then if $W_{\phi}$ has no fixed points it must hold that
$W_{\phi}(q)>q$ for all $q\in\mathbb{R}_{\geq 0}$. Otherwise, by the
intermediate value theorem, the function $W_{\phi}(p)-p$ must have a root and
hence $W$ must have a fixed point $q^{*}>0$. Considering the upper bound $U$
on $W_{\phi}$ from equation 28, assume that $U(q)>q$ for all
$q\in\mathbb{R}_{\geq 0}$. Then

$\displaystyle\frac{a^{2}}{\text{erf}\left(\frac{a}{\sqrt{2q}}\right)}+\sigma_{b}^{2}$
$\displaystyle>q,$ $\displaystyle
a^{2}+\sigma_{b}^{2}\text{erf}\left(\frac{a}{\sqrt{2q}}\right)$
$\displaystyle>q\text{erf}\left(\frac{a}{\sqrt{2q}}\right)$

However,
$\lim_{q\rightarrow\infty}a^{2}+\sigma_{b}^{2}\text{erf}\left(\frac{a}{\sqrt{2q}}\right)=a^{2}<\infty$
while
$\lim_{q\rightarrow\infty}q\text{erf}\left(\frac{a}{\sqrt{2q}}\right)=\infty$.
As $q\text{erf}\left(\frac{a}{\sqrt{2q}}\right)$ is continuous then there must
exist a $q\in\mathbb{R}_{\geq 0}$ such that $U(q)<q$, which is a
contradiction. We therefore conclude that $W_{\phi}(q)<U(q)<q$ for some
$q\in\mathbb{R}_{\geq 0}$ and that therefore $W_{\phi}$ must have a fixed
point $q^{*}>0$. ∎

###  B.2 Corollary 3.1.1

Corollary 3.1.1. Let $\phi$ be a scaled-bounded activation, see Definition
2.1, $\sigma_{b}^{2}>0$ and suppose

$\chi_{1}\vcentcolon=\sigma_{w}^{2}\mathbb{E}[\phi^{\prime}(\sqrt{q^{*}}Z)^{2}]=1,$

where $q^{*}>0$ is a fixed point of $W_{\phi}$, defined in (18). Then $q^{*}$
is a fixed point of the associated variance function $V_{\phi}$, defined in
(3).

###### Proof.

From Lemma 3.1 it holds that there exists $q^{*}>0$ such that
$W_{\phi}(q^{*})=q^{*}$. Inspecting (3) and (18) it follows that
$V_{\phi}(q^{*})=W_{\phi}(q^{*})=q^{*}$, therefore $q^{*}$ is a fixed point of
$V_{\phi}$. ∎

###  B.3 Lemma 3.2

Lemma 3.2. Let $\phi$ be a scaled-bounded activation, see Definition 2.1,
$\sigma_{b}^{2}>0$ and suppose

$\chi_{1}\vcentcolon=\sigma_{w}^{2}\mathbb{E}[\phi^{\prime}(\sqrt{q^{*}}Z)^{2}]=1,$

where $q^{*}>0$ is a fixed point of $W_{\phi}$. In addition, assume that all
inputs $\mathbf{x}$ are normalised so that $||\mathbf{x}||_{2}^{2}=q^{*}$.
Then the associated correlation map $R_{\phi,q*}$, defined in (5), is fixed at
each layer $l\in[L]$, satisfies $R_{\phi,q*}:[-1,1]\rightarrow[-1,1]$ and is
differentiable with

$R_{\phi,q^{*}}^{\prime}(\rho)=\sigma_{w}^{2}\mathbb{E}[\phi^{\prime}(U_{1})\phi^{\prime}(U_{2})]$

for all input correlations $\rho\in(-1,1)$.

###### Proof.

Recall that $U_{1}:=\sqrt{q^{*}}Z_{1}$ and $U_{2}:=\sqrt{q^{*}}(\rho
Z_{1}+\sqrt{1-\rho^{2}}Z_{2})$. From Corollary 3.1.1 we know that $q^{*}$ is a
fixed point of $V_{\phi}$ and therefore the variance of all inputs, given the
assumed normalisation, will remain fixed at $q^{*}$ for all layers of the
network. Since $q^{*}>0$, then $R_{\phi,q*}:[-1,1]\rightarrow[-1,1]$ by
construction as long as the correlation function (5) is finite for any
$\rho\in[-1,1]$. This follows by Cauchy-Schwarz,

$\displaystyle|\mathbb{E}[\phi(U_{1})\phi(U_{2})]|\leq\mathbb{E}[\phi(U_{1})^{2}]^{1/2}\mathbb{E}[\phi^{\prime}(U_{2})^{2}]^{1/2}<a^{2}k^{2}<\infty.$
(29)

It remains to be proved that $R_{\phi,q^{*}}$ is differentiable on $(-1,1)$
and to derive (21). To this end it suffices to show that
$H(\rho)\vcentcolon=\mathbb{E}[\phi(U_{1})\phi(U_{2})]$ is differentiable on
$(-1,1)$ and derive an expression for its derivative. We rewrite $H$ as
follows,

$\displaystyle H(\rho)$
$\displaystyle\vcentcolon=\int_{\mathbb{R}\times\mathbb{R}}\phi(u_{1})\phi(u_{2})d\gamma^{(2)}(z_{1},z_{2})$

recalling that $\gamma^{(2)}$ denotes the standard two dimensional Gaussian
measure. We proceed by applying Lemma A.2 in the context of the measure space
$(\mathbb{R}^{2},\mathcal{B}(\mathbb{R}^{2}),\gamma^{2})$, where
$\mathcal{B}(\mathbb{R}^{2})$ denotes the completion of the Borel
$\sigma$-algebra on $\mathbb{R}^{2}$, and the interval $(-1,1)$. First observe
that condition 1 of Lemma A.2 is satisfied as

$\int_{\mathbb{R}\times\mathbb{R}}|\phi(u_{1})\phi(u_{2})|d\gamma^{(2)}(z_{1},z_{2})<a^{2}k^{2}<\infty.$

For condition 2, by construction $(\phi\circ u_{2})(\rho,z_{1},z_{2})$ is non-
differentiable only on the set

$\bigcup_{i=1}^{|\mathcal{D}|}\\{(z_{1},z_{2}):\rho
z_{1}+\sqrt{1-\rho^{2}}z_{2}=\frac{d_{i}}{\sqrt{q^{*}}}\\}.$

This is the union of a finite number of one dimensional lines in
$\mathbb{R}^{2}$ and hence has measure 0. Therefore, for each $\rho\in(-1,1)$,
then $\frac{\partial\phi\circ u_{2}}{(\partial\rho)}(x,z_{1},z_{2})$ exists
almost everywhere and hence condition 2 is also satisfied. Finally, for
condition 3 note that this partial derivative can be expressed as

$\displaystyle\frac{\partial\phi\circ u_{2}}{\partial\rho}(\rho,z_{1},z_{2})$
$\displaystyle=\frac{\partial
u_{2}}{\partial\rho}(\rho,z_{1},z_{2})\phi^{\prime}(u_{2})$
$\displaystyle=\sqrt{q^{*}}\left(z_{1}-\frac{\rho}{\sqrt{1-\rho^{2}}}z_{2}\right)\phi^{\prime}(u_{2}),$

from which it follows

$\displaystyle\left|\frac{\partial\phi\circ
u_{2}}{\partial\rho}(\rho,z_{1},z_{2})\right|$
$\displaystyle\leq\sqrt{q^{*}}\left(|z_{1}|+\frac{|\rho|}{\sqrt{1-\rho^{2}}}|z_{2}|\right)|\phi^{\prime}(u_{2})|$
$\displaystyle<\sqrt{q^{*}}k\left(|z_{1}|+\frac{|\rho|}{\sqrt{1-\rho^{2}}}|z_{2}|\right).$

For any $\rho\in(-1,1)$ consider the open interval
$(-1+\delta(\rho),1-\delta(\rho))$, where $\delta(\rho):=\frac{|1-\rho|}{2}$.
It follows that

$\sup_{\rho\in(-1+\delta(\rho),1-\delta(\rho))}\left|\frac{\partial\phi\circ
u_{2}}{\partial\rho}(\rho,z_{1},z_{2})\right|\leq\sqrt{q^{*}}k\left(|z_{1}|+\frac{1-\delta(\rho)}{\sqrt{2\delta(\rho)-\delta(\rho)^{2}}}|z_{2}|\right)=\vcentcolon
g_{K}(z_{1},z_{2}).$

Letting
$\kappa_{\delta}\vcentcolon=\frac{1-\delta(\rho)}{\sqrt{2\delta(\rho)-\delta(\rho)^{2}}}$,
then by applying the Fubini-Tonelli theorem it follows that

$\displaystyle\int_{\mathbb{R}^{2}}g_{K}(z_{1},z_{2})d\gamma^{(2)}(z_{1},z_{2})$
$\displaystyle=\sqrt{q^{*}}k\left(\int_{\mathbb{R}^{2}}|z_{1}|d\gamma^{(2)}(z_{1},z_{2})+\kappa_{\delta}\int_{\mathbb{R}^{2}}|z_{2}|d\gamma^{(2)}(z_{1},z_{2})\right)$
$\displaystyle=\sqrt{q^{*}}k\sqrt{\frac{2}{\pi}}\left(1+\kappa_{\delta}\right)$
$\displaystyle<\infty.$

Hence for any $\rho\in(-1,1)$ there exists an open interval
$K=(-1+\delta(\rho),1-\delta(\rho))$ with $\rho\in K$ and an integrable
function $g_{K}(z_{1},z_{2}):\mathbb{R}^{2}\rightarrow\mathbb{R}$, such that

$\left|\frac{\partial(\phi\circ
u_{2})}{\partial\rho}(\rho,z_{1},z_{2})\right|\leq g_{K}(z_{1},z_{2}).$

We conclude then that condition 3 is also satisfied, and therefore, by Lemma
A.2, $H$ is differentiable on $(-1,1)$. The derivative of $H$ can be expressed
as follows,

$\displaystyle H^{\prime}(\rho)$
$\displaystyle=\int_{\mathbb{R}^{2}}\phi(u_{1})\frac{\partial
u_{2}}{\partial\rho}\phi^{\prime}(u_{2})d\gamma^{(2)}(z_{1},z_{2})$ (30)
$\displaystyle=\sqrt{q^{*}}\int_{\mathbb{R}}\int_{\mathbb{R}}\phi(u_{1})\left(z_{1}-\frac{\rho}{\sqrt{1-\rho^{2}}}z_{2}\right)\phi^{\prime}(u_{2})d\gamma(z_{1})d\gamma(z_{2})$
$\displaystyle=\sqrt{q^{*}}\int_{\mathbb{R}}\left(\int_{\mathbb{R}}z_{1}\phi(u_{1})\phi^{\prime}(u_{2})d\gamma(z_{1})\right)d\gamma(z_{2})-\frac{\rho\sqrt{q^{*}}}{\sqrt{1-\rho^{2}}}\int_{\mathbb{R}}\phi(u_{1})\left(\int_{\mathbb{R}}z_{2}\phi^{\prime}(u_{2})d\gamma(z_{2})\right)d\gamma(z_{1}),$

where the third equality in the above follows by applying the Fubini-Tonelli
theorem. We proceed first to derive an expression for $H^{\prime}(0)$ and then
to analyse $H^{\prime}(\rho)$ for $\rho\in(-1,0)\cup(0,1)$. In all that
follows, and as $\phi$ is odd, we let
$-d_{T}...<-d_{2}<-d_{1}<d_{1}<d_{2}..<d_{T}$ be the elements of $\mathcal{D}$
where $T\vcentcolon=|\mathcal{D}|/2$. From (30) it follows that

$\displaystyle H^{\prime}(0)$
$\displaystyle=\sqrt{q^{*}}\int_{\mathbb{R}}\left(\int_{\mathbb{R}}z_{1}\phi(\sqrt{q^{*}}z_{1})\phi^{\prime}(\sqrt{q^{*}z_{2}})d\gamma(z_{1})\right)d\gamma(z_{2})$
$\displaystyle=\sqrt{q^{*}}\int_{\mathbb{R}}\phi^{\prime}(\sqrt{q^{*}z_{2}})\left(\int_{\mathbb{R}}z_{1}\phi(\sqrt{q^{*}}z_{1})d\gamma(z_{1})\right)d\gamma(z_{2})$

By construction $\phi(\sqrt{q^{*}}z_{1})$ is differentiable at all points
other than

$\displaystyle s_{i}\vcentcolon=\frac{d_{i}}{\sqrt{q^{*}}},$ $\displaystyle
l_{i}\vcentcolon=-\frac{d_{i}}{\sqrt{q^{*}}}$

for $i\in[T]$. Applying Lemma A.3

$\displaystyle\int_{\mathbb{R}}z_{1}\phi(\sqrt{q^{*}}z_{1})d\gamma(z_{1})$
$\displaystyle=-\sum_{i}^{T}\left(\left[\frac{\exp(-\frac{1}{2}z_{1}^{2})}{\sqrt{2\pi}}\phi(\sqrt{q^{*}}z_{1})\right]_{l_{i}^{+}}^{l_{i}^{-}}+\left[\frac{\exp(-\frac{1}{2}z_{1}^{2})}{\sqrt{2\pi}}\phi(\sqrt{q^{*}}z_{1})\right]_{s_{i}^{+}}^{s_{i}^{-}}\right)$
$\displaystyle+\sqrt{q^{*}}\int_{\mathbb{R}}\phi^{\prime}(\sqrt{q^{*}}z_{1})d\gamma(z_{1})$
$\displaystyle=\sqrt{q^{*}}\int_{\mathbb{R}}\phi^{\prime}(\sqrt{q^{*}}z_{1})d\gamma(z_{1}).$

Here the second equality follows from the continuity of $\phi$. Therefore

$\displaystyle H^{\prime}(0)$
$\displaystyle=\sqrt{q^{*}}\int_{\mathbb{R}}\phi^{\prime}(\sqrt{q^{*}z_{2}})\left(\sqrt{q^{*}}\int_{\mathbb{R}}\phi^{\prime}(\sqrt{q^{*}}z_{1})d\gamma(z_{1})\right)d\gamma(z_{2})$
(31)
$\displaystyle=q^{*}\int_{\mathbb{R}^{2}}\phi^{\prime}(\sqrt{q^{*}}z_{1})\phi^{\prime}(\sqrt{q^{*}}z_{2})d\gamma^{(2)}(z_{1},z_{2})$

To derive an expression for $H^{\prime}(\rho)$ for any
$\rho\in(-1,0)\cup(0,1)$, we apply integration by parts to each of the inner
integrals in (30) using Lemma A.3. Starting with the second inner integral,
with $\rho\in(-1,0)\cup(0,1)$ and $z_{1}\in\mathbb{R}$ fixed, we define
$\psi:\mathbb{R}\rightarrow\mathbb{R}$ as the function
$z_{2}\mapsto(\phi^{\prime}\circ u_{2})(z_{1},z_{2},\rho)$. By construction
$\phi^{\prime}$ is continuously differentiable at all points other than the
following,

$\displaystyle p_{i}(z_{1})$
$\displaystyle\vcentcolon=\frac{d_{i}}{\sqrt{q^{*}(1-\rho^{2})}}-z_{1}\frac{\rho}{\sqrt{1-\rho^{2}}},$
$\displaystyle n_{i}(z_{1})$
$\displaystyle\vcentcolon=-\frac{d_{i}}{\sqrt{q^{*}(1-\rho^{2})}}-z_{1}\frac{\rho}{\sqrt{1-\rho^{2}}}$

for all $i\in[T]$. Applying Lemma A.3 it follows that

$\displaystyle\int_{\mathbb{R}}z_{2}\psi(z_{2})d\gamma(z_{2})$
$\displaystyle=-\sum_{i}^{T}\left(\left[\frac{\exp(-\frac{1}{2}z_{2}^{2})}{\sqrt{2\pi}}\psi(z_{2})\right]_{n_{i}^{+}}^{n_{i}^{-}}+\left[\frac{\exp(-\frac{1}{2}z_{2}^{2})}{\sqrt{2\pi}}\psi(z_{2})\right]_{p_{i}^{+}}^{p_{i}^{-}}\right)+\int_{\mathbb{R}}\psi^{\prime}(z_{2})d\gamma(z_{2})$
$\displaystyle=\vcentcolon-\kappa_{1}(z_{1})+\int_{\mathbb{R}}\psi^{\prime}(z_{2})d\gamma(z_{2}).$

Therefore the second integral on the final line of (30) can be expressed as

$\displaystyle\int_{\mathbb{R}}\phi(u_{1})\left(\int_{\mathbb{R}}z_{2}\phi^{\prime}(u_{2})d\gamma(z_{2})\right)d\gamma(z_{1})$
$\displaystyle=-\int_{\mathbb{R}}\phi(u_{1})\kappa_{1}(z_{1})d\gamma(z_{1})+\int_{\mathbb{R}}\phi(u_{1})\int_{\mathbb{R}}\psi^{\prime}(z_{2})d\gamma(z_{2})d\gamma(z_{1})$
$\displaystyle=-\int_{\mathbb{R}}\phi(u_{1})\kappa_{1}(z_{1})d\gamma(z_{1})+\sqrt{q^{*}(1-\rho^{2})}\int_{\mathbb{R}^{2}}\phi(u_{1})\phi^{\prime\prime}(u_{2})d\gamma^{(2)}(z_{1},z_{2}).$

Analysing $\kappa_{1}(z_{1})$, then by construction
$\psi(n_{i})=\phi^{\prime}(-d_{i})$ and $\psi(p_{i})=\phi^{\prime}(d_{i})$.
Furthermore, as $\phi$ is odd and continuous then $\phi^{\prime}$ is even and
as a result $\phi^{\prime}(-d_{i}^{-})=\phi^{\prime}(d_{i}^{+})$ and
$\phi^{\prime}(-d_{i}^{+})=\phi^{\prime}(d_{i}^{-})$. Therefore

$\displaystyle\kappa_{1}(z_{1})$
$\displaystyle=\sum_{i}^{T}\left(\left[\frac{\exp(-\frac{1}{2}z^{2})}{\sqrt{2\pi}}\psi(z)\right]_{n_{i}^{+}}^{n_{i}^{-}}+\left[\frac{\exp(-\frac{1}{2}z^{2})}{\sqrt{2\pi}}\psi(z)\right]_{p_{i}^{+}}^{p_{i}^{-}}\right)$
$\displaystyle=\sum_{i}^{T}\left(\frac{\exp(-\frac{1}{2}n_{i}^{2})}{\sqrt{2\pi}}\left(\psi(n_{i}^{-})-\psi(n_{i}^{+})\right)+\frac{\exp(-\frac{1}{2}p_{i}^{2})}{\sqrt{2\pi}}\left(\psi(p_{i}^{-})-\psi(p_{i}^{+})\right)\right)$
$\displaystyle=\sum_{i}^{T}\left(\frac{\exp(-\frac{1}{2}n_{i}^{2})}{\sqrt{2\pi}}\left(\phi^{\prime}(-d_{i}^{-})-\phi^{\prime}(-d_{i}^{+})\right)+\frac{\exp(-\frac{1}{2}p_{i}^{2})}{\sqrt{2\pi}}\left(\phi^{\prime}(d_{i}^{-})-\phi^{\prime}(d_{i}^{+})\right)\right)$
$\displaystyle=\sum_{i}^{T}\left(\frac{\exp(-\frac{1}{2}n_{i}^{2})}{\sqrt{2\pi}}\left(\phi^{\prime}(d_{i}^{+})-\phi^{\prime}(d_{i}^{-})\right)+\frac{\exp(-\frac{1}{2}p_{i}^{2})}{\sqrt{2\pi}}\left(\phi^{\prime}(d_{i}^{-})-\phi^{\prime}(d_{i}^{+})\right)\right)$
$\displaystyle=\frac{1}{\sqrt{2\pi}}\sum_{i}^{T}\left(\phi^{\prime}(d_{i}^{+})-\phi^{\prime}(d_{i}^{-})\right)\left(\exp(-\frac{1}{2}n_{i}^{2})-\exp(-\frac{1}{2}p_{i}^{2})\right).$

Expanding the integral involving $\kappa_{1}$ then

$\displaystyle\int_{\mathbb{R}}\phi(u_{1})\kappa_{1}(z_{1})d\gamma(z_{1})=\sum_{i}^{T}\left(\phi^{\prime}(d_{i}^{+})-\phi^{\prime}(d_{i}^{-})\right)\int_{\mathbb{R}}\phi(u_{1})\frac{\exp(-\frac{1}{2}n_{i}^{2})-\exp(-\frac{1}{2}p_{i}^{2})}{\sqrt{2\pi}}d\gamma(z_{1}).$

Observe that as

$\displaystyle
p_{i}(z_{1})^{2}=\frac{d_{i}^{2}}{q^{*}(1-\rho^{2})}-z_{1}\frac{2d_{i}\rho}{\sqrt{q^{*}}(1-\rho^{2})}+z_{1}^{2}\frac{\rho^{2}}{1-\rho^{2}}$
$\displaystyle
n_{i}(z_{1})^{2}=\frac{d_{i}^{2}}{q^{*}(1-\rho^{2})}+z_{1}\frac{2d_{i}\rho}{\sqrt{q^{*}}(1-\rho^{2})}+z_{1}^{2}\frac{\rho^{2}}{1-\rho^{2}}$

then clearly $n_{i}(z_{1})^{2}=p_{i}(-z_{1})^{2}$. Let
$\beta_{i}(z_{1})\vcentcolon=e^{-\frac{1}{2}n_{i}(z_{1})^{2}}-e^{-\frac{1}{2}p_{i}(z_{1})^{2}}$
for all $i\in[T]$, then these $\beta_{i}$ are odd functions as

$\displaystyle\beta_{i}(-z_{1})$
$\displaystyle=e^{-\frac{1}{2}n_{i}(-z_{1})^{2}}-e^{-\frac{1}{2}p_{i}(-z_{1})^{2}}$
$\displaystyle=e^{-\frac{1}{2}p_{i}(z_{1})^{2}}-e^{-\frac{1}{2}n_{i}(z_{1})^{2}}$
$\displaystyle=-\beta_{i}(z_{1}).$

As the product of two odd functions is odd then

$\displaystyle\int_{\mathbb{R}}\phi(u_{1})\kappa_{1}(z_{1})d\gamma(z_{1})=\frac{1}{\sqrt{2\pi}}\sum_{i}^{T}\left(\phi^{\prime}(d_{i}^{+})-\phi^{\prime}(d_{i}^{-})\right)\int_{\mathbb{R}}\phi(u_{1})\beta_{i}(z_{1})d\gamma(z_{1})=0$

and so

$\int_{\mathbb{R}}\phi(u_{1})\left(\int_{\mathbb{R}}z_{2}\phi^{\prime}(u_{2})d\gamma(z_{2})\right)d\gamma(z_{1})=\sqrt{q^{*}(1-\rho^{2})}\int_{\mathbb{R}^{2}}\phi(u_{1})\phi^{\prime\prime}(u_{2})d\gamma^{(2)}(z_{1},z_{2}).$
(32)

Now we turn our attention to the first inner integral on the final line of
(30). With $\rho\in(-1,0)\cup(0,1)$ and $z_{1}\in\mathbb{R}$ fixed, and
defining for typographical ease $\varphi:\mathbb{R}\rightarrow\mathbb{R}$ as
the function $z_{1}\mapsto(\phi^{\prime}\circ u_{2})(z_{1},z_{2},\rho)$, then
$\phi(u_{1})\phi^{\prime}(u_{2})$ is continuously differentiable at all points
other than the following,

$\displaystyle
o_{i}(z_{2})\vcentcolon=\frac{d_{i}}{\rho\sqrt{q^{*}}}-z_{2}\frac{\sqrt{1-\rho^{2}}}{\rho}$
$\displaystyle
e_{i}(z_{2})\vcentcolon=-\frac{d_{i}}{\rho\sqrt{q^{*}}}-z_{2}\frac{\sqrt{1-\rho^{2}}}{\rho}$
$\displaystyle s_{i}\vcentcolon=\frac{d_{i}}{\sqrt{q^{*}}}$ $\displaystyle
l_{i}\vcentcolon=-\frac{d_{i}}{\sqrt{q^{*}}}$

for all $i\in[T]$. Applying Lemma A.3,

$\displaystyle\int_{\mathbb{R}}z_{1}\phi(\sqrt{q^{*}}z_{1})\varphi(z_{1})d\gamma(z_{1})$
$\displaystyle=-\sum_{i}^{T}\left(\left[\frac{\exp(-\frac{1}{2}z_{1}^{2})}{\sqrt{2\pi}}\phi(\sqrt{q^{*}}z_{1})\varphi(z_{1})\right]_{e_{i}^{+}}^{e_{i}^{-}}+\left[\frac{\exp(-\frac{1}{2}z_{1}^{2})}{\sqrt{2\pi}}\phi(\sqrt{q^{*}}z_{1})\varphi(z_{1})\right]_{o_{i}^{+}}^{o_{i}^{-}}\right)$
$\displaystyle-\sum_{i}^{T}\left(\left[\frac{\exp(-\frac{1}{2}z_{1}^{2})}{\sqrt{2\pi}}\phi(\sqrt{q^{*}}z_{1})\varphi(z_{1})\right]_{l_{i}^{+}}^{l_{i}^{-}}+\left[\frac{\exp(-\frac{1}{2}z_{1}^{2})}{\sqrt{2\pi}}\phi(\sqrt{q^{*}}z_{1})\varphi(z_{1})\right]_{s_{i}^{+}}^{s_{i}^{-}}\right)$
$\displaystyle+\int_{\mathbb{R}}\phi(\sqrt{q^{*}}z_{1})\varphi^{\prime}(z_{1})+\sqrt{q^{*}}\phi^{\prime}(\sqrt{q^{*}}z_{1})\varphi(z_{1})d\gamma(z_{1}).$

As both $\phi$ and $\varphi$ are continuous at $l_{i}$ and $s_{i}$ for all
$i\in[T]$, then the left and right limits of
$\phi(\sqrt{q^{*}}z_{1})\varphi(z_{1})$ at these points are equal. Therefore

$\displaystyle\int_{\mathbb{R}}z_{1}\phi(\sqrt{q^{*}}z_{1})\varphi(z_{1})d\gamma(z_{1})$
$\displaystyle=-\sum_{i}^{T}\left(\left[\frac{\exp(-\frac{1}{2}z_{1}^{2})}{\sqrt{2\pi}}\phi(\sqrt{q^{*}}z_{1})\varphi(z_{1})\right]_{e_{i}^{+}}^{e_{i}^{-}}+\left[\frac{\exp(-\frac{1}{2}z_{1}^{2})}{\sqrt{2\pi}}\phi(\sqrt{q^{*}}z_{1})\varphi(z_{1})\right]_{o_{i}^{+}}^{o_{i}^{-}}\right)$
$\displaystyle+\int_{\mathbb{R}}\phi(\sqrt{q^{*}}z_{1})\varphi^{\prime}(z_{1})+\sqrt{q^{*}}\phi^{\prime}(\sqrt{q^{*}}z_{1})\varphi(z_{1})d\gamma(z_{1})$
$\displaystyle=\vcentcolon-\kappa_{2}(z_{2})+\int_{\mathbb{R}}\phi(\sqrt{q^{*}}z_{1})\varphi^{\prime}(z_{1})+\sqrt{q^{*}}\phi^{\prime}(\sqrt{q^{*}}z_{1})\varphi(z_{1})d\gamma(z_{1}).$

Analysing $\kappa_{2}$, then by construction
$\varphi(e_{i})=\phi^{\prime}(-d_{i})$ and
$\varphi(o_{i})=\phi^{\prime}(d_{i})$. Similar to before it follows that

$\displaystyle\kappa_{2}(z_{2})$
$\displaystyle=\sum_{i=1}^{T}\left(\left[\frac{\exp(-\frac{1}{2}z_{1}^{2})}{\sqrt{2\pi}}\phi(\sqrt{q^{*}}z_{1})\varphi(z_{1})\right]_{e_{i}^{+}}^{e_{i}^{-}}+\left[\frac{\exp(-\frac{1}{2}z_{1}^{2})}{\sqrt{2\pi}}\phi(\sqrt{q^{*}}z_{1})\varphi(z_{1})\right]_{o_{i}^{+}}^{o_{i}^{-}}\right)$
$\displaystyle=\sum_{i=1}^{T}\left(\frac{\exp(-\frac{1}{2}e_{i}^{2})}{\sqrt{2\pi}}\phi(\sqrt{q^{*}}e_{i})\left(\varphi(-d_{i}^{-})-\varphi(-d_{i}^{+})\right)+\frac{\exp(-\frac{1}{2}o_{i}^{2})}{\sqrt{2\pi}}\phi(\sqrt{q^{*}}o_{i})\left(\varphi(d_{i}^{-})-\varphi(d_{i}^{+})\right)\right)$
$\displaystyle=\sum_{i=1}^{T}\left(\frac{\exp(-\frac{1}{2}e_{i}^{2})}{\sqrt{2\pi}}\phi(\sqrt{q^{*}}e_{i})\left(\varphi(d_{i}^{+})-\varphi(d_{i}^{-})\right)+\frac{\exp(-\frac{1}{2}o_{i}^{2})}{\sqrt{2\pi}}\phi(\sqrt{q^{*}}o_{i})\left(\varphi(d_{i}^{-})-\varphi(d_{i}^{+})\right)\right)$
$\displaystyle=\frac{1}{\sqrt{2\pi}}\sum_{i=1}^{T}\left(\left(\varphi(d_{i}^{+})-\varphi(d_{i}^{-})\right)\left(\exp(-\frac{1}{2}e_{i}^{2})\phi(\sqrt{q^{*}}e_{i})-\exp(-\frac{1}{2}o_{i}^{2})\phi(\sqrt{q^{*}}o_{i})\right)\right)$

Note that $o_{i}(-z_{2})=-e_{i}(z_{2})$, therefore
$\phi(\sqrt{q^{*}}o_{i}(-z_{2}))=\phi(-\sqrt{q^{*}}e_{i}(z_{2}))=-\phi(\sqrt{q^{*}}e_{i}(z_{2}))$
as $\phi$ is odd. Likewise $e_{i}(-z_{2})=-o_{i}(z_{2})$ and so
$\phi(\sqrt{q^{*}}e_{i}(-z_{2}))=\phi(-\sqrt{q^{*}}o_{i}(z_{2}))=-\phi(\sqrt{q^{*}}o_{i}(z_{2}))$.
Furthermore as

$\displaystyle o_{i}(z_{2})^{2}$
$\displaystyle=\frac{d_{i}^{2}}{\rho^{2}q^{*}}-z_{2}\frac{\sqrt{1-\rho^{2}}}{\rho}+z_{2}^{2}\frac{1-\rho^{2}}{\rho^{2}},$
$\displaystyle e_{i}(z_{2})^{2}$
$\displaystyle=\frac{d_{i}^{2}}{\rho^{2}q^{*}}+z_{2}\frac{\sqrt{1-\rho^{2}}}{\rho}+z_{2}^{2}\frac{1-\rho^{2}}{\rho^{2}}$

then $o_{i}(-z_{2})^{2}=e_{i}(z_{2})^{2}$. Analogous to $\beta_{i}$, we now
define
$\Gamma_{i}(z_{2})\vcentcolon=\exp(-\frac{1}{2}e_{i}^{2})\phi(\sqrt{q^{*}}e_{i})-\exp(-\frac{1}{2}o_{i}^{2})\phi(\sqrt{q^{*}}o_{i})$.
The fact that $\Gamma_{i}(z_{2})$ is odd follows from

$\displaystyle\Gamma_{i}(-z_{2})$
$\displaystyle=\exp(-\frac{1}{2}e_{i}(-z_{2})^{2})\phi(\sqrt{q^{*}}e_{i}(-z_{2})))-\exp(-\frac{1}{2}o_{i}(-z_{2})^{2})\phi(\sqrt{q^{*}}o_{i}(-z_{2}))$
$\displaystyle=-\exp(-\frac{1}{2}o_{i}(z_{2})^{2})\phi(\sqrt{q^{*}}o_{i}(z_{2})))+\exp(-\frac{1}{2}e_{i}(z_{2})^{2})\phi(\sqrt{q^{*}}e_{i}(z_{2}))$
$\displaystyle=-\Gamma_{i}(z_{2}).$

Analysing the first integral on the last line of 30), then as $\Gamma_{i}$ is
odd for each $i\in[T]$

$\displaystyle\int_{\mathbb{R}}\left(\int_{\mathbb{R}}z_{1}\phi(u_{1})\phi^{\prime}(u_{2})d\gamma(z_{1})\right)d\gamma(z_{2})$
(33)
$\displaystyle=-\int_{\mathbb{R}}\kappa_{2}(z_{2})d\gamma(z_{2})+\int_{\mathbb{R}}\int_{\mathbb{R}}\phi(\sqrt{q^{*}}z_{1})\varphi^{\prime}(z_{1})+\sqrt{q^{*}}\phi^{\prime}(\sqrt{q^{*}}z_{1})\varphi(z_{1})d\gamma(z_{1})d\gamma(z_{2})$
$\displaystyle=-\frac{1}{\sqrt{2\pi}}\sum_{i=1}^{T}\left(\varphi(d_{i}^{+})-\varphi(d_{i}^{-})\right)\int_{\mathbb{R}}\Gamma_{i}(z_{2})d\gamma(z_{2})+\int_{\mathbb{R}^{2}}\phi(\sqrt{q^{*}}z_{1})\varphi^{\prime}(z_{1})d\gamma^{(2)}(z_{1},z_{2})$
$\displaystyle+\sqrt{q^{*}}\int_{\mathbb{R}^{2}}\phi^{\prime}(\sqrt{q^{*}}z_{1})\varphi(z_{1})d\gamma^{(2)}(z_{1},z_{2})$
$\displaystyle=\sqrt{q^{*}}\rho\int_{\mathbb{R}^{2}}\phi(u_{1})\phi^{\prime\prime}(u_{2})d\gamma^{(2)}(z_{1},z_{2})+\sqrt{q^{*}}\int_{\mathbb{R}^{2}}\phi^{\prime}(u_{1})\phi^{\prime}(u_{2})d\gamma^{(2)}(z_{1},z_{2}).$

Substituting (33) and (32) into (30) it follows that

$\displaystyle H^{\prime}(\rho)$
$\displaystyle=\sqrt{q^{*}}\left(\int_{\mathbb{R}}\left(\int_{\mathbb{R}}z_{1}\phi(u_{1})\phi^{\prime}(u_{2})d\gamma(z_{1})\right)d\gamma(z_{2})-\frac{\rho}{\sqrt{1-\rho^{2}}}\int_{\mathbb{R}}\phi(u_{1})\left(\int_{\mathbb{R}}z_{2}\phi^{\prime}(u_{2})d\gamma(z_{2})\right)d\gamma(z_{1})\right)$
$\displaystyle=q^{*}\left(\rho\int_{\mathbb{R}^{2}}\phi(u_{1})\phi^{\prime\prime}(u_{2})d\gamma^{(2)}(z_{1},z_{2})+\int_{\mathbb{R}^{2}}\phi^{\prime}(u_{1})\phi^{\prime}(u_{2})d\gamma(z_{1},z_{2})-\rho\int_{\mathbb{R}^{2}}\phi(u_{1})\phi^{\prime\prime}(u_{2})d\gamma^{(2)}(z_{1},z_{2})\right)$
$\displaystyle=q^{*}\int_{\mathbb{R}^{2}}\phi^{\prime}(u_{1})\phi^{\prime}(u_{2})d\gamma^{(2)}(z_{1},z_{2}).$

It therefore follows for all $\rho\in(-1,1)$ that

$\displaystyle
R_{\phi,q^{*}}^{\prime}(\rho)=\sigma_{w}^{2}\int_{\mathbb{R}^{2}}\phi^{\prime}(u_{1})\phi^{\prime}(u_{2})d\gamma^{(2)}(z_{1},z_{2})$
$\displaystyle=\sigma_{w}^{2}\mathbb{E}[\phi^{\prime}(U_{1})\phi^{\prime}(U_{2})]$

as claimed. ∎

###  B.4 Lemma 3.3

Lemma 3.3. Under the same conditions and assumptions as in Lemma 3.2, it holds
that

$\max_{\rho\in[0,1]}|R_{\phi,q^{*}}(\rho)-\rho|=\frac{\sigma_{b}^{2}}{q^{*}}.$

###### Proof.

Observe that

$\displaystyle R_{\phi,q^{*}}(0)$
$\displaystyle=\frac{\mathbb{E}[\phi(\sqrt{q^{*}}Z_{1})\phi(\sqrt{q^{*}}Z_{2})]}{q^{*}\mathbb{E}[\phi^{{}^{\prime}}(\sqrt{q^{*}}Z)^{2}]}+\frac{\sigma_{b}^{2}}{q^{*}}$
$\displaystyle=\frac{\mathbb{E}[\phi(\sqrt{q^{*}}Z_{1})]^{2}}{q^{*}\mathbb{E}[\phi^{{}^{\prime}}(\sqrt{q^{*}}Z)^{2}]}+\frac{\sigma_{b}^{2}}{q^{*}}.$

As $\phi$ is odd then $\mathbb{E}[\phi(\sqrt{q^{*}}Z_{1})]=0$ and therefore

$R_{\phi,q^{*}}(0)=\frac{\sigma_{b}^{2}}{q^{*}}>0.$

Lemma 3.3 follows then as long as
$|R_{\phi,q^{*}}(\rho)-\rho|<R_{\phi,q^{*}}(0)$ for all $\rho\in[0,1]$. As

$\displaystyle
R_{\phi,q^{*}}(1)=\frac{\sigma_{w}^{2}}{q^{*}}\mathbb{E}[\phi(\sqrt{q^{*}}Z_{1})^{2}]+\frac{\sigma_{b}^{2}}{q^{*}}=\frac{V(q^{*})}{q^{*}}=1$

then $|R_{\phi,q^{*}}(1)-1|=0<R_{\phi,q^{*}}(0)$. All that remains to show is
that the inequality holds for $\rho\in(0,1)$. We proceed using an approach
similar to that used to prove Proposition 3 in [10]. Using Lemma 3.2 then for
any $\rho\in(0,1)$ we have

$\displaystyle R_{\phi,q^{*}}^{\prime}(\rho)$
$\displaystyle=\frac{\mathbb{E}[\phi^{\prime}(U_{1})\phi^{\prime}(U_{2})]}{\mathbb{E}[\phi^{{}^{\prime}}(\sqrt{q^{*}}Z)^{2}]}$
$\displaystyle\leq\frac{\mathbb{E}[\phi^{\prime}(\sqrt{q^{*}}Z_{1})^{2}]^{\frac{1}{2}}\mathbb{E}[\phi^{\prime}(U_{2})^{2}]^{\frac{1}{2}}}{\mathbb{E}[\phi^{{}^{\prime}}(\sqrt{q^{*}}Z)^{2}]}$
$\displaystyle=\frac{\mathbb{E}[\phi^{\prime}(U_{2})^{2}]^{\frac{1}{2}}}{\mathbb{E}[\phi^{{}^{\prime}}(\sqrt{q^{*}}Z)^{2}]^{\frac{1}{2}}}$
$\displaystyle=1.$

The inequality on the second line of the above follows from Cauchy-Schwarz and
the equalities on the third and fourth lines are due to the fact that
$Z_{1},Z,U_{2}\sim\mathcal{N}(0,1)$ are all identically distributed. Note that
equality holds iff either $\rho=0$ or there exists an $\alpha\in\mathbb{R}$
such that $\phi^{\prime}(U_{1})=\alpha\phi^{\prime}(U_{2})$. Since $Z_{1}$ and
$Z_{2}$ are i.i.d. this can only occur if $\phi^{\prime}$ is a constant, which
in turn would imply that $\phi$ must be linear. However, by construction
linear functions are clearly not scaled-bounded activations and therefore for
any $\rho\in(0,1)$ it holds that

$R_{\phi,q^{*}}^{\prime}(\rho)<1.$

For $\rho\in(0,1)$ then integrating both sides of the above inequality and
applying the fundamental theorem of calculus we have

$\int_{0}^{\rho}R_{\phi,q^{*}}^{\prime}(t)dt<\int_{0}^{\rho}1dt\implies
R_{\phi,q^{*}}(\rho)-\rho<R_{\phi,q^{*}}(0)$

and

$\int_{\rho}^{1}R_{\phi,q^{*}}^{\prime}(t)dt<\int_{\rho}^{1}1dt\implies\rho-
R_{\phi,q^{*}}(\rho)<0.$

As $R_{\phi,q^{*}}(0)>0$ then we conclude for $\rho\in(0,1)$ that
$|R_{\phi,q^{*}}(\rho)-\rho|<R_{\phi,q^{*}}(0)$. Therefore

$|R_{\phi,q^{*}}(\rho)-\rho|\leq
R_{\phi,q^{*}}(0)=\frac{\sigma_{b}^{2}}{q^{*}}$

for all $\rho\in[0,1]$ as claimed. ∎

###  B.5 Lemma 3.4

Lemma 3.4. Under the same conditions and assumptions as in Lemma 3.2, and
defining $y:=\frac{\sigma_{b}^{2}}{a^{2}}$, then

$\Lambda(y)<\frac{a}{\sqrt{q^{*}}}<\left(\frac{8}{\pi}\right)^{1/6}y^{-1/3},$

where $\Lambda(y)$ is defined as

$\left(W_{0}\left(\frac{2}{\pi}\left(\sqrt{\frac{8}{\pi}}\exp\left(-\frac{\left(\sqrt{W_{0}\left(\frac{2}{\pi}y^{-2}\right)}\right)^{2}}{2}\right)\left(\frac{1}{\left(\sqrt{W_{0}\left(\frac{2}{\pi}y^{-2}\right)}\right)}+\left(\sqrt{W_{0}\left(\frac{2}{\pi}y^{-2}\right)}\right)\right)\right)^{-2}\right)\right)^{1/2}$

and $W_{0}$ denotes the principal branch of the Lambert W function.

###### Proof.

To derive the upper bound on $a/\sqrt{q^{*}}$ we study lower bounds for
$V_{\phi}(q)$. To this end we first lower bound
$\mathbb{E}[\phi(\sqrt{q}Z)^{2}]$. The fact that $\phi$ is odd implies
$\phi^{2}$ is even, therefore

$\displaystyle\mathbb{E}[\phi(\sqrt{q}Z)^{2}]$
$\displaystyle=2\left(\int_{0}^{a/\sqrt{q}}k^{2}qz^{2}d\gamma(z)+\int_{a/\sqrt{q}}^{\infty}\phi(\sqrt{q}z)^{2}d\gamma(z)\right)$
$\displaystyle\geq 2\left(k^{2}q\int_{0}^{a/\sqrt{q}}z^{2}d\gamma(z)\right)$
$\displaystyle=qk^{2}\text{erf}\left(\frac{a}{\sqrt{2q}}\right)-\sqrt{\frac{2}{\pi}}ak^{2}\sqrt{q}\exp\left(-\frac{a^{2}}{2q}\right).$

To now upper bound $\mathbb{E}[\phi^{{}^{\prime}}(\sqrt{q*}Z)^{2}]$ we use the
fact that $|\phi^{\prime}(z)|\leq k$, which implies

$\displaystyle\mathbb{E}[\phi^{{}^{\prime}}(\sqrt{q^{*}}Z)^{2}]$
$\displaystyle\leq 2k^{2}\int_{0}^{\infty}d\gamma(z)$ $\displaystyle=k^{2}.$

As $q^{*}=V_{\phi}(q^{*})$, it therefore follows that

$q^{*}\geq
q^{*}\text{erf}\left(\frac{a}{\sqrt{2q^{*}}}\right)-\sqrt{\frac{2}{\pi}}a\sqrt{q^{*}}\exp\left(-\frac{a^{2}}{2q^{*}}\right)+\sigma_{b}^{2}.$

Rearranging and dividing by $a^{2}$ gives

$\frac{q^{*}}{a^{2}}\text{erfc}\left(\frac{a}{\sqrt{2q^{*}}}\right)\geq\frac{\sigma_{b}^{2}}{a^{2}}-\frac{\sqrt{q^{*}}}{a}\sqrt{\frac{2}{\pi}}\exp\left(-\frac{a^{2}}{2q^{*}}\right).$

For typographical ease we now substitute $x=\frac{a}{\sqrt{q^{*}}}$ and
multiply by $x^{2}$,

$\text{erfc}\left(\frac{x}{\sqrt{2}}\right)\geq
x^{2}\frac{\sigma_{b}^{2}}{a^{2}}-x\sqrt{\frac{2}{\pi}}\exp\left(-\frac{x^{2}}{2}\right).$

It is known that
$\text{erfc}\left(\frac{x}{\sqrt{2}}\right)\leq\sqrt{\frac{2}{\pi}}\frac{1}{x}\exp\left(-\frac{x^{2}}{2}\right)\leq\sqrt{\frac{2}{\pi}}\frac{1}{x}$
(see e.g., [31]). Additionally it holds that
$\exp\left(-\frac{x^{2}}{2}\right)<\frac{1}{x^{2}}$, this follows from the
fact that $x>2\ln(x)$, which can be proved using elementary calculus. Using
these results we formulate the following inequality,

$x^{2}\frac{\sigma_{b}^{2}}{a^{2}}-x\sqrt{\frac{2}{\pi}}\frac{1}{x^{2}}<\sqrt{\frac{2}{\pi}}\frac{1}{x},$

which simplifies to

$\frac{a}{\sqrt{q^{*}}}=\vcentcolon
x<\left(\frac{8}{\pi}\right)^{1/6}\frac{a^{2/3}}{\sigma_{b}^{2/3}}=\left(\frac{8}{\pi}\right)^{1/6}y^{-1/3}.$

as claimed.

The derivation of the lower bound is more involved, we proceed by deriving an
upper bound on $V_{\phi}(q)$ which will allow us to upper bound $q^{*}$.
First, as $\phi$ is upper bounded by $k^{2}a^{2}$ then

$\displaystyle\mathbb{E}[\phi(\sqrt{q}Z)^{2}]$ $\displaystyle\leq
2\left(\int_{0}^{a/\sqrt{q}}k^{2}qz^{2}d\gamma(z)+k^{2}a^{2}\int_{a/\sqrt{q}}^{\infty}d\gamma(z)\right)$
$\displaystyle=k^{2}q\text{erf}\left(\frac{a}{\sqrt{2q}}\right)-k^{2}\sqrt{\frac{2}{\pi}}a\sqrt{q}\exp\left(-\frac{a^{2}}{2q}\right)+k^{2}a^{2}\text{erfc}\left(\frac{a}{\sqrt{2q}}\right).$

A simple lower bound for $\mathbb{E}[\phi^{{}^{\prime}}(\sqrt{q^{*}}Z)^{2}]$
is as follows,

$\displaystyle\mathbb{E}[\phi^{{}^{\prime}}(\sqrt{q^{*}}Z)^{2}]$
$\displaystyle\geq 2k^{2}\int_{0}^{a/\sqrt{q^{*}}}\gamma(z)$
$\displaystyle=k^{2}\text{erf}\left(\frac{a}{\sqrt{2q^{*}}}\right).$

As a result of these inequalities we may formulate the following inequality,

$q^{*}\leq\frac{k^{2}q^{*}\text{erf}\left(\frac{a}{\sqrt{2q^{*}}}\right)-k^{2}\sqrt{\frac{2}{\pi}}a\sqrt{q^{*}}\exp\left(-\frac{a^{2}}{q^{*}}\right)+k^{2}a^{2}\text{erfc}\left(\frac{a}{\sqrt{2q^{*}}}\right)}{k^{2}\text{erf}\left(\frac{a}{\sqrt{2q^{*}}}\right)}+\sigma_{b}^{2}.$

Rearranging this expression gives

$q^{*}\text{erf}\left(\frac{a}{\sqrt{2q^{*}}}\right)\leq
q^{*}\text{erf}\left(\frac{a}{\sqrt{2q^{*}}}\right)-\sqrt{\frac{2}{\pi}}a\sqrt{q^{*}}\exp\left(-\frac{a^{2}}{2q^{*}}\right)+a^{2}\text{erfc}\left(\frac{a}{\sqrt{2q^{*}}}\right)+\sigma_{b}^{2}\text{erf}\left(\frac{a}{\sqrt{2q^{*}}}\right).$

Further simplification leads to

$\displaystyle\sqrt{\frac{2}{\pi}}a\sqrt{q^{*}}\exp\left(-\frac{a^{2}}{2q^{*}}\right)$
$\displaystyle\leq(a^{2}-\sigma_{b}^{2})\text{erfc}\left(\frac{a}{\sqrt{2q^{*}}}\right)+\sigma_{b}^{2}$
$\displaystyle<a^{2}\text{erfc}\left(\frac{a}{\sqrt{2q^{*}}}\right)+\sigma_{b}^{2}.$

Dividing by $a^{2}$ and making the substitution $x=\frac{a}{\sqrt{q^{*}}}$ we
arrive at the key inequality

$\sqrt{\frac{2}{\pi}}\frac{1}{x}\exp\left(-\frac{x^{2}}{2}\right)<\text{erfc}\left(\frac{x}{\sqrt{2}}\right)+\frac{\sigma_{b}^{2}}{a^{2}}.$
(34)

Recalling that our objective is to lower bound the quantity
$\frac{a}{\sqrt{q^{*}}}$, then by construction any value of $x$ which
satisfies the above inequality is a viable candidate. The challenging aspect
now is to find a non-trivial candidate, in particular one that scales
appropriately with $a$ and $\sigma_{b}^{2}$. For typographical ease we now
define
$g(x):=\sqrt{\frac{2}{\pi}}\frac{1}{x}\exp\left(-\frac{x^{2}}{2}\right)$ and
$h(x):=\text{erfc}\left(\frac{x}{\sqrt{2}}\right)+\frac{\sigma_{b}^{2}}{a^{2}}$.
The derivatives of each of these function are as follows,

$\displaystyle g^{\prime}(x)$
$\displaystyle=-\sqrt{\frac{2}{\pi}}\exp\left(-\frac{x^{2}}{2}\right)\left(\frac{1}{x^{2}}+1\right)$
$\displaystyle h^{\prime}(x)$
$\displaystyle=-\sqrt{\frac{2}{\pi}}\exp\left(-\frac{x^{2}}{2}\right),$

note that $g^{\prime}(x)>h^{\prime}(x)$ for all $x\in\mathbb{R}$. We now
identify a candidate lower bound by considering the linear equation $u(x)$
defined by being tangent to $g(x)$ at
$x=\delta:=\sqrt{W_{0}\left(\frac{2}{\pi}y^{-2}\right)}$. Here $W_{0}$ refers
to the principle branch of the Lambert W function (see e.g., [32]). The line
$u(x)$ therefore has the form

$u(x)=-\sqrt{\frac{2}{\pi}}\exp\left(-\frac{\delta^{2}}{2}\right)\left(\frac{1}{\delta^{2}}+1\right)x+\beta.$

Observe that as $u(\delta)=y=g(\delta)$ by construction, then

$\displaystyle\beta$
$\displaystyle=y+\sqrt{\frac{2}{\pi}}\exp\left(-\frac{\delta^{2}}{2}\right)\left(\frac{1}{\delta}+\delta\right)$
$\displaystyle=\sqrt{\frac{2}{\pi}}\frac{1}{\delta}\exp\left(-\frac{\delta^{2}}{2}\right)+\sqrt{\frac{2}{\pi}}\exp\left(-\frac{\delta^{2}}{2}\right)\left(\frac{1}{\delta}+\delta\right)$
$\displaystyle=\sqrt{\frac{2}{\pi}}\exp\left(-\frac{\delta^{2}}{2}\right)\left(\frac{2}{\delta}+\delta\right)$

We now identify a candidate lower bound $\gamma$ as the solution of
$g(\gamma)=\beta$, the solution of which can be expressed using the Lambert W
function, $\gamma=\sqrt{W_{0}\left(\frac{2}{\pi}\beta^{-2}\right)}$. Consider
the following compositions of functions,
$\Phi_{g}(y):=g(\gamma(\beta(\delta(y))))$ and
$\Psi_{h}(y):=h(\gamma(\beta(\delta(y))))$. To ensure that
$\gamma\leq\frac{a}{\sqrt{q^{*}}}$, i.e., that $\gamma$ is indeed a valid
lower bound, then it must hold that $\Phi_{g}(y)\geq\Phi_{h}(y)$. We refer the
reader to Figure 9 for a graphical description of this approach to identifying
a candidate lower bound. Our task then is to inspect the range of
$y\in\mathbb{R}_{\geq 0}$ values for which this inequality holds true. First
we note that $\Phi_{f}(0)=\Phi_{h}(0)=0$. Second we observe $h,g,\gamma,\beta$
and $\delta$ are differentiable functions on $\mathbb{R}_{>0}$. As

$\displaystyle\frac{d\delta}{dy}$
$\displaystyle=-\frac{\sqrt{W_{0}\left(\frac{2}{\pi}y^{-2}\right)}}{y\left(W_{0}\left(\frac{2}{\pi}y^{-2}\right)+1\right)}<0,$
$\displaystyle\frac{d\beta}{d\delta}$
$\displaystyle=-\frac{\sqrt{\frac{2}{\pi}}\exp\left(-\frac{\delta^{2}}{2}\right)\left(\delta^{4}+\delta^{2}+2\right)}{\delta^{2}}<0,$
$\displaystyle\frac{d\gamma}{d\beta}$
$\displaystyle=-\frac{\sqrt{W_{0}\left(\frac{2}{\pi}\beta^{-2}\right)}}{\beta\left(W_{0}\left(\frac{2}{\pi}\beta^{-2}\right)+1\right)}<0$

then
$\frac{d\gamma}{dy}=\frac{d\gamma}{d\beta}\frac{d\beta}{d\delta}\frac{d\delta}{dy}<0$.
Applying the chain rule

$\displaystyle\Phi_{g}^{\prime}(y)$
$\displaystyle=\frac{dg}{d\gamma}\frac{d\gamma}{dy}$
$\displaystyle=-\sqrt{\frac{2}{\pi}}\exp\left(-\frac{\gamma^{2}}{2}\right)\left(\frac{1}{\gamma^{2}}+1\right)\frac{d\gamma}{dy}$
$\displaystyle=\sqrt{\frac{2}{\pi}}\exp\left(-\frac{\gamma^{2}}{2}\right)\left(\frac{1}{\gamma^{2}}+1\right)\left|\frac{d\gamma}{dy}\right|$

and

$\displaystyle\Phi_{h}^{\prime}(y)$
$\displaystyle=\frac{dh}{d\gamma}\frac{d\gamma}{dy}$
$\displaystyle=\sqrt{\frac{2}{\pi}}\exp\left(-\frac{\gamma^{2}}{2}\right)\left|\frac{d\gamma}{dy}\right|.$

The fact that

$\displaystyle\sqrt{\frac{2}{\pi}}\exp\left(-\frac{\gamma^{2}}{2}\right)\left(\frac{1}{\gamma^{2}}+2\right)\left|\frac{d\gamma}{dy}\right|>\sqrt{\frac{2}{\pi}}\exp\left(-\frac{\gamma^{2}}{2}\right)\left|\frac{d\gamma}{dy}\right|,$

which can be further simplified to

$\left(\frac{1}{\gamma^{2}}+1\right)>1,$

holds for all $\gamma\in\mathbb{R}$ implies that
$\Phi_{g}^{\prime}(y)>\Phi_{h}^{\prime}(y)$ for all $y\in\mathbb{R}_{>0}$.
Therefore

$\displaystyle\int_{0}^{y}\Phi_{g}^{\prime}(t)-\Phi_{h}^{\prime}(t)dt$
$\displaystyle>0$

for all $y\in\mathbb{R}_{>0}$. Now applying the fundamental theorem of
calculus

$\displaystyle\int_{0}^{y}\Phi_{g}^{\prime}(t)-\Phi_{h}^{\prime}(t)dt$
$\displaystyle=\Phi_{g}(y)-\Phi_{h}(y)+\Phi_{h}(0)-\Phi_{f}(0)$
$\displaystyle=\Phi_{g}(y)-\Phi_{h}(y)$

we conclude that $\Phi_{g}(y)>\Phi_{h}(y)$ for all $y\in\mathbb{R}_{>0}$. As a
result, $\gamma$ is a valid lower bound for $\frac{a^{2}}{\sqrt{q^{*}}}$ as
long as $\frac{\sigma_{b}^{2}}{a^{2}}>0$. Finally, to recover the statement of
the theorem, we define the composite function $\Lambda:\mathbb{R}_{\geq
0}\rightarrow\mathbb{R}_{\geq 0}$ as
$\Lambda(y)\vcentcolon=\gamma(\beta(\delta(y)))$. ∎

Figure 9: illustration of proof approach to finding a lower bound $\gamma$ for
$\frac{a}{\sqrt{q^{*}}}$ for a given value of $y=\sigma_{b}^{2}/a^{2}$ (in
this example $y=0.25$). The red line indicates the value on the x-axis defined
to be $\gamma$. This is computed by first identifying where the line $u(x)$,
which lies tangent to the point where $g$ hits the asymptotic limit of $h(x)$
which in turn corresponds to the value of $y$, intercepts the y-axis, denoted
as $\beta$. We define $\gamma$ then to be the point at which the horizontal
line at $\beta$ (the upper dashed brown line) intercepts $g(x)$. To ensure
that $\beta$ intercepts $g(x)$ above the critical point of the inequality
$g(x)<h(x)$, it is necessary to check that $h(\gamma)<g(\gamma)$ is true.

###  B.6 Convergence of preactivation correlations with depth for a range of
activation functions

(a)  Figure 10: evolution of numerically computed preactivation correlations
from two initial correlations, 0.01 and 0.6, across a range of activation
functions and fixed point values.