---
abstract: 'We study randomly initialized residual networks using mean field theory and the theory of difference equations. Classical feedforward neural networks, such as those with tanh activations, exhibit exponential behavior on the average when propagating inputs forward or gradients backward. The exponential forward dynamics causes rapid collapsing of the input space geometry, while the exponential backward dynamics causes drastic vanishing or exploding gradients. We show, in contrast, that by adding skip connections, the network will, depending on the nonlinearity, adopt subexponential forward and backward dynamics, and in many cases in fact polynomial. The exponents of these polynomials are obtained through analytic methods and proved and verified empirically to be correct. In terms of the “edge of chaos” hypothesis, these subexponential and polynomial laws allow residual networks to “hover over the boundary between stability and chaos,” thus preserving the geometry of the input space and the gradient information flow. In our experiments, for each activation function we study here, we initialize residual networks with different hyperparameters and train them on MNIST. Remarkably, our [*initialization time*]{} theory can accurately predict [*test time*]{} performance of these networks, by tracking either the expected amount of gradient explosion or the expected squared distance between the images of two input vectors. Importantly, we show, theoretically as well as empirically, that common initializations such as the Xavier or the He schemes are not optimal for residual networks, because [*the optimal initialization variances depend on the depth*]{}. Finally, we have made mathematical contributions by deriving several new identities for the kernels of powers of ReLU functions by relating them to the zeroth Bessel function of the second kind.'
author:
- |
    Greg Yang[^1]\
    [Microsoft Research AI]{}\
    `gregyang@microsoft.com`\
    Samuel S. Schoenholz\
    [Google Brain]{}\
    `schsam@google.com`\
bibliography:
- 'neural\_dynamics.bib'
title: 'Mean Field Residual Networks: On the Edge of Chaos'
---

Introduction
============

Previous works [@poole_exponential_2016; @daniely_toward_2016; @schoenholz_deep_2017] have shown that randomly initialized neural networks exhibit a spectrum of behavior with depth, from stable to chaotic, which depends on the variance of the initializations: the cosine distance of two input vectors converges exponentially fast with depth to a fixed point in \[0, 1\]; if this fixed point is 1, then the behavior is stable; if this fixed point is 0, then the behavior is chaotic. It has been argued in many prior works [@bertschinger_real-time_2004; @poole_exponential_2016] that effective computation can only be supported by a dynamical behavior that is on [**the edge of chaos**]{}. Too much stability prevents the neural network from telling apart two different inputs. While some chaotic behavior can increase the expressivity of a network, too much chaos makes the neural network think two similar inputs are very different. At the same time, the same initialization variances also control how far gradient information can be propagated through the network; the networks with chaotic forward dynamics will tend to suffer from exploding gradients, while networks with stable forward dynamics will tend to suffer from vanishing gradients.

These works have focused on vanilla (fully connected) feedforward networks. Here we consider residual networks [@he_deep_2016; @he_identity_2016] (with fully-connected layers and without batchnorm), which are a family of recently proposed neural network architectures that has achieved state-of-the-art performance on image recognition tasks, beating all other approaches by a large margin. The main innovation of this family of architectures is the addition of a passthrough (identity) connection from the previous layer to the next, such that the usual nonlinearity computes the “residual” between the next-layer activation and the previous-layer activation.

In this work, we seek to characterize randomly initialized residual networks. One of our main results is that random residual networks for many nonlinearities such as $\tanh$ [**live on the edge of chaos**]{}, in that the cosine distance of two input vectors will converge to a fixed point at a polynomial rate, rather than an exponential rate, as with vanilla tanh networks. Thus a typical residual network will slowly cross the stable-chaotic boundary with depth, hovering around this boundary for many layers. In addition, for most of the nonlinearities considered here, the mean field estimate of the gradient grows subexponentially with depth. In fact, for $\alpha$-ReLU, the $\alpha$th-power of ReLU, for $\alpha < 1$, the gradient grows only polynomially. These theoretical results provide some theoretical justification for why residual networks work so well in practice. In our experiments, we are also able to predict surprisingly well the relative performances of [*trained*]{} residual networks based only on their initialization hyperparameters, in a variety of settings. In particular, we find that the quality of initialization for tanh resnets is determined by [*trainability*]{} (how much gradient explosion on average) while that for ($\alpha$-)ReLU resnets is determined by expressivity (how far can two different input vectors be pulled apart) (see \[sec:experiments\]). To the best of our knowledge, this is the first time that a quantity other than gradient explosion/vanishing has been found to control the quality of initialization. We establish theoretically and empirically that the best initialization variances for residual networks depend on the depth of the network (contrary to the feedforward case [@schoenholz_deep_2017]), so that **common initialization schemes like Xavier [@glorot_understanding_2010] or He [@he_delving_2015] cannot be optimal**. In fact, even the rationale of He initialization is incorrect for ReLU residual networks because it tries to control gradient dynamics rather than expressivity. However we want to emphasize that we study a simplified model of residual networks in this work, with no batchnorm or convolutional layers, so that these results are not necessarily indicative of the MSRA residual network used in practice [@he_deep_2016].

In the body of this paper, we give account of general intuition and/or proof strategy when appropriate for our theoretical results, but we relegate all formal statements and proofs to the appendix.

Background {#sec:background}
==========

Consider a vanilla feedforward neural network of $L$ layers, with each layer $l$ having $\p N l$ neurons; here layer 0 is the input layer. For the ease of presentation we assume all hidden layer widths are the same $\p N l = N$ for all $l > 0$. Let $\p x 0 = (\p x 0_1, \ldots, \p x 0_{\p N 0})$ be the input vector to the network, and let $\p x l$ for $l > 0$ be the activation of layer $l$. Then a neural network is given by the equations $$\begin{aligned}
        \p x l _i &= \phi(\p h l _i), & \p h l _i &= \sum_{j=1}^N \p w l _{ij} \p x {l-1}_j + \p b l_i
    \end{aligned}$$ where

$\p h l$ is the pre-activation at layer $l$,

$\p w l$ is the weight matrix,

$\p b l$ is the bias vector, and

$\phi$ is a nonlinearity, for example $\tanh$ or ReLU, which is applied coordinatewise to its input.

To lighten up notation, we suppress the explicit layer numbers $l$ and write $$\begin{aligned}
        x_i &= \phi(h_i), & h_i &= \sum_j w_{ij} {\underline}x_j + b_i
    \end{aligned}$$ where $\bullet$ implicitly denotes $\p \bullet l$, and ${\underline}\bullet$ denotes $\p \bullet {l-1}$ (and analogously, ${\overline}\bullet$ denotes $\p \bullet {l+1}$).

A series of papers [@poole_exponential_2016; @raghu_expressive_2016; @schoenholz_deep_2017] investigated the “average behavior” of random neural networks sampled via $\p w l_{ij} \sim {\mathcal{N}}(0, \sigma_w^2/N), \p b l_i \sim {\mathcal{N}}(0, \sigma_b^2)$, for fixed parameters $\sigma_w$ and $\sigma_b$, independent of $l$. Consider the expectation of $\f 1 N \sum_{i=1}^N x^2_i$, the normalized squared length of $x$, over the sampling of $w$ and $b$. @poole_exponential_2016 showed that this quantity converges to a fixed point exponentially fast for sigmoid nonlinearities. Now suppose we propagate two different vectors $\p x 0$ and $(\p x 0)'$ through the network. @poole_exponential_2016 also showed that the expectation of the normalized dot product $\f 1 N \sum_{i=1}^N x_i x'_i$ converges exponentially fast to a fixed point. The ratio between the normalized squared length and the normalized dot product is the cosine distance between $x$ and $x'$. Thus these two exponential convergence results show that the cosine distance converges exponentially fast to a fixed point as well. Intuitively, this means that a vanilla feedforward network “forgets” the geometry of the input space “very quickly,” after only a few layers.

In addition, @schoenholz_deep_2017, under certain independence assumptions, showed that the expected normalized squared norm of the gradient also vanishes or explodes in an exponential fashion with depth, with the ”half-life” controlled by $\sigma_w$ and $\sigma_b$. They verified that this theoretical ”half-life” correlates in practice with the maximal number of layers that are admissible to good performance.

At the same time, @daniely_toward_2016 published work of similar nature, but phrased in the language of reproducing kernel Hilbert spaces, and provided high probability estimates that are meaningful for the case when the width $N$ is finite and the depth is logarithmic in $N$. However, they essentially fixed the variance parameters $\sigma_\bullet$, and furthermore, their framework (for example the notion of a “skeleton”) does not immediately generalize to the residual network case.

In this work, we show that residual networks have very different dynamics from vanilla feedforward networks. In most cases, the cosine distance convergence rate and the gradient growth rate are subexponential in a residual network, and in most cases, these rates may be polynomial.

Preliminaries
=============

Residual networks were first introduced by [@he_deep_2016] and later refined by [@he_identity_2016], and they are now commonplace among deployed neural systems. The key innovation there is the addition of a shortcut connection from the previous layer to the next. We define the following idealized architectures for ease of analysis. Note that we only consider fully-connected affine layers instead of convolutional layers. A [**reduced residual network (RRN)**]{} has the recurrence $$\begin{aligned}
        x_i &= \phi(h_i) + {\underline}x, & h_i &= \sum_j w_{ij} {\underline}x_j + b_i.
    \end{aligned}$$ A [**(full) residual network (FRN)**]{} in addition has an affine connection given by weights $v$ and biases $a$ from the nonlinearity $\phi(h)$ to the next layer: $$\begin{aligned}
    x_i &= \sum_j v_{ij} \phi( h_j) + {\underline}x_i + a_i, &  h_i &= \sum_j w_{ij} {\underline}x_j + b_i
    \end{aligned}$$

We are interested in the “average behavior” of these network when the weights and biases, $\p w l_{ij}, \p b l_i, \p v l_{ij}$, and $\p a l_i$ are sampled i.i.d. from Gaussian distributions resp. with standard deviations $\sigma_w, \sigma_b, \sigma_v,$ and $\sigma_a$, independent from $l$. Here we take the variance of $\p w l_{ij}$ to be $\sigma_w^2/N$ so that the variance of each $h_i$ is $\sigma_w^2$, assuming each ${\underline}x_j$ is fixed (similarity for $\p v l_{ij}$). Such an initialization scheme is standard in practice.

We make several key “physical assumptions” to make theoretical computations tractable:

\[ass:symAct\] (a) We assume $\la (\p h l _i)^2 \ra = \la (\p h l _j)^2 \ra$ and $\la (\p x 0 _i)^2 \ra = \la (\p x 0 _j)^2 \ra$ for any $i, j, l$. (b) We also assume that the gradient $\pd E/\pd \p x l _i$ with respect to the loss function $E$ satisfies $\la (\pd E/\pd \p x l _i)^2 \ra = \la (\pd E/ \pd \p x l _j)^2 \ra$ for any $i, j, l$.

One can see that \[ass:symAct\](a) is satisfied if the input $\p x 0 \in \{\pm 1\}^N$ and \[ass:symAct\](b) is satisfied if \[ass:gradInd\] below is true and the gradient at the last layer $\pd E/\pd x L \in \{\pm 1 \}^N$. But in general it is justified both empirically and theoretically as an approximation, because $(\p h l _i)^2 - (\p h l _j)^2$ stays about constant with $l$, but $(\p h l _i)^2$ and $(\p h l _j)^2$ grow rather quickly at the same pace with $l$ (as will be seen later in calculations), so that their additive difference becomes negligible; similarly for $(\p x l _i)^2$ and $(\pd E/\pd \p h l _i)^2$.

\[ass:gradInd\] (a) We assume the we use a different set of weights for backpropagation than those used to compute the network outputs, but sampled i.i.d. from the same distributions. (b) For any loss function $E$, we assume that the gradient at layer $l$, $\pd E/\pd \p x l _i$, is independent from all activations $\p h {l} _j$ and $\p x {l-1} _j$ from the previous layer.

\[ass:gradInd\](a) was first made in [@schoenholz_deep_2017] for computing the mean field theory of gradients for feedforward tanh networks. This is similar to the practice of feedback alignment [@lillicrap_random_2016]. Even though we are the first to explicitly formulate \[ass:gradInd\](b), in fact it was already applied implicitly in the gradient calculations of [@schoenholz_deep_2017]. Note that a priori \[ass:gradInd\](b) is not true, as $\pd E/\pd \p x l _i$ depends on $\dot \phi(\p h {l+1} _k)$ for every $k$, which depend on $\p h {l} _j$ for each $j$, and which depends on $\p x {l-1} _k$ for every $k$. Nevertheless, in practice both subassumptions hold very well.

Now we define the central quantities studied in this paper. Inevitably, our paper involves a large amount of notation that may be confusing for the first-time reader. We have included a glossary of symbols (\[tab:glossary\]) to ameliorate notation confusion.

\[defn:length\] Fix an input $\p x 0$. Define the [**length quantities**]{} $\p {\mathbf{q}}l := \la (\p h l_1)^2 \ra$ and $\p {\mathbf{p}}l := \la (\p x l_1)^2 \ra$ for $l > 0$ and $\p {\mathbf{p}}0 = \|\p x 0\|^2/N$. Here the expectations $\la \bullet \ra$ are taken over all random initialization of weights and biases for all layers $l$, as $N \to \infty$ (large width limit).

Note that in our definition, the index $1$ does not matter by \[ass:symAct\].

\[defn:corr\] Fix two inputs $\p x 0$ and $\p x 0{}'$. We write $\bullet'$ to denote a quantity $\bullet$ with respect to the input ${\p x 0}'$. Then define [**the correlation quantities**]{} $\p {\boldsymbol{\oldgamma}}l:= \la \p h l_1 \p h l_1{}' \ra$ and $\p {\boldsymbol{\oldlambda}}l:= \la \p x l_1 \p x l_1{}'\ra$ for $l > 0$ and $\p {\boldsymbol{\oldgamma}}0 = \p x 0 \cdot \p x 0 {}' / N$, where the expectations $\la \bullet \ra$ are taken over all random initialization of weights and biases for all layers $l$, as $N \to \infty$ (large width limit). Again, here the index $1$ does not matter by \[ass:symAct\]. By [**metric expressivity**]{}, we mean $\p {\mathbf{s}}l := \f 1 {2N} \la \|\p x l - \p x l {}'\|^2\ra = \f 1 {2N} (\la \|\p x l\|^2\ra + \la \|\p x l {}'\|^2 \ra - 2 \la \p x l \cdot \p x l {}'\ra) = \f 1 2 (\p {\mathbf{p}}l + \p {\mathbf{p}}l {}') - \p {\boldsymbol{\oldgamma}}l$. Additionally, define [**the cosine distance quantities**]{} $\p {\mathbf{e}}l := \p {\boldsymbol{\oldgamma}}l / \sqrt{\p {\mathbf{p}}l \p {\mathbf{p}}l {}'}$ and $\p {\mathbf{c}}l := \p {\boldsymbol{\oldlambda}}l / \sqrt{\p {\mathbf{q}}l \p {\mathbf{q}}l{}'}$, and we will also call $\p {\mathbf{e}}l$ [**angular expressivity**]{}.

In this paper, for the ease of presentation, we assume $\p {\mathbf{p}}0 = \p {\mathbf{p}}0 {}'$. Then, as we will see, $\p {\mathbf{p}}l = \p {\mathbf{p}}l{}', \p {\mathbf{q}}l = \p {\mathbf{q}}l{}'$ for all $l$, and as a result, $\p {\mathbf{e}}l = \p {\boldsymbol{\oldgamma}}l / \p {\mathbf{p}}l$ and $\p {\mathbf{s}}l = \p {\mathbf{p}}l - \p {\boldsymbol{\oldgamma}}l = (1 - \p {\mathbf{e}}l) \p {\mathbf{p}}l$.

\[defn:grad\] Fix an input $\p x 0$ and a gradient vector $(\pd E/ \pd{\p x L_i})_i$ of some loss function $E$ with respect to the last layer $\p x L$. Then define [**the gradient quantities**]{} $\p {{\boldsymbol{\oldchi}}}l:= \la (\pd E/\pd \p x l _1)^2 \ra, \p {\boldsymbol{\oldchi}}l _\bullet := \la (\pd E/\pd \p \bullet l _1)^2 \ra$ for $\bullet = a, b$, and $\p {\boldsymbol{\oldchi}}l _\bullet := \la (\pd E/\pd \p \bullet l _{11})^2 \ra$ for $\bullet = w, v$. Here the expectations are taken with \[ass:gradInd\] in mind, over both random initialization of forward and backward weights and biases, as $N \to \infty$ (large width limit). Again, the index $1$ or $11$ does not matter by \[ass:symAct\].

#### Asymptotic notations.

The expressions $f = O(g) \iff g = \Omega(f)$ have their typical meanings, and $f = \Theta(g)$ iff $f = O(g), g = O(f)$. We take $f(x) = \tilde O(g(x)) \iff g(x) = \tilde \Omega(f(x))$ to mean $f(x) = O(g\log^k x)$ for some $k \in \Z$ (this is slightly different from the standard usage of $\tilde O$), and $f = \tilde\Theta(g) \iff f = \tilde O(g) \And g = \tilde O(f).$ We introduce a new notation: $f = {{\check\Theta}}(g)$ if $f(x) = O(g(x) \cdot x^{\epsilon})$ and $f(x) = \Omega(g(x) \cdot x^{-{\epsilon}})$, as $x \to \infty$, for any ${\epsilon}> 0$. All asymptotic notations are sign-less, i.e. can indicate either positive or negative quantities, unless stated otherwise.

Overview
========

The primary reason we may say anything about the average behavior of any of the above quantities is the central limit theorem: every time the activations of the previous layer pass through an affine layer whose weights are sampled i.i.d., the output is a sum of a large number of random variables, and thus follows approximately Gaussian distributions. The mean and variance of these distributions can be computed by keeping track of the mean and variances of the activations in the previous layer.

In what follows, we use this technique to derive recurrence equations governing ${\mathbf{p}}, {\mathbf{q}}, {\boldsymbol{\oldgamma}}, {\boldsymbol{\oldlambda}}, {{\boldsymbol{\oldchi}}}$ for different architectures and different activation functions. We use these equations to investigate the dynamics of ${\mathbf{e}}$ and ${\mathbf{s}}$, the key quantities in the forward pass, and the dynamics of ${{\boldsymbol{\oldchi}}}$, the key quantity in the backward pass.

The cosine distance ${\mathbf{e}}$ in some sense measures the angular geometry of two vectors. If ${\mathbf{e}}= 1$, then the vectors are parallel; if ${\mathbf{e}}= 0$, then they are orthogonal. Just as in [@poole_exponential_2016] and [@schoenholz_deep_2017], we will show that in all of the architectures and activations we consider in this paper, $\p {\mathbf{e}}l$ converges to a fixed point ${\mathbf{e}}^*$ as $l \to \infty$ [^2]. Thus, on the average, as vectors propagate through network, the geometry of the original input space, for example, linear separability, is “forgotten” by residual networks as well as by vanilla networks. But we will prove and verify experimentally that, while @poole_exponential_2016 and [@schoenholz_deep_2017] showed that the convergence rate to ${\mathbf{e}}^*$ is exponential in a vanilla network, the convergence rate is rather only polynomial in residual networks, for tanh and $\alpha$-ReLU (\[defn:alphaReLU\]) nonlinearities; see \[thm:dalethRecReduced\], \[thm:eDynamicsFullResTanh\], \[thm:ReLUSquaredConvergence\], and \[thm:alphaReLUeConvergence\]. This slow convergence preserves geometric information in the input space, and allows a typical residual network to “hover over the edge of chaos”: Even when the cosine distance $\p{\mathbf{e}}l$ converges to 0, corresponding to “chaos”, (resp. 1, corresponding to “stability”), for the number of layers usually seen in practice, $\p {\mathbf{e}}l$ will reside well away from 0 (resp. 1).

Similarly, the quantity ${\mathbf{s}}$ measures the metric geometry of two vectors. The evolution of $\p {\mathbf{s}}l$ with $l$ tells us the ability of the average network to separate two input points in terms of Euclidean distance. Again, for tanh and $\alpha$-ReLU ($\alpha < 1$) nonlinearities, ${\mathbf{s}}$ varies only polynomially with $l$.

On the other hand, $\p {{\boldsymbol{\oldchi}}}l$ measures the size of gradient at layer $l$, and through it we track the dynamics of gradient backpropagation, be it explosion or vanishing. In contrast to vanilla tanh networks, which can experience both of these two phenomenon depending on the initialization variances, typical residual networks cannot have vanishing gradient, in the sense of vanishing $\p {{\boldsymbol{\oldchi}}}l$ as $l \to 1$; see \[thm:dalethRecReduced\] and \[thm:dalethRecFull\]. Furthermore, while vanilla tanh networks exhibit exponentially vanishing or exploding gradients, all of the activation/architecture pairings considered here, except the full residual network with ReLU, have subexponential gradient dynamics. While tanh residual networks (reduced or full) has $\p {{\boldsymbol{\oldchi}}}0 \approx \exp(\Theta(\sqrt l)) \p {{\boldsymbol{\oldchi}}}l$ (\[thm:dalethExpSqrtTanhFullRes\]), $\alpha$-ReLU residual networks for $\alpha < 1$ have $\p {{\boldsymbol{\oldchi}}}0 \approx {\mathsf{poly}}(l) \p {{\boldsymbol{\oldchi}}}l$ (\[thm:dalethDynamicsAlphaReLU\]). Instead of $\pd E/\pd x_i$, we may also consider the size of gradients of actual trainable parameters. For tanh and $\alpha$-ReLU with $\alpha < 1$, they are still subexponential and polynomial (\[thm:alphaReLUAllGradients\]). On the other hand, while $\p {{\boldsymbol{\oldchi}}}0 = \exp(\Theta(l))\p {{\boldsymbol{\oldchi}}}l$ for a ReLU resnet, its weight gradients have size independent of layer, within $O(1)$ (\[thm:alphaReLUAllGradients\])! This is the only instance in this paper of gradient norm being completely preserved across layers.

The above overviews the theoretical portion of this paper. Through experiments, we discover that we can very accurately predict whether one random initialization leads to better performance than another on the test set, after training, by leveraging this theory we build. Residual networks of different nonlinearities have different [*controlling quantities*]{}: for resnets with tanh, the optimal initialization is obtained by controlling the gradient explosion $\p {{\boldsymbol{\oldchi}}}0 / \p {{\boldsymbol{\oldchi}}}L$; whereas for ReLU and $\alpha$-ReLU, the optimal initialization is obtained by maximizing ${\mathbf{s}}$ without running into numerical issues (with floating point computation). See \[sec:experiments\] for details.

Over the course of our investigation of $\alpha$-ReLU, we derived several new identities involving the associated kernel functions, first defined in [@cho_kernel_2009], which relate them to the zeroth Bessel functions (\[lemma:JalphaBessel,lemma:LAlphaRec,lemma:JalphaRec,lemma:JalphaGrad\]).

Theoretical Results
===================

In what follows in the main text, we assume $\sigma_\bullet > 0$ for all $\bullet = w, v, b, a$; in the appendix, the formal statement of each main theorem will contain results for other cases. We are interested in the two major categories of nonlinearities used today: tanh-like and rectified units. We make the following formal definitions as a foundation for further consideration.

We say a function $\phi$ is [**tanh-like**]{} if $\phi$ is antisymmetric ($\phi(-x) = -\phi(x)$), $|\phi(x)| \le 1$ for all $x$, $\phi(x) \ge 0, \forall x \ge 0$, and $\phi(x)$ monotonically increases to 1 as $x \to \infty$.

\[defn:alphaReLU\] Define the $\alpha$-ReLU $\psi_\alpha(x) = x^\alpha$ if $x > 0$ and 0 otherwise. [^3]

  Antisymmetric/RRN                                                                            Any/FRN                                                                                 
  -------------------------------------------------------------------------------------------- --------------------------------------------------------------------------------------- --
                                                                                                                                                                                       
  Theorems \[thm:p\_q\_linear\], \[thm:lambda\_gamma\_recurrence\], \[thm:dalethRecReduced\]   Theorems \[thm:fullResPQRec\], \[thm:full\_res\_l\_g\_recurr\], \[thm:dalethRecFull\]   

  : Main Recurrences[]{data-label="tab:recurrences"}

By applying the central limit theorem as described in the last section, we derive a set of recurrences for different activation/architecture pairs, shown in \[tab:recurrences\] (see appendix for proofs). They leverage certain integral transforms [^4] as in the following

\[defn:integralTransform\] Define the transforms ${\mathrm{V}}$ and ${\mathrm{W}}$ by ${\mathrm{V}}\phi(q) := \EV[\phi(z)^2: z \sim {\mathcal{N}}(0, q)]$ and ${\mathrm{W}}\phi(\rho, \nu) := \EV[\phi(z)\phi(z'): (z, z') \sim {\mathcal{N}}(0, \begin{pmatrix}\rho & \nu \\ \nu & \rho \end{pmatrix})]$.

These recurrences are able to track the corresponding quantities in practice very well. For example, \[fig:theory\_tracks\_pratice\] compares theory vs experiments for the tanh/FRN pair. The agreement is very good for tanh/RRN (not shown, but similar to the case of tanh/FRN with $\sigma_v = 1$ and $\sigma_a = 0$) and $\alpha$-ReLU/FRN as well (see \[fig:alphaReLUTheoryVsEmpirics\]).

As mentioned in previous sections, we seek to characterize the long term/high depth behavior of all of the quantities defined in \[sec:background\]. To do so, we solve for the asymptotics of the recurrences in \[tab:recurrences\], where $\phi$ is instantiated with tanh or $\alpha$-ReLU. Our main dynamics results are summarized in \[tab:dynamics\].

                                        Tanh/RRN                                                Tanh/FRN                                             ReLU/FRN                                             $\alpha$-ReLU/FRN, $\alpha<1$                                                           
  ------------------------------------- ------------------------------------------------------- ---------------------------------------------------- ---------------------------------------------------- --------------------------------------------------------------------------------------- --
  $\p {\mathbf{p}}l$                    $\Theta(l)$, \[thm:p\_q\_linear\]                       $\Theta(l)$, \[thm:pIsLinearTanh\]                   $\exp(\Theta(l))$, \[thm:pDynamicAlphaReLU\]         $\Theta(l^{1/(1-\alpha)})$, \[thm:pDynamicAlphaReLU\]                                   
  $\p {\mathbf{s}}l$                    $\Theta(l)$, \[thm:edynamics\]                          $\Theta(l)$, \[thm:eDynamicsFullResTanh\]            $\exp(\Theta(l))$, \[thm:ReLUSquaredConvergence\]    $\Theta(l^{1/(1-\alpha)})$, \[thm:alphaReLUeConvergence\]                               
  $\p {\mathbf{e}}l - {\mathbf{e}}^*$   ${{\check\Theta}}(l^{\f2 \pi -1})$, \[thm:edynamics\]   ${\mathsf{poly}}(l)$, \[thm:eDynamicsFullResTanh\]   $\Theta(l^{-2})$, \[thm:ReLUSquaredConvergence\]     ${\mathsf{poly}}(l)$, \[thm:alphaReLUeConvergence\]                                     
  $\p {{\boldsymbol{\oldchi}}}l$        $\exp(\Theta(\sqrt l))$, \[thm:dalethExpSqrtTanh\]      $\exp(\Theta(\sqrt l))$, \[thm:dalethRecFull\]       $\exp(\Theta(l))$, \[thm:dalethDynamicsAlphaReLU\]   $\Theta(l^{\f{\alpha^2}{(1-\alpha)(2 \alpha - 1)}})$, \[thm:dalethDynamicsAlphaReLU\]   

  : Summary of Main Dynamics Results. Note that while $\p{{\boldsymbol{\oldchi}}}l$ is exponential for ReLU/FRN, the gradients with respect to weight parameters have norms (${\boldsymbol{\oldchi}}_w$ and ${\boldsymbol{\oldchi}}_v$) constant in $l$ (\[thm:alphaReLUAllGradients\]). Also, the $\p {{\boldsymbol{\oldchi}}}l$ entry for $\alpha$-ReLU is for $\alpha \in (3/4, 1)$ only[]{data-label="tab:dynamics"}

Tanh
----

![Our equations predict the relevant quantities very well in practice. These plots make the comparison between prediction and measurements for the full resnet with tanh activation, with $\sigma_v^2 = 1.5$, $\sigma_a^2 = .5$, $\sigma_w^2 = 1.69$, $\sigma_b^2 = .49$. Left-to-right: [**(a)**]{} $\p {\mathbf{p}}l$ and $\p {\boldsymbol{\oldgamma}}l$ against layer $l$ for 200 layers. [**(b)**]{} $\p {\mathbf{e}}l = \p {\boldsymbol{\oldgamma}}l /\p {\mathbf{p}}l$ against $l$ for 200 layers. Both (a) and (b) trace out curves for different initial conditions. [**(c)**]{} Different gradient quantities against $l$ for 50 layers. From left to right the layer number $l$ decreases, following the direction of backpropagation. Notice that the gradient increases in norm as $l \to 1$. All three figures exhibit smooth curves, which are theoretical estimates, and irregular curves with shades around them, which indicate empirical means and standard deviations (both of which taken in regular scale, not log scale). (a) and (b) are made with 20 runs of resnets of width 1000. (c) is made with 25 runs of resnets of width 250.[]{data-label="fig:theory_tracks_pratice"}](graphics/tanh_full_res_theory_vs_exp-chi.pdf){height=".17\textheight"}

#### Forward dynamics.

When $\phi = \tanh$, $\p {\mathbf{p}}l$ and $\p {\mathbf{q}}l$ increase as $\Theta(l)$ in either RRN or FRN (\[thm:p\_q\_linear\]), as one might expect by observing that ${\mathrm{V}}\tanh( {\mathbf{q}}) \to 1$ as ${\mathbf{q}}\to \infty$ so that, for example in the RRN case, the recurrence ${\mathbf{p}}= {\mathrm{V}}\tanh( {\mathbf{q}}) + {\underline}{\mathbf{p}}$ becomes ${\mathbf{p}}= 1 + {\underline}{\mathbf{p}}$. This is confirmed graphically by the black lines of the leftmost chart of \[fig:theory\_tracks\_pratice\]. We carefully verify that this intuition is correct in its proof in the appendix, and find that in fact $\p {\mathbf{p}}l \sim l$ in the RRN case and $\p {\mathbf{p}}l \sim (\sigma_v^2 + \sigma_a^2)l$ in the FRN case.

What about $\p {\boldsymbol{\oldgamma}}l$? The middle chart of \[fig:theory\_tracks\_pratice\] shows that over time, $\p {\mathbf{e}}l = \p {\boldsymbol{\oldgamma}}l / \p {\mathbf{p}}l$ contracts toward the center of the interval $[0, 1]$, but from the looks of it, it is not clear whether there is a stable fixed point ${\mathbf{e}}^*$ of ${\mathbf{e}}$ or not. We prove that, in fact, [**all trajectories of ${\mathbf{e}}$ not starting at 1 do converge to a single fixed point, but only at a polynomial rate**]{}, in both the RRN and FRN cases (\[thm:p\_q\_linear\] and \[thm:full\_res\_l\_g\_recurr\]); we can even explicitly compute the fixed point and the rate of convergence: For FRN, there is a [**unique stable fixed point**]{} ${\mathbf{e}}^* < 1$ determined by the equation $${\mathbf{e}}^* = \f 1 {\sigma_v^2 + \sigma_a^2}[\sigma_v^2 \f 2 \pi \arcsin\lp {\mathbf{e}}^* \rp + \sigma_a^2],$$ and $|{\mathbf{e}}^* - \p {\mathbf{e}}l|$ decreases like $l^{-\delta^*}$, where $$\delta^* := 1 - \f 2 \pi \f 1 {\sqrt{1 - ({\mathbf{e}}^*)^2}} \f{\sigma_v^2 }{\sigma_v^2 + \sigma_a^2}.$$ Since ${\mathbf{e}}^* < 1$, ${\mathbf{s}}= (1 - {\mathbf{e}}) {\mathbf{p}}= \Theta({\mathbf{p}}) = \Theta(l).$ The case of RRN can be viewed as a special case of the above, setting $\sigma_v^2 = 1$ and $\sigma_a^2 = 0$, which yields ${\mathbf{e}}^* = 0$ and $\delta^* = 1 - \f 2 \pi$. We observe that both ${\mathbf{e}}^*$ and $\delta^*$ only depend on the ratio $\rho := \sigma_a/\sigma_v$, so in \[fig:edelta\_plot\] we graph these two quantities as a function of $\rho$. ${\mathbf{e}}^*$ and $\delta^*$ both increase with $\rho$ and asymptotically approach 1 and $\nicefrac 1 2$ respectively from below. When $\rho = \sigma_a = 0$, ${\mathbf{e}}^* = 0$ and $\delta^* = 1 - \f 2 \pi$. Thus the rate of convergence at its [**slowest**]{} for tanh/FRN is $\delta^* = 1 - \f 2 \pi \approx 0.36338$, where asymptotically the network tends toward a [**chaotic regime**]{} ${\mathbf{e}}^* = 0$, corresponding to a large weight variance and a small bias variance; it at its [**fastest**]{} is $\delta^* = \nicefrac 1 2$, where asymptotically the network tends toward a [**stable regime**]{} ${\mathbf{e}}^* = 1$, corresponding to a large bias variance and small weight variance. We verify $\delta^*$ by comparing $\p {\mathbf{e}}l - \p {\mathbf{e}}{l-1}$ to $l^{-\delta^* - 1}$ in log-log scale. If $\p {\mathbf{e}}l = \Theta(l^{-\delta^*})$, then $\p {\mathbf{e}}l - \p {\mathbf{e}}{l-1} = \Theta(l^{-\delta^* - 1})$ and should obtain the same slope as $l^{-\delta^* - 1}$ as $l \to \infty$. The middle figure of \[fig:edelta\_plot\] ascertains that this is indeed the case, starting around layer number 400.

![Left-to-right: [**(a)**]{} Plots of ${\mathbf{e}}^*$ and $\delta^*$ against $\sigma_a/\sigma_v$. [**(b)**]{} In log-log scale: the dashed line is $l^{-\delta^* - 1}$, and the colored lines are $\p {\mathbf{e}}l - \p {\mathbf{e}}{l-1}$ for different initial conditions $\p {\mathbf{e}}0$. That they become parallel at about $l = 400$ on verifies that $\p {\mathbf{e}}l = \Theta(l^{-\delta^*})$. [^5] [**(c)**]{} In log-log scale: The dashed line is $\mathcal A \sqrt l$ ($\mathcal A$ given in \[thm:dalethExpSqrtTanhFullRes\]), and the colored lines are $\log(\p \bullet 1/\p \bullet l)$ for $\bullet = {{\boldsymbol{\oldchi}}}, {\boldsymbol{\oldchi}}_b, {\boldsymbol{\oldchi}}_w$. That they all converge together starting around $l=1000$ indicates that the approximation in \[thm:dalethExpSqrtTanhFullRes\] is very good for large $l$.[]{data-label="fig:edelta_plot"}](graphics/edelta_plot.pdf "fig:"){height=".16\textheight" width=".3\textwidth"} ![Left-to-right: [**(a)**]{} Plots of ${\mathbf{e}}^*$ and $\delta^*$ against $\sigma_a/\sigma_v$. [**(b)**]{} In log-log scale: the dashed line is $l^{-\delta^* - 1}$, and the colored lines are $\p {\mathbf{e}}l - \p {\mathbf{e}}{l-1}$ for different initial conditions $\p {\mathbf{e}}0$. That they become parallel at about $l = 400$ on verifies that $\p {\mathbf{e}}l = \Theta(l^{-\delta^*})$. [^6] [**(c)**]{} In log-log scale: The dashed line is $\mathcal A \sqrt l$ ($\mathcal A$ given in \[thm:dalethExpSqrtTanhFullRes\]), and the colored lines are $\log(\p \bullet 1/\p \bullet l)$ for $\bullet = {{\boldsymbol{\oldchi}}}, {\boldsymbol{\oldchi}}_b, {\boldsymbol{\oldchi}}_w$. That they all converge together starting around $l=1000$ indicates that the approximation in \[thm:dalethExpSqrtTanhFullRes\] is very good for large $l$.[]{data-label="fig:edelta_plot"}](graphics/verify_deltastar_n_grad-chi.pdf "fig:"){height=".16\textheight"}

#### Backward dynamics.

Finally, we show that the gradient is approximated by $$\begin{aligned}
\p {{\boldsymbol{\oldchi}}}{m} &= \exp(\mathcal A(\sqrt{l} - \sqrt{m}) + O(\log l - \log m))\p {{\boldsymbol{\oldchi}}}{l} \label{eqn:tanhGradEst}\tag{$\star$}\end{aligned}$$ where $\mathcal A = \f 4 3 \sqrt{\f 2 \pi} \sigma_w$ in the RRN case and $
\mathcal A = \f 4 3 \sqrt{\f 2 \pi} \f{\sigma_v^2 \sigma_w}{\sqrt{\sigma_v^2 + \sigma_a^2}}$ in the FRN case (\[thm:dalethExpSqrtTanh\] and \[thm:dalethExpSqrtTanhFullRes\]). The rightmost plot of \[fig:edelta\_plot\] verifies that indeed, for large $l \ge 1000$, this is a very good approximation. This demonstrates that the mean field assumption of independent backpropagation weights is very practical and convenient even for residual networks.

Note that in the FRN case, the constant $\mathcal A$ can be decomposed into $\mathcal A = \f 4 3 \sqrt{\f 2 \pi} \cdot \sigma_v \cdot \sigma_w \cdot (1 + \sigma_a^2/\sigma_v^2)^{-1/2}$. Consider the ratio $\rho := \sigma_a/\sigma_v$. If $\rho \gg 1$, then ${\mathbf{e}}^* \approx 1$ (\[fig:jjj\_vs\_id\_main\]), meaning that the typical network essentially computes a constant function, and thus unexpressive; at the same time, large $\rho$ makes $\mathcal A$ small, and thus ameliorating the gradient explosion problem, making the network more trainable. On the other hand, if $\rho \ll 1$, then ${\mathbf{e}}^* \approx 0$ (\[fig:jjj\_vs\_id\_main\]), the typical network can tease out the finest differences between any two input vectors, and a final linear layer on top of such a network should be able to express a wide variety of functions [@poole_exponential_2016]; at the same time, small $\rho$ increases $\mathcal A$, worsening the gradient explosion problem, making the network less trainable. This is the same expressivity-trainability tradeoff discussed in [@schoenholz_deep_2017].

$\alpha$-ReLU
-------------

#### Forward dynamics.

As with the tanh case, to deduce the asymptotic behavior of random $\alpha$-ReLU resnets, we need to understand the transforms ${\mathrm{V}}\psi_\alpha$ and ${\mathrm{W}}\psi_\alpha$. Fortunately, ${\mathrm{V}}\psi_\alpha$ has a closed form, and ${\mathrm{W}}\psi_\alpha$ has been studied before [@cho_kernel_2009]. In particular, if $\alpha > -\f 1 2$, then ${\mathrm{V}}\psi_\alpha( {\mathbf{q}}) = {\mathsf{c}}_\alpha {\mathbf{q}}^{\alpha}$, where ${\mathsf{c}}_\alpha$ is a constant with a closed form given by \[lemma:VtPsiAlpha\]. In addition, by [@cho_kernel_2009], we know that ${\mathrm{W}}\psi_\alpha( {\mathbf{q}}, {\mathbf{c}}{\mathbf{q}}) = {\mathrm{V}}\psi_\alpha( {\mathbf{q}}) {\mathbb{J}}_\alpha({\mathbf{c}})$ for ${\mathbb{J}}_\alpha$ given in \[sec:AlphaReluForwardProofs\]. \[fig:jjj\_vs\_id\_main\] shows a comparison of ${\mathbb{J}}_\alpha$ for different $\alpha$s along with the identity function.

Substituting in ${\mathsf{c}}_\alpha {\mathbf{q}}^\alpha$ for ${\mathrm{V}}\psi_\alpha$, we get a difference equation ${\mathbf{p}}- {\underline}{\mathbf{p}}= \sigma_v^2 {\mathsf{c}}_\alpha (\sigma_w^2 {\underline}{\mathbf{p}}+ \sigma_b^2)^\alpha + \sigma_a^2$ governing the evolution of ${\mathbf{p}}$. This should be reminiscent of the differential equation $\dot P(l) = C P(l)^\alpha$, which has solution $\propto l^{1/(1-\alpha)}$ for $\alpha < 1$, and $\propto \exp(Cl)$ when $\alpha = 1$. And indeed, the solutions $\p {\mathbf{p}}l$ to these difference equations behave asymptotically exactly like so (\[thm:pDynamicAlphaReLU\]). Thus [**ReLU behaves very explosively compared to $\alpha$-ReLU with $\alpha<1$**]{}. In fact, in simulations, for $\sigma_w^2 = 1.69$ and $\sigma_v^2 = 1.5$, the ReLU resnets overflows into `inf`s after around 100 layers, while there’s no problem from any other kind of networks we consider.

Regardless, [**$\alpha$-ReLU for all $\alpha$ massages $\p {\mathbf{e}}l$ toward a fixed point ${\mathbf{e}}^*$ that depends on $\alpha$**]{}. [ When $\phi = \psi_1$, the standard ReLU, $\p {\mathbf{e}}l$ converges to 1 asymptotically as $C l^{-2}$ for an explicit constant $C$ depending on $\sigma_v$ and $\sigma_w$ only (\[thm:ReLUSquaredConvergence\]), so that ${\mathbf{s}}= (1 - {\mathbf{e}}){\mathbf{p}}= \Theta(l^{-2}\exp(\Theta(l))) = \exp(\Theta(l)).$ When $\phi = \psi_\alpha$ for $\alpha < 1$, then $\p {\mathbf{e}}l$ converges to the nonunit fixed point ${\mathbf{e}}^*$ of ${\mathbb{J}}_\alpha$ at a rate of ${{\check\Theta}}(l^{-\mu})$, where $\mu = (1-\dot {\mathbb{J}}_\alpha({\mathbf{e}}^*))/(1-\alpha)$ is independent of the variances (\[thm:alphaReLUeConvergence\]), so that ${\mathbf{s}}= \Theta({\mathbf{p}})$.]{} These rates are verified in \[fig:alphaReLUVerifyExponents\].

#### Backward dynamics.

Finally, we have also characterized the rate of gradient growth for any $\alpha \in (\f 3 4, 1]$. [^7] [**In the case of $\alpha = 1$, the dynamics of ${{\boldsymbol{\oldchi}}}$ is exponential**]{}, the same as that of ${\mathbf{p}}$, $\p{{\boldsymbol{\oldchi}}}{l-m} = \p {{\boldsymbol{\oldchi}}}{l} B^m$ where $B =\f 1 2 \sigma_v^2 \sigma_w^2 + 1$. [**For $\alpha \in (\f 3 4, 1)$, the dynamics is polynomial**]{}, but with different exponent in general from that of the forward pass: $\p{{\boldsymbol{\oldchi}}}{l-m} = \Theta(1) \p {{\boldsymbol{\oldchi}}}{l} (l/(l-m))^R$ for $R = \f{\alpha^2}{(1-\alpha)(2 \alpha - 1)}$, where the constants in $\Theta(1)$ do not depend on $l$ or $m$. This exponent $R$ is minimized on $\alpha \in [\f 3 4, 1)$ at $\alpha = \nicefrac 3 4$, where $R = \nicefrac 9 2$ (but on $\alpha \in (\f 1 2, 1)$ it is minimized at $\alpha = \nicefrac 2 3$, where $R = 4$); see \[fig:backprop\_exponent\_alpha-relu\]. These exponents are verified empirically in \[fig:alphaReLUVerifyExponents\].

Looking only at ${{\boldsymbol{\oldchi}}}$ and the gradients against the biases, it seems that ReLU suffers from a dramatic case of exploding gradients. But in fact, because ${{\boldsymbol{\oldchi}}}$ gains a factor of $B$ moving backwards while ${\mathbf{p}}$ loses a factor of $B$, the gradient norm ${\boldsymbol{\oldchi}}_w^{(l-m)}$ (and similarly for ${\boldsymbol{\oldchi}}_v^{(l-m)}$) is independent of how far, $m$, the gradient has been propagated (\[thm:alphaReLUAllGradients\]) — this is certainly the best gradient preservation among all of the models considered in this paper. Thus strangely, random ReLU FRN exhibits both the best (constant for $v$ and $w$) and the worse (exponential for $a$ and $b$) gradient dynamics. This begs the question, then, is this a better deal than other $\alpha$-ReLU for which for any learnable parameter we have at most a polynomial blowup with depth in its gradient? Our experiments (discussed below) show that $\alpha$-ReLU is useful to the extent that smaller $\alpha$ avoids numerical issues with exponentiating forward and backward dynamics, but the best performance is given by the largest $\alpha$ that avoids them (\[fig:tanhHeatmaps\](c, d)); in fact, the metric expressivity ${\mathbf{s}}$, determines performance, not gradient explosion (see $\alpha$-ReLU experiments).

Experimental Results {#sec:experiments}
====================

![From left to right, top to bottom: **(a)** and **(b)**: $\sigma_w^2$, $L$, and test set accuracy of a grid of tanh reduced (left) and full (right) resnets trained on MNIST. Color indicates performance, with ligher colors indicating higher accuracy on test set. Other than the values on the axes, we have fixed $\sigma_b^2 = \sigma_a^2 = \f 1 2$ and $\sigma_v^2 = 1$. The white dotted lines are given by $\sigma_w^2 L = C$, where $C = 170$ on the left and $C = 145$ on the right. We see that both dotted lines accurately predict the largest optimal $\sigma_w$ for each depth $L$. **(c)** Varying the ratio $\sigma_a^2/\sigma_v^2$ while fixing $\sigma_v/\sqrt{1 + \sigma_a^2/\sigma_v^2}$, and thus fixing $\mathcal A$, the leading constant of $\log \p {{\boldsymbol{\oldchi}}}0 / \p {{\boldsymbol{\oldchi}}}L$. **(d)** in log-log scale: Heatmap gives the test accuracies of ReLU FRN for varying $\sigma_w^2$ and $L$. Curves give level sets for the log ratios $\log \p {\mathbf{s}}L / \p {\mathbf{s}}0 \approx \log \p {\mathbf{p}}L / \p {\mathbf{p}}0 \approx \log \p{{\boldsymbol{\oldchi}}}0 / \p {{\boldsymbol{\oldchi}}}L = L \log(1 + \sigma_v^2 \sigma_w^2/2)$. **(e)** Red heatmap shows the test accuracies of a grid of $\alpha$-ReLU FRN with varying $\alpha$ and $L$ as shown, but with all $\sigma_\bullet$s fixed. The white dashed curve gives a typical contour line of $L^R = \text{const}$, where $R = \f {\alpha^2}{(1-\alpha)(2\alpha-1)}.$ The yellow-to-blue curves form a set of level curves for $\p {\mathbf{s}}l = \p {\mathbf{p}}l - \p {\boldsymbol{\oldgamma}}l = \text{const}$, with yellow curves corresponding to higher levels. []{data-label="fig:tanhHeatmaps"}](graphics/tanh_simple_resnet_trainability "fig:"){height="\height"} ![From left to right, top to bottom: **(a)** and **(b)**: $\sigma_w^2$, $L$, and test set accuracy of a grid of tanh reduced (left) and full (right) resnets trained on MNIST. Color indicates performance, with ligher colors indicating higher accuracy on test set. Other than the values on the axes, we have fixed $\sigma_b^2 = \sigma_a^2 = \f 1 2$ and $\sigma_v^2 = 1$. The white dotted lines are given by $\sigma_w^2 L = C$, where $C = 170$ on the left and $C = 145$ on the right. We see that both dotted lines accurately predict the largest optimal $\sigma_w$ for each depth $L$. **(c)** Varying the ratio $\sigma_a^2/\sigma_v^2$ while fixing $\sigma_v/\sqrt{1 + \sigma_a^2/\sigma_v^2}$, and thus fixing $\mathcal A$, the leading constant of $\log \p {{\boldsymbol{\oldchi}}}0 / \p {{\boldsymbol{\oldchi}}}L$. **(d)** in log-log scale: Heatmap gives the test accuracies of ReLU FRN for varying $\sigma_w^2$ and $L$. Curves give level sets for the log ratios $\log \p {\mathbf{s}}L / \p {\mathbf{s}}0 \approx \log \p {\mathbf{p}}L / \p {\mathbf{p}}0 \approx \log \p{{\boldsymbol{\oldchi}}}0 / \p {{\boldsymbol{\oldchi}}}L = L \log(1 + \sigma_v^2 \sigma_w^2/2)$. **(e)** Red heatmap shows the test accuracies of a grid of $\alpha$-ReLU FRN with varying $\alpha$ and $L$ as shown, but with all $\sigma_\bullet$s fixed. The white dashed curve gives a typical contour line of $L^R = \text{const}$, where $R = \f {\alpha^2}{(1-\alpha)(2\alpha-1)}.$ The yellow-to-blue curves form a set of level curves for $\p {\mathbf{s}}l = \p {\mathbf{p}}l - \p {\boldsymbol{\oldgamma}}l = \text{const}$, with yellow curves corresponding to higher levels. []{data-label="fig:tanhHeatmaps"}](graphics/tanh_resnet_trainability "fig:"){height="\height"} ![From left to right, top to bottom: **(a)** and **(b)**: $\sigma_w^2$, $L$, and test set accuracy of a grid of tanh reduced (left) and full (right) resnets trained on MNIST. Color indicates performance, with ligher colors indicating higher accuracy on test set. Other than the values on the axes, we have fixed $\sigma_b^2 = \sigma_a^2 = \f 1 2$ and $\sigma_v^2 = 1$. The white dotted lines are given by $\sigma_w^2 L = C$, where $C = 170$ on the left and $C = 145$ on the right. We see that both dotted lines accurately predict the largest optimal $\sigma_w$ for each depth $L$. **(c)** Varying the ratio $\sigma_a^2/\sigma_v^2$ while fixing $\sigma_v/\sqrt{1 + \sigma_a^2/\sigma_v^2}$, and thus fixing $\mathcal A$, the leading constant of $\log \p {{\boldsymbol{\oldchi}}}0 / \p {{\boldsymbol{\oldchi}}}L$. **(d)** in log-log scale: Heatmap gives the test accuracies of ReLU FRN for varying $\sigma_w^2$ and $L$. Curves give level sets for the log ratios $\log \p {\mathbf{s}}L / \p {\mathbf{s}}0 \approx \log \p {\mathbf{p}}L / \p {\mathbf{p}}0 \approx \log \p{{\boldsymbol{\oldchi}}}0 / \p {{\boldsymbol{\oldchi}}}L = L \log(1 + \sigma_v^2 \sigma_w^2/2)$. **(e)** Red heatmap shows the test accuracies of a grid of $\alpha$-ReLU FRN with varying $\alpha$ and $L$ as shown, but with all $\sigma_\bullet$s fixed. The white dashed curve gives a typical contour line of $L^R = \text{const}$, where $R = \f {\alpha^2}{(1-\alpha)(2\alpha-1)}.$ The yellow-to-blue curves form a set of level curves for $\p {\mathbf{s}}l = \p {\mathbf{p}}l - \p {\boldsymbol{\oldgamma}}l = \text{const}$, with yellow curves corresponding to higher levels. []{data-label="fig:tanhHeatmaps"}](graphics/tanh_resnet_sa_sb "fig:"){height="\height"}\
![From left to right, top to bottom: **(a)** and **(b)**: $\sigma_w^2$, $L$, and test set accuracy of a grid of tanh reduced (left) and full (right) resnets trained on MNIST. Color indicates performance, with ligher colors indicating higher accuracy on test set. Other than the values on the axes, we have fixed $\sigma_b^2 = \sigma_a^2 = \f 1 2$ and $\sigma_v^2 = 1$. The white dotted lines are given by $\sigma_w^2 L = C$, where $C = 170$ on the left and $C = 145$ on the right. We see that both dotted lines accurately predict the largest optimal $\sigma_w$ for each depth $L$. **(c)** Varying the ratio $\sigma_a^2/\sigma_v^2$ while fixing $\sigma_v/\sqrt{1 + \sigma_a^2/\sigma_v^2}$, and thus fixing $\mathcal A$, the leading constant of $\log \p {{\boldsymbol{\oldchi}}}0 / \p {{\boldsymbol{\oldchi}}}L$. **(d)** in log-log scale: Heatmap gives the test accuracies of ReLU FRN for varying $\sigma_w^2$ and $L$. Curves give level sets for the log ratios $\log \p {\mathbf{s}}L / \p {\mathbf{s}}0 \approx \log \p {\mathbf{p}}L / \p {\mathbf{p}}0 \approx \log \p{{\boldsymbol{\oldchi}}}0 / \p {{\boldsymbol{\oldchi}}}L = L \log(1 + \sigma_v^2 \sigma_w^2/2)$. **(e)** Red heatmap shows the test accuracies of a grid of $\alpha$-ReLU FRN with varying $\alpha$ and $L$ as shown, but with all $\sigma_\bullet$s fixed. The white dashed curve gives a typical contour line of $L^R = \text{const}$, where $R = \f {\alpha^2}{(1-\alpha)(2\alpha-1)}.$ The yellow-to-blue curves form a set of level curves for $\p {\mathbf{s}}l = \p {\mathbf{p}}l - \p {\boldsymbol{\oldgamma}}l = \text{const}$, with yellow curves corresponding to higher levels. []{data-label="fig:tanhHeatmaps"}](graphics/reluHeatmapContour "fig:"){height="\height"} ![From left to right, top to bottom: **(a)** and **(b)**: $\sigma_w^2$, $L$, and test set accuracy of a grid of tanh reduced (left) and full (right) resnets trained on MNIST. Color indicates performance, with ligher colors indicating higher accuracy on test set. Other than the values on the axes, we have fixed $\sigma_b^2 = \sigma_a^2 = \f 1 2$ and $\sigma_v^2 = 1$. The white dotted lines are given by $\sigma_w^2 L = C$, where $C = 170$ on the left and $C = 145$ on the right. We see that both dotted lines accurately predict the largest optimal $\sigma_w$ for each depth $L$. **(c)** Varying the ratio $\sigma_a^2/\sigma_v^2$ while fixing $\sigma_v/\sqrt{1 + \sigma_a^2/\sigma_v^2}$, and thus fixing $\mathcal A$, the leading constant of $\log \p {{\boldsymbol{\oldchi}}}0 / \p {{\boldsymbol{\oldchi}}}L$. **(d)** in log-log scale: Heatmap gives the test accuracies of ReLU FRN for varying $\sigma_w^2$ and $L$. Curves give level sets for the log ratios $\log \p {\mathbf{s}}L / \p {\mathbf{s}}0 \approx \log \p {\mathbf{p}}L / \p {\mathbf{p}}0 \approx \log \p{{\boldsymbol{\oldchi}}}0 / \p {{\boldsymbol{\oldchi}}}L = L \log(1 + \sigma_v^2 \sigma_w^2/2)$. **(e)** Red heatmap shows the test accuracies of a grid of $\alpha$-ReLU FRN with varying $\alpha$ and $L$ as shown, but with all $\sigma_\bullet$s fixed. The white dashed curve gives a typical contour line of $L^R = \text{const}$, where $R = \f {\alpha^2}{(1-\alpha)(2\alpha-1)}.$ The yellow-to-blue curves form a set of level curves for $\p {\mathbf{s}}l = \p {\mathbf{p}}l - \p {\boldsymbol{\oldgamma}}l = \text{const}$, with yellow curves corresponding to higher levels. []{data-label="fig:tanhHeatmaps"}](graphics/alphaReluGrid_GammaGapOverlay "fig:"){height="\height"}

Our experiments show a dichotomy of what matters in initialization: for tanh resnets, quality of an initialization is determined by how much gradient explosion there is (measured by $\p {{\boldsymbol{\oldchi}}}0 / \p {{\boldsymbol{\oldchi}}}L$); for ($\alpha$-)ReLU resnets, it is determined by how expressive the random network is (measured by the metric expressivity $\p {\mathbf{s}}L$). We hypothesize this is because in tanh resnets, the gradient dynamics is much more explosive than the expressivity dynamics ($\exp(\Theta(\sqrt l))$ vs $\Theta(l)$), whereas for ReLU it’s somewhat the opposite (${\boldsymbol{\oldchi}}_w, {\boldsymbol{\oldchi}}_v = \Theta(1)$ vs ${\mathbf{s}}= \exp(\Theta(l))$).

#### Tanh, vary $\sigma_w$.

We train a grid of reduced and full tanh resnets on MNIST, varying the variance $\sigma_w^2$ and the number of layers (for FRN we fix $\sigma_v = 1$). The results are indicated in \[fig:tanhHeatmaps\](a, b). We see that in either model, deeper resnets favor much smaller $\sigma_w$ than shallower ones. The white dotted lines in \[fig:tanhHeatmaps\](a, b) confirm our theory: according to \[eqn:tanhGradEst\], for the same gradient ratio $R = \p {{\boldsymbol{\oldchi}}}0 / \p {{\boldsymbol{\oldchi}}}L$, we want $\log R \approx \sigma_w \sqrt L$. Indeed, the white dotted lines in \[fig:tanhHeatmaps\](a, b) trace out such a level curve and it remarkably pinpoints the largest $\sigma_w$ that gives the optimal test set accuracy for each depth $L$. Why isn’t the best initialization given by $R = 1 \iff \sigma_w = 0$? We believe that when $L$ and/or $\sigma_w$ is small, gradient dynamics no longer dominates the initialization quality because it has “less room to explode,” and expressivity issues start to dampen the test time performance.

#### Tanh, vary $\sigma_a^2/\sigma_v^2$.

As suggested in the analysis of \[eqn:tanhGradEst\], the ratio $\rho^2 = \sigma_a^2/\sigma_v^2$ determines the fixed point ${\mathbf{e}}^*$ and its convergence rate by itself while also contributes to the rate of gradient explosion in tanh FRN. We seek to isolate its effect on forward dynamics by varying $\sigma_v$ with $\rho$ such that $\sigma_v/\sqrt{1 + \rho^2}$ is kept constant, so that the leading term of the log gradient ratio is kept approximately equal for each $L$ and $\rho$. \[fig:tanhHeatmaps\](c) shows the test accuracies of a grid of tanh FRN initialized with such an ensemble of $\sigma_\bullet$s. What stands out the most is that performance is maximized essentially around a fixed value of $L$ regardless of $\rho$, which shows that indeed gradient dynamics determines the initialization quality in tanh resnets. There is also a minor increase in performance with increasing $\rho$ regardless of $L$; this is counterintuitive as increasing $\rho$ means “decreasing expressivity.” It is currently not clear what accounts for this effect.

#### ReLU, vary $\sigma_w$

We train a grid of ReLU FRN on MNIST, varying $\sigma_w^2 \in [0, 1.5]$ while fixing $\sigma_v^2 = 1, \sigma_a^2 = \sigma_b^2 = \f 1 2$. The resulting test set accuracies are shown in \[fig:tanhHeatmaps\](d). The dark upper region signifies failure of training caused by numerical issues with exploding activation and gradient norms: This corresponds to the region where $\p {\mathbf{p}}L$, which is a measure of the mean magnitude of an neuronal activation in layer $L$, becomes too big. We see that the best test accuracies are given by depths just below where these numerical issues occur. However, if we were to predict that the optimal init is the one minimizing $\p {\boldsymbol{\oldchi}}0 / \p {\boldsymbol{\oldchi}}L \ge 1$, then we would be wrong — in fact it is exactly the opposite. In this case, the dynamics of $\p {\mathbf{s}}l, \p {\mathbf{p}}l$, and $\p {\boldsymbol{\oldchi}}0 / \p {\boldsymbol{\oldchi}}l$ are approximately the same (all $\exp(\Theta(l))$ with the same hidden constants), and optimal performance corresponds to the highest $\p {\mathbf{s}}L$, $\p {\mathbf{p}}L$, and $\p {\boldsymbol{\oldchi}}0 / \p {\boldsymbol{\oldchi}}L$ without running into `inf`s.

#### $\alpha$-ReLU, vary $\alpha$.

We similarly trained a grid of $\alpha$-ReLU FRN on MNIST, varying only $\alpha$ and the depth, fixing all $\sigma_\bullet$. \[fig:tanhHeatmaps\](e) shows their test accuracies. We see similar behavior to ReLU, where when the net is too deep, numerical issues doom the training (black upper right corner), but the best performance is given by $L$ just below where this problem occurs. In this case, if we were to predict optimality based on minimizing gradient explosion, we would be again wrong, and furthermore, the contour plot of $\p {\boldsymbol{\oldchi}}0 / \p {\boldsymbol{\oldchi}}L$ (white dashed line) now gives no information at all on the test set accuracy. In contrast, the contours for $\p {\mathbf{s}}l$ succeeds remarkably well at this prediction (yellow/green lines).[^8] By interpolation, this suggests that indeed in the ReLU case, it is expressivity, not trainability, which determines performance at test time.

In all of our experiments, we did not find ${\mathbf{e}}$ dynamics to be predictive of neural network performance.

Conclusion
==========

In this paper, we have extended the mean field formalism developed by [@poole_exponential_2016; @raghu_expressive_2016; @schoenholz_deep_2017] to residual networks, a class of models closer to practice than classical feedforward neural networks as were investigated earlier. We proved and verified that in both the forward and backward passes, most of the residual networks discussed here do not collapse their input space geometry or the gradient information exponentially. We found our theory incredibly predictive of test time performance despite saying nothing about the dynamics of training. In addition, we overwhelmingly find, through theory and experiments, that an optimal initialization scheme must take into account the depth of the residual network. The reason that Xavier [@glorot_understanding_2010] or He [@he_delving_2015] scheme are not the best for residual networks is in fact not that their statistical assumptions are fragile — theirs are similar to our mean field theoretic assumptions, and they hold up in experiments for large width — but rather that their structural assumptions on the network break very badly on residual nets.

#### Open Problems.

Our work thus have shown that optimality of initialization schemes can be very unstable with respect to architecture. We hope this work will form a foundation toward a mathematically grounded initialization scheme for state-of-the-art architectures like the original He et al. residual network. To do so, there are still two major components left to study out of the following three:

Residual/skip connection

Batchnorm

Convolutional layers.

Recurrent architectures and attention mechanisms are also still mostly unexplored in terms of mean field theory. Furthermore, many theoretical questions still yet to be resolved; the most important with regard to mean field theory is: why can we make \[ass:symAct,ass:gradInd\] and still be able to make accurate predictions? We hope to make progress on these problems in the future and encourage readers to take part in this effort.

Acknowledgments {#acknowledgments .unnumbered}
===============

Thanks to Jeffrey Ling for early exploration experiments and help with the initial draft. Thanks to Felix Wong for offering his wisdom and experience working in statistical physics.

Additional Figures
==================

In figures appearing in the appendix, $\olddaleth$ means ${\boldsymbol{\oldchi}}$ (due to legacy reasons).

[m[.1]{}m[0.6]{}]{} $\alpha=1$ & ![Empirical vs theoretical dynamics for $\p {\mathbf{p}}l, \p {\mathbf{e}}l$, and different gradient quantities for $\alpha$-ReLU, with format similar to \[fig:theory\_tracks\_pratice\]. We refer to each figure on each row from left to right as (a), (b), and (c). Note that in the $\alpha=1$ case, figure (a) ($\p {\mathbf{p}}l$ and $\p {\boldsymbol{\oldgamma}}l$ for different initial values) has log scale y-axis and (a) and (b) have x-axis ranging from 1 to 50, while for other $\alpha$, (a) has normal y-axis and (a) and (b) have x-axis ranging from 1 to 200. We do so because the norm of the activation vector in a typical ReLU resnet blows up into `NaN` at around layer 90, while this is not a problem for $\alpha < 1$. Our theoretical predictions track the average of empirical values closely for forward quantities $\p {\mathbf{p}}l, \p {\boldsymbol{\oldgamma}}l,$ and $\p {\mathbf{e}}\l$ for all $\alpha$, but variance is extremely large for $\p {\mathbf{e}}l$ at $\alpha = 1$; it also predicts the average gradient norm accurately for $\alpha = 1$ to $\alpha = .7$ (despite the fact that we should not expect so for $\alpha \le .75$ due to exploding variance (\[thm:dalethInfVarAlphaReLU\])), although variance is large for $\alpha = 1$ at earlier layers (i.e. later layers w.r.t backpropagation). However it [*consistently and significantly overestimates*]{} the average gradient norm for $\alpha = .6$ to $\alpha = .5$, where the variance is so large that one standard deviation below the mean results in negative values. All plots are made with parameters $\sigma_v^2 = 1.5, \sigma_a^2 = .5, \sigma_w^2 = 1.69, \sigma_b^2 = .49$; only $\alpha$ is varied. All figures exhibit smooth curves, which are theoretical estimates, and irregular curves with shades around them, which indicate empirical means and standard deviations (both of which taken in regular scale, not log scale). For each $\alpha$, figures (a) and (b) are made with 20 runs of resnets of width 1000. (c) is made with 25 runs of resnets of width 250.[]{data-label="fig:alphaReLUTheoryVsEmpirics"}](graphics/relu_full_res_theory_vs_exp "fig:"){height=".08\textheight"}\
$\alpha=.9$ & ![Empirical vs theoretical dynamics for $\p {\mathbf{p}}l, \p {\mathbf{e}}l$, and different gradient quantities for $\alpha$-ReLU, with format similar to \[fig:theory\_tracks\_pratice\]. We refer to each figure on each row from left to right as (a), (b), and (c). Note that in the $\alpha=1$ case, figure (a) ($\p {\mathbf{p}}l$ and $\p {\boldsymbol{\oldgamma}}l$ for different initial values) has log scale y-axis and (a) and (b) have x-axis ranging from 1 to 50, while for other $\alpha$, (a) has normal y-axis and (a) and (b) have x-axis ranging from 1 to 200. We do so because the norm of the activation vector in a typical ReLU resnet blows up into `NaN` at around layer 90, while this is not a problem for $\alpha < 1$. Our theoretical predictions track the average of empirical values closely for forward quantities $\p {\mathbf{p}}l, \p {\boldsymbol{\oldgamma}}l,$ and $\p {\mathbf{e}}\l$ for all $\alpha$, but variance is extremely large for $\p {\mathbf{e}}l$ at $\alpha = 1$; it also predicts the average gradient norm accurately for $\alpha = 1$ to $\alpha = .7$ (despite the fact that we should not expect so for $\alpha \le .75$ due to exploding variance (\[thm:dalethInfVarAlphaReLU\])), although variance is large for $\alpha = 1$ at earlier layers (i.e. later layers w.r.t backpropagation). However it [*consistently and significantly overestimates*]{} the average gradient norm for $\alpha = .6$ to $\alpha = .5$, where the variance is so large that one standard deviation below the mean results in negative values. All plots are made with parameters $\sigma_v^2 = 1.5, \sigma_a^2 = .5, \sigma_w^2 = 1.69, \sigma_b^2 = .49$; only $\alpha$ is varied. All figures exhibit smooth curves, which are theoretical estimates, and irregular curves with shades around them, which indicate empirical means and standard deviations (both of which taken in regular scale, not log scale). For each $\alpha$, figures (a) and (b) are made with 20 runs of resnets of width 1000. (c) is made with 25 runs of resnets of width 250.[]{data-label="fig:alphaReLUTheoryVsEmpirics"}](graphics/9relu_full_res_theory_vs_exp "fig:"){height=".08\textheight"}\
$\alpha=.8$ & ![Empirical vs theoretical dynamics for $\p {\mathbf{p}}l, \p {\mathbf{e}}l$, and different gradient quantities for $\alpha$-ReLU, with format similar to \[fig:theory\_tracks\_pratice\]. We refer to each figure on each row from left to right as (a), (b), and (c). Note that in the $\alpha=1$ case, figure (a) ($\p {\mathbf{p}}l$ and $\p {\boldsymbol{\oldgamma}}l$ for different initial values) has log scale y-axis and (a) and (b) have x-axis ranging from 1 to 50, while for other $\alpha$, (a) has normal y-axis and (a) and (b) have x-axis ranging from 1 to 200. We do so because the norm of the activation vector in a typical ReLU resnet blows up into `NaN` at around layer 90, while this is not a problem for $\alpha < 1$. Our theoretical predictions track the average of empirical values closely for forward quantities $\p {\mathbf{p}}l, \p {\boldsymbol{\oldgamma}}l,$ and $\p {\mathbf{e}}\l$ for all $\alpha$, but variance is extremely large for $\p {\mathbf{e}}l$ at $\alpha = 1$; it also predicts the average gradient norm accurately for $\alpha = 1$ to $\alpha = .7$ (despite the fact that we should not expect so for $\alpha \le .75$ due to exploding variance (\[thm:dalethInfVarAlphaReLU\])), although variance is large for $\alpha = 1$ at earlier layers (i.e. later layers w.r.t backpropagation). However it [*consistently and significantly overestimates*]{} the average gradient norm for $\alpha = .6$ to $\alpha = .5$, where the variance is so large that one standard deviation below the mean results in negative values. All plots are made with parameters $\sigma_v^2 = 1.5, \sigma_a^2 = .5, \sigma_w^2 = 1.69, \sigma_b^2 = .49$; only $\alpha$ is varied. All figures exhibit smooth curves, which are theoretical estimates, and irregular curves with shades around them, which indicate empirical means and standard deviations (both of which taken in regular scale, not log scale). For each $\alpha$, figures (a) and (b) are made with 20 runs of resnets of width 1000. (c) is made with 25 runs of resnets of width 250.[]{data-label="fig:alphaReLUTheoryVsEmpirics"}](graphics/8relu_full_res_theory_vs_exp "fig:"){height=".08\textheight"}\
$\alpha=.7$ & ![Empirical vs theoretical dynamics for $\p {\mathbf{p}}l, \p {\mathbf{e}}l$, and different gradient quantities for $\alpha$-ReLU, with format similar to \[fig:theory\_tracks\_pratice\]. We refer to each figure on each row from left to right as (a), (b), and (c). Note that in the $\alpha=1$ case, figure (a) ($\p {\mathbf{p}}l$ and $\p {\boldsymbol{\oldgamma}}l$ for different initial values) has log scale y-axis and (a) and (b) have x-axis ranging from 1 to 50, while for other $\alpha$, (a) has normal y-axis and (a) and (b) have x-axis ranging from 1 to 200. We do so because the norm of the activation vector in a typical ReLU resnet blows up into `NaN` at around layer 90, while this is not a problem for $\alpha < 1$. Our theoretical predictions track the average of empirical values closely for forward quantities $\p {\mathbf{p}}l, \p {\boldsymbol{\oldgamma}}l,$ and $\p {\mathbf{e}}\l$ for all $\alpha$, but variance is extremely large for $\p {\mathbf{e}}l$ at $\alpha = 1$; it also predicts the average gradient norm accurately for $\alpha = 1$ to $\alpha = .7$ (despite the fact that we should not expect so for $\alpha \le .75$ due to exploding variance (\[thm:dalethInfVarAlphaReLU\])), although variance is large for $\alpha = 1$ at earlier layers (i.e. later layers w.r.t backpropagation). However it [*consistently and significantly overestimates*]{} the average gradient norm for $\alpha = .6$ to $\alpha = .5$, where the variance is so large that one standard deviation below the mean results in negative values. All plots are made with parameters $\sigma_v^2 = 1.5, \sigma_a^2 = .5, \sigma_w^2 = 1.69, \sigma_b^2 = .49$; only $\alpha$ is varied. All figures exhibit smooth curves, which are theoretical estimates, and irregular curves with shades around them, which indicate empirical means and standard deviations (both of which taken in regular scale, not log scale). For each $\alpha$, figures (a) and (b) are made with 20 runs of resnets of width 1000. (c) is made with 25 runs of resnets of width 250.[]{data-label="fig:alphaReLUTheoryVsEmpirics"}](graphics/7relu_full_res_theory_vs_exp "fig:"){height=".08\textheight"}\
$\alpha=.6$ & ![Empirical vs theoretical dynamics for $\p {\mathbf{p}}l, \p {\mathbf{e}}l$, and different gradient quantities for $\alpha$-ReLU, with format similar to \[fig:theory\_tracks\_pratice\]. We refer to each figure on each row from left to right as (a), (b), and (c). Note that in the $\alpha=1$ case, figure (a) ($\p {\mathbf{p}}l$ and $\p {\boldsymbol{\oldgamma}}l$ for different initial values) has log scale y-axis and (a) and (b) have x-axis ranging from 1 to 50, while for other $\alpha$, (a) has normal y-axis and (a) and (b) have x-axis ranging from 1 to 200. We do so because the norm of the activation vector in a typical ReLU resnet blows up into `NaN` at around layer 90, while this is not a problem for $\alpha < 1$. Our theoretical predictions track the average of empirical values closely for forward quantities $\p {\mathbf{p}}l, \p {\boldsymbol{\oldgamma}}l,$ and $\p {\mathbf{e}}\l$ for all $\alpha$, but variance is extremely large for $\p {\mathbf{e}}l$ at $\alpha = 1$; it also predicts the average gradient norm accurately for $\alpha = 1$ to $\alpha = .7$ (despite the fact that we should not expect so for $\alpha \le .75$ due to exploding variance (\[thm:dalethInfVarAlphaReLU\])), although variance is large for $\alpha = 1$ at earlier layers (i.e. later layers w.r.t backpropagation). However it [*consistently and significantly overestimates*]{} the average gradient norm for $\alpha = .6$ to $\alpha = .5$, where the variance is so large that one standard deviation below the mean results in negative values. All plots are made with parameters $\sigma_v^2 = 1.5, \sigma_a^2 = .5, \sigma_w^2 = 1.69, \sigma_b^2 = .49$; only $\alpha$ is varied. All figures exhibit smooth curves, which are theoretical estimates, and irregular curves with shades around them, which indicate empirical means and standard deviations (both of which taken in regular scale, not log scale). For each $\alpha$, figures (a) and (b) are made with 20 runs of resnets of width 1000. (c) is made with 25 runs of resnets of width 250.[]{data-label="fig:alphaReLUTheoryVsEmpirics"}](graphics/6relu_full_res_theory_vs_exp "fig:"){height=".08\textheight"}\
$\alpha=.55$ & ![Empirical vs theoretical dynamics for $\p {\mathbf{p}}l, \p {\mathbf{e}}l$, and different gradient quantities for $\alpha$-ReLU, with format similar to \[fig:theory\_tracks\_pratice\]. We refer to each figure on each row from left to right as (a), (b), and (c). Note that in the $\alpha=1$ case, figure (a) ($\p {\mathbf{p}}l$ and $\p {\boldsymbol{\oldgamma}}l$ for different initial values) has log scale y-axis and (a) and (b) have x-axis ranging from 1 to 50, while for other $\alpha$, (a) has normal y-axis and (a) and (b) have x-axis ranging from 1 to 200. We do so because the norm of the activation vector in a typical ReLU resnet blows up into `NaN` at around layer 90, while this is not a problem for $\alpha < 1$. Our theoretical predictions track the average of empirical values closely for forward quantities $\p {\mathbf{p}}l, \p {\boldsymbol{\oldgamma}}l,$ and $\p {\mathbf{e}}\l$ for all $\alpha$, but variance is extremely large for $\p {\mathbf{e}}l$ at $\alpha = 1$; it also predicts the average gradient norm accurately for $\alpha = 1$ to $\alpha = .7$ (despite the fact that we should not expect so for $\alpha \le .75$ due to exploding variance (\[thm:dalethInfVarAlphaReLU\])), although variance is large for $\alpha = 1$ at earlier layers (i.e. later layers w.r.t backpropagation). However it [*consistently and significantly overestimates*]{} the average gradient norm for $\alpha = .6$ to $\alpha = .5$, where the variance is so large that one standard deviation below the mean results in negative values. All plots are made with parameters $\sigma_v^2 = 1.5, \sigma_a^2 = .5, \sigma_w^2 = 1.69, \sigma_b^2 = .49$; only $\alpha$ is varied. All figures exhibit smooth curves, which are theoretical estimates, and irregular curves with shades around them, which indicate empirical means and standard deviations (both of which taken in regular scale, not log scale). For each $\alpha$, figures (a) and (b) are made with 20 runs of resnets of width 1000. (c) is made with 25 runs of resnets of width 250.[]{data-label="fig:alphaReLUTheoryVsEmpirics"}](graphics/55relu_full_res_theory_vs_exp "fig:"){height=".08\textheight"}\
$\alpha=.51$ & ![Empirical vs theoretical dynamics for $\p {\mathbf{p}}l, \p {\mathbf{e}}l$, and different gradient quantities for $\alpha$-ReLU, with format similar to \[fig:theory\_tracks\_pratice\]. We refer to each figure on each row from left to right as (a), (b), and (c). Note that in the $\alpha=1$ case, figure (a) ($\p {\mathbf{p}}l$ and $\p {\boldsymbol{\oldgamma}}l$ for different initial values) has log scale y-axis and (a) and (b) have x-axis ranging from 1 to 50, while for other $\alpha$, (a) has normal y-axis and (a) and (b) have x-axis ranging from 1 to 200. We do so because the norm of the activation vector in a typical ReLU resnet blows up into `NaN` at around layer 90, while this is not a problem for $\alpha < 1$. Our theoretical predictions track the average of empirical values closely for forward quantities $\p {\mathbf{p}}l, \p {\boldsymbol{\oldgamma}}l,$ and $\p {\mathbf{e}}\l$ for all $\alpha$, but variance is extremely large for $\p {\mathbf{e}}l$ at $\alpha = 1$; it also predicts the average gradient norm accurately for $\alpha = 1$ to $\alpha = .7$ (despite the fact that we should not expect so for $\alpha \le .75$ due to exploding variance (\[thm:dalethInfVarAlphaReLU\])), although variance is large for $\alpha = 1$ at earlier layers (i.e. later layers w.r.t backpropagation). However it [*consistently and significantly overestimates*]{} the average gradient norm for $\alpha = .6$ to $\alpha = .5$, where the variance is so large that one standard deviation below the mean results in negative values. All plots are made with parameters $\sigma_v^2 = 1.5, \sigma_a^2 = .5, \sigma_w^2 = 1.69, \sigma_b^2 = .49$; only $\alpha$ is varied. All figures exhibit smooth curves, which are theoretical estimates, and irregular curves with shades around them, which indicate empirical means and standard deviations (both of which taken in regular scale, not log scale). For each $\alpha$, figures (a) and (b) are made with 20 runs of resnets of width 1000. (c) is made with 25 runs of resnets of width 250.[]{data-label="fig:alphaReLUTheoryVsEmpirics"}](graphics/51relu_full_res_theory_vs_exp "fig:"){height=".08\textheight"}

[m[.1]{}m[0.4]{}]{} $\alpha=.9$ & ![We verify the exponents of the forward and backward dynamics for $\alpha$-ReLU FRN. For each row, the figures are labeled (a) and (b) from left to right. The format is the same as in \[fig:jjj\_vs\_id\_main\]. All figures are in log-log scale. [**(a)**]{} We exhibit our theoretical dynamics of the cosine distance $\p {\mathbf{e}}l$ based on the recurrences \[thm:fullResPQRec\] and \[thm:full\_res\_l\_g\_recurr\] for different initial conditions $\p {\mathbf{e}}0$. We draw $|\p {\mathbf{e}}l - \p {\mathbf{e}}{l-1}|$ for each of these dynamics in colored solid lines. We predict that each dynamic is ${{\check\Theta}}(l^{-\mu})$, where $\mu = (1 - \dot{\mathbb{J}}_\alpha({\mathbf{e}}^*))/(1-\alpha)$, and the dashed line gives $l^{-\mu-1}$ (\[thm:alphaReLUeConvergence\]), shifted vertically to better compare the slope in log scale (i.e. the exponent of the polynomial dynamics). (See footnote \[footnote:plotDelta\] for why we plot the dynamics this way). We see that the our asymptotic prediction is very accurate for the sequence of $\p {\mathbf{e}}l$ that starts with $\p {\mathbf{e}}0 = 0.99$, the closest to ${\mathbf{e}}^*$ for each $\alpha$, while other lines only slowly converge to the same exponent (which is the slope in the log-log plot). This is to be expected based on the proof of \[thm:alphaReLUeConvergence\]. For $\alpha = .9$, the $\p {\mathbf{e}}0 = .99$ line upticks at around $10^3$ and then turn into `NaN`s due to numerical instability. [**(b)**]{} Colored lines are $\p \bullet 0 / \p \bullet l$ for $\bullet = {{\boldsymbol{\oldchi}}}, {\boldsymbol{\oldchi}}_b, {\boldsymbol{\oldchi}}_w$ (we are not taking logs in addition to plotting in log-log scale like in \[fig:jjj\_vs\_id\]). The dashed lines are our asymptotic predictions for the dynamics with corresponding colors, based on \[thm:alphaReLUAllGradients\], again shifted appropriately to easily compare slope visually. We see that for every alpha our asymptotic predictions are highly accurate. For both (a) and (b), we did not show $\alpha = 1$ case as ReLU FRN runs into numerical issues quickly (i.e. with even for 100 layers) because of exponential explosions in $\p {\mathbf{p}}l$ and $\p {{\boldsymbol{\oldchi}}}l$ as predicted by \[thm:pDynamicAlphaReLU,thm:dalethDynamicsAlphaReLU\], so we cannot expect to empirically verify the precise predicted asymptotics. All plots are made with parameters $\sigma_v^2 = 1.5, \sigma_a^2 = .5, \sigma_w^2 = 1.69, \sigma_b^2 = .49$; only $\alpha$ is varied. []{data-label="fig:alphaReLUVerifyExponents"}](graphics/9relu_verify_deltastar_n_grad "fig:"){height=".08\textheight"}\
$\alpha=.8$ & ![We verify the exponents of the forward and backward dynamics for $\alpha$-ReLU FRN. For each row, the figures are labeled (a) and (b) from left to right. The format is the same as in \[fig:jjj\_vs\_id\_main\]. All figures are in log-log scale. [**(a)**]{} We exhibit our theoretical dynamics of the cosine distance $\p {\mathbf{e}}l$ based on the recurrences \[thm:fullResPQRec\] and \[thm:full\_res\_l\_g\_recurr\] for different initial conditions $\p {\mathbf{e}}0$. We draw $|\p {\mathbf{e}}l - \p {\mathbf{e}}{l-1}|$ for each of these dynamics in colored solid lines. We predict that each dynamic is ${{\check\Theta}}(l^{-\mu})$, where $\mu = (1 - \dot{\mathbb{J}}_\alpha({\mathbf{e}}^*))/(1-\alpha)$, and the dashed line gives $l^{-\mu-1}$ (\[thm:alphaReLUeConvergence\]), shifted vertically to better compare the slope in log scale (i.e. the exponent of the polynomial dynamics). (See footnote \[footnote:plotDelta\] for why we plot the dynamics this way). We see that the our asymptotic prediction is very accurate for the sequence of $\p {\mathbf{e}}l$ that starts with $\p {\mathbf{e}}0 = 0.99$, the closest to ${\mathbf{e}}^*$ for each $\alpha$, while other lines only slowly converge to the same exponent (which is the slope in the log-log plot). This is to be expected based on the proof of \[thm:alphaReLUeConvergence\]. For $\alpha = .9$, the $\p {\mathbf{e}}0 = .99$ line upticks at around $10^3$ and then turn into `NaN`s due to numerical instability. [**(b)**]{} Colored lines are $\p \bullet 0 / \p \bullet l$ for $\bullet = {{\boldsymbol{\oldchi}}}, {\boldsymbol{\oldchi}}_b, {\boldsymbol{\oldchi}}_w$ (we are not taking logs in addition to plotting in log-log scale like in \[fig:jjj\_vs\_id\]). The dashed lines are our asymptotic predictions for the dynamics with corresponding colors, based on \[thm:alphaReLUAllGradients\], again shifted appropriately to easily compare slope visually. We see that for every alpha our asymptotic predictions are highly accurate. For both (a) and (b), we did not show $\alpha = 1$ case as ReLU FRN runs into numerical issues quickly (i.e. with even for 100 layers) because of exponential explosions in $\p {\mathbf{p}}l$ and $\p {{\boldsymbol{\oldchi}}}l$ as predicted by \[thm:pDynamicAlphaReLU,thm:dalethDynamicsAlphaReLU\], so we cannot expect to empirically verify the precise predicted asymptotics. All plots are made with parameters $\sigma_v^2 = 1.5, \sigma_a^2 = .5, \sigma_w^2 = 1.69, \sigma_b^2 = .49$; only $\alpha$ is varied. []{data-label="fig:alphaReLUVerifyExponents"}](graphics/8relu_verify_deltastar_n_grad "fig:"){height=".08\textheight"}\
$\alpha=.7$ & ![We verify the exponents of the forward and backward dynamics for $\alpha$-ReLU FRN. For each row, the figures are labeled (a) and (b) from left to right. The format is the same as in \[fig:jjj\_vs\_id\_main\]. All figures are in log-log scale. [**(a)**]{} We exhibit our theoretical dynamics of the cosine distance $\p {\mathbf{e}}l$ based on the recurrences \[thm:fullResPQRec\] and \[thm:full\_res\_l\_g\_recurr\] for different initial conditions $\p {\mathbf{e}}0$. We draw $|\p {\mathbf{e}}l - \p {\mathbf{e}}{l-1}|$ for each of these dynamics in colored solid lines. We predict that each dynamic is ${{\check\Theta}}(l^{-\mu})$, where $\mu = (1 - \dot{\mathbb{J}}_\alpha({\mathbf{e}}^*))/(1-\alpha)$, and the dashed line gives $l^{-\mu-1}$ (\[thm:alphaReLUeConvergence\]), shifted vertically to better compare the slope in log scale (i.e. the exponent of the polynomial dynamics). (See footnote \[footnote:plotDelta\] for why we plot the dynamics this way). We see that the our asymptotic prediction is very accurate for the sequence of $\p {\mathbf{e}}l$ that starts with $\p {\mathbf{e}}0 = 0.99$, the closest to ${\mathbf{e}}^*$ for each $\alpha$, while other lines only slowly converge to the same exponent (which is the slope in the log-log plot). This is to be expected based on the proof of \[thm:alphaReLUeConvergence\]. For $\alpha = .9$, the $\p {\mathbf{e}}0 = .99$ line upticks at around $10^3$ and then turn into `NaN`s due to numerical instability. [**(b)**]{} Colored lines are $\p \bullet 0 / \p \bullet l$ for $\bullet = {{\boldsymbol{\oldchi}}}, {\boldsymbol{\oldchi}}_b, {\boldsymbol{\oldchi}}_w$ (we are not taking logs in addition to plotting in log-log scale like in \[fig:jjj\_vs\_id\]). The dashed lines are our asymptotic predictions for the dynamics with corresponding colors, based on \[thm:alphaReLUAllGradients\], again shifted appropriately to easily compare slope visually. We see that for every alpha our asymptotic predictions are highly accurate. For both (a) and (b), we did not show $\alpha = 1$ case as ReLU FRN runs into numerical issues quickly (i.e. with even for 100 layers) because of exponential explosions in $\p {\mathbf{p}}l$ and $\p {{\boldsymbol{\oldchi}}}l$ as predicted by \[thm:pDynamicAlphaReLU,thm:dalethDynamicsAlphaReLU\], so we cannot expect to empirically verify the precise predicted asymptotics. All plots are made with parameters $\sigma_v^2 = 1.5, \sigma_a^2 = .5, \sigma_w^2 = 1.69, \sigma_b^2 = .49$; only $\alpha$ is varied. []{data-label="fig:alphaReLUVerifyExponents"}](graphics/7relu_verify_deltastar_n_grad "fig:"){height=".08\textheight"}\
$\alpha=.6$ & ![We verify the exponents of the forward and backward dynamics for $\alpha$-ReLU FRN. For each row, the figures are labeled (a) and (b) from left to right. The format is the same as in \[fig:jjj\_vs\_id\_main\]. All figures are in log-log scale. [**(a)**]{} We exhibit our theoretical dynamics of the cosine distance $\p {\mathbf{e}}l$ based on the recurrences \[thm:fullResPQRec\] and \[thm:full\_res\_l\_g\_recurr\] for different initial conditions $\p {\mathbf{e}}0$. We draw $|\p {\mathbf{e}}l - \p {\mathbf{e}}{l-1}|$ for each of these dynamics in colored solid lines. We predict that each dynamic is ${{\check\Theta}}(l^{-\mu})$, where $\mu = (1 - \dot{\mathbb{J}}_\alpha({\mathbf{e}}^*))/(1-\alpha)$, and the dashed line gives $l^{-\mu-1}$ (\[thm:alphaReLUeConvergence\]), shifted vertically to better compare the slope in log scale (i.e. the exponent of the polynomial dynamics). (See footnote \[footnote:plotDelta\] for why we plot the dynamics this way). We see that the our asymptotic prediction is very accurate for the sequence of $\p {\mathbf{e}}l$ that starts with $\p {\mathbf{e}}0 = 0.99$, the closest to ${\mathbf{e}}^*$ for each $\alpha$, while other lines only slowly converge to the same exponent (which is the slope in the log-log plot). This is to be expected based on the proof of \[thm:alphaReLUeConvergence\]. For $\alpha = .9$, the $\p {\mathbf{e}}0 = .99$ line upticks at around $10^3$ and then turn into `NaN`s due to numerical instability. [**(b)**]{} Colored lines are $\p \bullet 0 / \p \bullet l$ for $\bullet = {{\boldsymbol{\oldchi}}}, {\boldsymbol{\oldchi}}_b, {\boldsymbol{\oldchi}}_w$ (we are not taking logs in addition to plotting in log-log scale like in \[fig:jjj\_vs\_id\]). The dashed lines are our asymptotic predictions for the dynamics with corresponding colors, based on \[thm:alphaReLUAllGradients\], again shifted appropriately to easily compare slope visually. We see that for every alpha our asymptotic predictions are highly accurate. For both (a) and (b), we did not show $\alpha = 1$ case as ReLU FRN runs into numerical issues quickly (i.e. with even for 100 layers) because of exponential explosions in $\p {\mathbf{p}}l$ and $\p {{\boldsymbol{\oldchi}}}l$ as predicted by \[thm:pDynamicAlphaReLU,thm:dalethDynamicsAlphaReLU\], so we cannot expect to empirically verify the precise predicted asymptotics. All plots are made with parameters $\sigma_v^2 = 1.5, \sigma_a^2 = .5, \sigma_w^2 = 1.69, \sigma_b^2 = .49$; only $\alpha$ is varied. []{data-label="fig:alphaReLUVerifyExponents"}](graphics/6relu_verify_deltastar_n_grad "fig:"){height=".08\textheight"}\
$\alpha=.55$ & ![We verify the exponents of the forward and backward dynamics for $\alpha$-ReLU FRN. For each row, the figures are labeled (a) and (b) from left to right. The format is the same as in \[fig:jjj\_vs\_id\_main\]. All figures are in log-log scale. [**(a)**]{} We exhibit our theoretical dynamics of the cosine distance $\p {\mathbf{e}}l$ based on the recurrences \[thm:fullResPQRec\] and \[thm:full\_res\_l\_g\_recurr\] for different initial conditions $\p {\mathbf{e}}0$. We draw $|\p {\mathbf{e}}l - \p {\mathbf{e}}{l-1}|$ for each of these dynamics in colored solid lines. We predict that each dynamic is ${{\check\Theta}}(l^{-\mu})$, where $\mu = (1 - \dot{\mathbb{J}}_\alpha({\mathbf{e}}^*))/(1-\alpha)$, and the dashed line gives $l^{-\mu-1}$ (\[thm:alphaReLUeConvergence\]), shifted vertically to better compare the slope in log scale (i.e. the exponent of the polynomial dynamics). (See footnote \[footnote:plotDelta\] for why we plot the dynamics this way). We see that the our asymptotic prediction is very accurate for the sequence of $\p {\mathbf{e}}l$ that starts with $\p {\mathbf{e}}0 = 0.99$, the closest to ${\mathbf{e}}^*$ for each $\alpha$, while other lines only slowly converge to the same exponent (which is the slope in the log-log plot). This is to be expected based on the proof of \[thm:alphaReLUeConvergence\]. For $\alpha = .9$, the $\p {\mathbf{e}}0 = .99$ line upticks at around $10^3$ and then turn into `NaN`s due to numerical instability. [**(b)**]{} Colored lines are $\p \bullet 0 / \p \bullet l$ for $\bullet = {{\boldsymbol{\oldchi}}}, {\boldsymbol{\oldchi}}_b, {\boldsymbol{\oldchi}}_w$ (we are not taking logs in addition to plotting in log-log scale like in \[fig:jjj\_vs\_id\]). The dashed lines are our asymptotic predictions for the dynamics with corresponding colors, based on \[thm:alphaReLUAllGradients\], again shifted appropriately to easily compare slope visually. We see that for every alpha our asymptotic predictions are highly accurate. For both (a) and (b), we did not show $\alpha = 1$ case as ReLU FRN runs into numerical issues quickly (i.e. with even for 100 layers) because of exponential explosions in $\p {\mathbf{p}}l$ and $\p {{\boldsymbol{\oldchi}}}l$ as predicted by \[thm:pDynamicAlphaReLU,thm:dalethDynamicsAlphaReLU\], so we cannot expect to empirically verify the precise predicted asymptotics. All plots are made with parameters $\sigma_v^2 = 1.5, \sigma_a^2 = .5, \sigma_w^2 = 1.69, \sigma_b^2 = .49$; only $\alpha$ is varied. []{data-label="fig:alphaReLUVerifyExponents"}](graphics/55relu_verify_deltastar_n_grad "fig:"){height=".08\textheight"}\
$\alpha=.51$ & ![We verify the exponents of the forward and backward dynamics for $\alpha$-ReLU FRN. For each row, the figures are labeled (a) and (b) from left to right. The format is the same as in \[fig:jjj\_vs\_id\_main\]. All figures are in log-log scale. [**(a)**]{} We exhibit our theoretical dynamics of the cosine distance $\p {\mathbf{e}}l$ based on the recurrences \[thm:fullResPQRec\] and \[thm:full\_res\_l\_g\_recurr\] for different initial conditions $\p {\mathbf{e}}0$. We draw $|\p {\mathbf{e}}l - \p {\mathbf{e}}{l-1}|$ for each of these dynamics in colored solid lines. We predict that each dynamic is ${{\check\Theta}}(l^{-\mu})$, where $\mu = (1 - \dot{\mathbb{J}}_\alpha({\mathbf{e}}^*))/(1-\alpha)$, and the dashed line gives $l^{-\mu-1}$ (\[thm:alphaReLUeConvergence\]), shifted vertically to better compare the slope in log scale (i.e. the exponent of the polynomial dynamics). (See footnote \[footnote:plotDelta\] for why we plot the dynamics this way). We see that the our asymptotic prediction is very accurate for the sequence of $\p {\mathbf{e}}l$ that starts with $\p {\mathbf{e}}0 = 0.99$, the closest to ${\mathbf{e}}^*$ for each $\alpha$, while other lines only slowly converge to the same exponent (which is the slope in the log-log plot). This is to be expected based on the proof of \[thm:alphaReLUeConvergence\]. For $\alpha = .9$, the $\p {\mathbf{e}}0 = .99$ line upticks at around $10^3$ and then turn into `NaN`s due to numerical instability. [**(b)**]{} Colored lines are $\p \bullet 0 / \p \bullet l$ for $\bullet = {{\boldsymbol{\oldchi}}}, {\boldsymbol{\oldchi}}_b, {\boldsymbol{\oldchi}}_w$ (we are not taking logs in addition to plotting in log-log scale like in \[fig:jjj\_vs\_id\]). The dashed lines are our asymptotic predictions for the dynamics with corresponding colors, based on \[thm:alphaReLUAllGradients\], again shifted appropriately to easily compare slope visually. We see that for every alpha our asymptotic predictions are highly accurate. For both (a) and (b), we did not show $\alpha = 1$ case as ReLU FRN runs into numerical issues quickly (i.e. with even for 100 layers) because of exponential explosions in $\p {\mathbf{p}}l$ and $\p {{\boldsymbol{\oldchi}}}l$ as predicted by \[thm:pDynamicAlphaReLU,thm:dalethDynamicsAlphaReLU\], so we cannot expect to empirically verify the precise predicted asymptotics. All plots are made with parameters $\sigma_v^2 = 1.5, \sigma_a^2 = .5, \sigma_w^2 = 1.69, \sigma_b^2 = .49$; only $\alpha$ is varied. []{data-label="fig:alphaReLUVerifyExponents"}](graphics/51relu_verify_deltastar_n_grad "fig:"){height=".08\textheight"}\

  Symbol                                    Meaning                                                                                              Ref
  ----------------------------------------- ---------------------------------------------------------------------------------------------------- ----------------------------------
  $\sigma_\bullet$                          standard deviation of trainable parameter $\bullet$                                                  
  $\p x l$                                  activation vector/input vector                                                                       
  $\p h l$                                  hidden vector                                                                                        
  $N$                                       width (same across all layers)                                                                       
  $\p {\mathbf{p}}l$                        m.n. squared length of activation vector $\p x l$                                                    \[defn:length\]
  $\p {\mathbf{q}}l$                        m.n. squared length of hidden vector $\p h l$                                                        \[defn:length\]
  $\p {\boldsymbol{\oldgamma}}l$            m.n. dot product $\p x l \cdot \p x l{}'$                                                            \[defn:corr\]
  $\p {\boldsymbol{\oldlambda}}l$           m.n. dot product $\p h l \cdot \p h l{}'$                                                            \[defn:corr\]
  $\p {\mathbf{s}}l$                        m.n. squared distance $\|\p x l - \p x l {}'\|^2$                                                    \[defn:corr\]
  $\p {\mathbf{e}}l$                        cosine distance $\p {\boldsymbol{\oldgamma}}l / \sqrt{\p {\mathbf{p}}l \p {\mathbf{p}}l {}'}$        \[defn:corr\]
  ${\mathbf{e}}^*$                          limit value of $\p {\mathbf{e}}l$ as $l \to \infty$                                                  
  $\p {\mathbf{c}}l$                        cosine distance $\p {\boldsymbol{\oldlambda}}l / \sqrt{\p {\mathbf{q}}l \p {\mathbf{q}}l {}'}$       \[defn:corr\]
  $\p {\boldsymbol{\oldchi}}l$              m.n. gradient squared norm w.r.t. $\p x l$                                                           \[defn:grad\]
  $\p {{\boldsymbol{\oldchi}}_\bullet} l$   m.n. gradient squared norm w.r.t. trainable parameter $\bullet$                                      \[defn:grad\]
  $\phi$                                    variable nonlinearity $\R \to \R$                                                                    
  $\psi_\alpha$                             $\alpha$-ReLU                                                                                        \[defn:alphaReLU\]
  ${\mathrm{V}}$                            variance integral transform                                                                          \[defn:integralTransform\]
  ${\mathrm{W}}$                            covariance integral transform                                                                        \[defn:integralTransform\]
  $\delta^*$                                $\p {\mathbf{e}}l$ converges like $\Theta(l^{-\delta^*})$ in tanh FRN                                \[thm:eDynamicsFullResTanh\]
  $\mathcal A$                              leading coeff of $\log \p {\boldsymbol{\oldchi}}0 / \p {\boldsymbol{\oldchi}}L$ in tanh FRN          \[thm:dalethExpSqrtTanhFullRes\]
  $R$                                       $\log \p {\boldsymbol{\oldchi}}0 / \p {\boldsymbol{\oldchi}}L \sim R \log L$ for $(\alpha<1)$-ReLU   \[thm:dalethDynamicsAlphaReLU\]
  ${\mathbb{J}}_\alpha$                     kernel function of $\alpha$-ReLU                                                                     \[lemma:basicJalpha\]

  : Glossary of Symbols. “Mean normalized” is abbreviated “m.n.”[]{data-label="tab:glossary"}

A Listing of Main Theorems
==========================

Tanh
----

### Reduced Residual Network

[lemma]{}[pqrecurrence]{}\[lemma:p\_q\_recurrence\] Suppose $\phi$ is antisymmetric. Then in an RRN, ${\mathbf{p}}$ and ${\mathbf{q}}$ satisfy the recurrence $$\begin{aligned}
    {\mathbf{q}}&= \sigma_w^2 {\underline}{\mathbf{p}}+ \sigma_b^2\\
    {\mathbf{p}}&= {\mathrm{V}}\phi( {\mathbf{q}}) + {\underline}{\mathbf{p}}.\end{aligned}$$

[thm]{}[pqlinear]{}\[thm:p\_q\_linear\] Suppose $\phi$ is tanh-like. Assume RRN architecture.

-   If $\sigma_w = 0$, then $\p {\mathbf{p}}l = l {\mathrm{V}}\phi( \sigma_b^2) + \p {\mathbf{p}}0$ and $\p {\mathbf{q}}l = \sigma_b^2$.

-   If $\sigma_w > 0$, $\lim_{l \to \infty} \p {\mathbf{p}}l/ l = 1$ and $\lim_{l \to \infty} \p {\mathbf{q}}l /(\sigma_w^2 l) = 1$. If $\phi = \tanh$, then we can obtain more terms of the asymptotic expansions: $$\begin{aligned}
        \p {\mathbf{p}}l &= l - 2 C \sigma_w^{-1} l^{1/2} - C^2 \sigma_w^{-2} \log l + O(1)\\
        \p {\mathbf{q}}l &= \sigma_w^2 l - 2 C \sigma_w l^{1/2} - C^2 \log l + O(1)
        \end{aligned}$$ as $l \to \infty$, where $C = \sqrt{2 / \pi}$.

[thm]{}[lambdagammarecurrence]{}\[thm:lambda\_gamma\_recurrence\] Suppose $\phi$ is antisymmetric. Then in an RRN, ${\boldsymbol{\oldlambda}}$ and ${\boldsymbol{\oldgamma}}$ satisfy the recurrence $$\begin{aligned}
    {\boldsymbol{\oldlambda}}&= \sigma_w^2 {\underline}{\boldsymbol{\oldgamma}}+ \sigma_b^2\\
    {\boldsymbol{\oldgamma}}&= {\mathrm{W}}\phi({\mathbf{q}}, {\boldsymbol{\oldlambda}}) + {\underline}{\boldsymbol{\oldgamma}}.\end{aligned}$$

[thm]{}[edynamics]{}\[thm:edynamics\] Suppose $\phi$ is a tanh-like nonlinearity in an RRN. Assume $\p {\mathbf{e}}0 < 1$.

-   If $\sigma_w = 0$, then $\p {\boldsymbol{\oldgamma}}l = l {\mathrm{W}}\phi( \sigma_b^2, \sigma_b^2) + \p {\boldsymbol{\oldgamma}}0 = l {\mathrm{V}}\phi( \sigma_b^2) + \p {\boldsymbol{\oldgamma}}0$ and $\p {\boldsymbol{\oldlambda}}l = \sigma_b^2$, so that $\p {\mathbf{e}}l \to 1$ and $1 - \p {\mathbf{e}}l = \Theta(l^{-1})$. As a result, $\p {\mathbf{s}}l = \p {\mathbf{p}}l (1 - \p {\mathbf{e}}l) = \Theta(1).$

-   If $\sigma_w > 0$, then $\p {\boldsymbol{\oldgamma}}l = {{\check\Theta}}(l^{\f 2\pi})$, and $\p {\mathbf{e}}l \to 0$ like ${{\check\Theta}}(l^{\f 2 \pi - 1})$. Thus $\p {\mathbf{s}}l = \Theta(\p {\mathbf{p}}l) = \Theta(l).$

[thm]{}[dalethRecReduced]{}\[thm:dalethRecReduced\] For any nonlinearity $\phi$ in an RRN, under assumptions \[ass:symAct\] and \[ass:gradInd\], whenever $\dot \phi^2(\zeta)$ has finite variance for Gaussian variable $\zeta$, $$\begin{aligned}
{\underline}{{\boldsymbol{\oldchi}}}&= (\sigma_w^2  {\mathrm{V}}\dot \phi( {\mathbf{q}}) + 1){{\boldsymbol{\oldchi}}},&
{\boldsymbol{\oldchi}}_b &= {{\boldsymbol{\oldchi}}}{\mathrm{V}}\dot \phi( {\mathbf{q}}),&
{\boldsymbol{\oldchi}}_w &= {{\boldsymbol{\oldchi}}}{\mathrm{V}}\dot \phi( {\mathbf{q}}) {\underline}{\mathbf{p}}. \end{aligned}$$

[thm]{}[dalethExpSqrtTanh]{}\[thm:dalethExpSqrtTanh\] For $\phi = \tanh$ in an RRN,

-   If $\sigma_w = 0$, $\p {{\boldsymbol{\oldchi}}}m = \p {{\boldsymbol{\oldchi}}}l$ for all $l, m$.

-   If $\sigma_w > 0$, $$\begin{aligned}
        \log(\p {{\boldsymbol{\oldchi}}}{m}/\p {{\boldsymbol{\oldchi}}}l) &= \mathcal A (\sqrt l - \sqrt m) + \mathcal B (\log l - \log m) + O(1)
        \end{aligned}$$ where $\mathcal A = \f 4 3 \sqrt{\f 2 \pi} \sigma_w$ and $\mathcal B = \f 4 {3\pi} - \sigma_w^2 \f 4 {9\pi}$.

[thm]{}[dalethExpSqrtTanhAllGrad]{}\[thm:dalethExpSqrtTanhAllGrad\] Suppose $\phi = \tanh$. Then in an RRN

-   If $\sigma_w = 0$, $\p {\boldsymbol{\oldchi}}l _b = \p {{\boldsymbol{\oldchi}}}L {\mathrm{V}}\dot \phi( \sigma_b^2)$ and $\p {\boldsymbol{\oldchi}}l _w = \p {{\boldsymbol{\oldchi}}}L {\mathrm{V}}\dot \phi( \sigma_b^2) ((l-1) {\mathrm{V}}\phi( \sigma_b^2) + \p {\mathbf{p}}0),$ where $L$ is the last layer.

-   If $\sigma_w > 0$, $$\begin{aligned}
        \log(\p{\boldsymbol{\oldchi}}m _b / \p {\boldsymbol{\oldchi}}l _b) &= \mathcal A (\sqrt l - \sqrt m) + \mathcal B_b (\log l - \log m) + O(1)\\
        \log(\p{\boldsymbol{\oldchi}}m _w / \p {\boldsymbol{\oldchi}}l _w) &= \mathcal A (\sqrt l - \sqrt m) + \mathcal B_w (\log l - \log m) + O(1)
        \end{aligned}$$ where $\mathcal A = \f 4 3 \sqrt{\f 2 \pi} \sigma_w$ (same as $\mathcal A$ in \[thm:dalethExpSqrtTanh\]) and $\mathcal B_b = \mathcal B + \f 1 2, \mathcal B_w = \mathcal B - \f 1 2$, with $\mathcal B = \f 4 {3\pi} - \sigma_w^2 \f 4 {9\pi}$ (same as $\mathcal B$ in \[thm:dalethExpSqrtTanh\]).

### Full Residual Network

[thm]{}[fullResPQRec]{} \[thm:fullResPQRec\] For any nonlinearity $\phi$ in an FRN, $$\begin{aligned}
{\mathbf{q}}&= \sigma_w^2 {\underline}{\mathbf{p}}+ \sigma_b^2\\
{\mathbf{p}}&= \sigma_v^2 {\mathrm{V}}\phi( {\mathbf{q}}) + \sigma_a^2 + {\underline}{\mathbf{p}}\end{aligned}$$

[thm]{}[pIsLinearTanh]{}\[thm:pIsLinearTanh\] Suppose $\phi$ is tanh-like. Assume the FRN architecture.

-   If $\sigma_w = 0$, then $\p {\mathbf{p}}l = (\sigma_v^2 {\mathrm{V}}\phi( \sigma_b^2) + \sigma_a^2)l +\p {\mathbf{p}}0$, and $\p {\mathbf{q}}l = \sigma_b^2$.

-   If $\sigma_w > 0$, then $\p {\mathbf{p}}l = b_0 l + b_1 l^{1/2} + b_2 \log l + O(1)$, where $$\begin{aligned}
        b_0 &= \sigma_v^2 + \sigma_a^2\\
        b_1 &= \f{-2C \sigma_v^2 \sigma_w^{-1}}{\sqrt{\sigma_v^2 + \sigma_a^2}}\\
        b_2 &= \f{-C^2 \sigma_v^4 \sigma_w^{-2}}{(\sigma_v^2 + \sigma_a^2)^2}
        \end{aligned}$$ and $C = \sqrt{\f 2 \pi}$. Additionally, $\p {\mathbf{q}}l = \sigma_w^2 b_0 l + \sigma_w^2 b_1 l^{1/2} + \sigma_w^2 b_2 \log l + O(1)$.

![Empirical verification of \[thm:pIsLinearTanh\].](graphics/tanh_p_asymptotic_expansion.pdf){height=".16\textheight"}

[thm]{}[LGRecFullRes]{}\[thm:full\_res\_l\_g\_recurr\] For any nonlinearity $\phi$, in an FRN $$\begin{aligned}
{\boldsymbol{\oldlambda}}&= \sigma_w^2 {\underline}{\boldsymbol{\oldgamma}}+ \sigma_b^2\\
{\boldsymbol{\oldgamma}}&= \sigma_v^2 {\mathrm{W}}\phi({\mathbf{q}}, {\boldsymbol{\oldlambda}}) + \sigma_a^2 + {\underline}{\boldsymbol{\oldgamma}}\end{aligned}$$

[thm]{}[eDynamicsFullResTanh]{} \[thm:eDynamicsFullResTanh\] Assume $\phi = \tanh$ in an FRN. Suppose $\p {\mathbf{e}}0 < 1$.

-   If $\sigma_w = 0$, then $\p {\boldsymbol{\oldlambda}}l = \sigma_b^2$ and $\p {\boldsymbol{\oldgamma}}l = l (\sigma_v^2 {\mathrm{W}}\phi( \sigma_b^2, \sigma_b^2) + \sigma_a^2) + \p {\boldsymbol{\oldgamma}}0 = l (\sigma_v^2 {\mathrm{V}}\phi( \sigma_b^2) + \sigma_a^2) + \p {\boldsymbol{\oldgamma}}0$. Thus $\p {\mathbf{e}}l \to 1$ and $1 - \p {\mathbf{e}}l = \Theta(l^{-1})$. As a result, $\p {\mathbf{s}}l = \p {\mathbf{p}}l (1 - \p {\mathbf{e}}l) = \Theta(1).$

-   If $\sigma_w > 0$, then $\p {\mathbf{e}}l$ converges to the unique fixed point ${\mathbf{e}}^* \not = 1$ determined by the equation $${\mathbf{e}}^* = \f 1 {\sigma_v^2 + \sigma_a^2}[\sigma_v^2 \f 2 \pi \arcsin\lp {\mathbf{e}}^* \rp + \sigma_a^2].$$ Furthermore, $\p {\mathbf{e}}l$ converges to ${\mathbf{e}}^*$ polynomially: $|\p {\mathbf{e}}l - {\mathbf{e}}^*|$ is ${{\check\Theta}}(l^{-\delta^*})$, where $$\delta^* := 1 - \f 2 \pi \f 1 {\sqrt{1 - ({\mathbf{e}}^*)^2}} \f{\sigma_v^2 }{\sigma_v^2 + \sigma_a^2} \in [\f 2 \pi - 1, \f 1 2)$$ Since ${\mathbf{e}}^* < 1$, $\p {\mathbf{s}}l = \Theta(\p {\mathbf{p}}l) = \Theta(l).$

[thm]{}[dalethRecFull]{} \[thm:dalethRecFull\] For any nonlinearity $\phi$ in an FRN, under assumptions \[ass:symAct\] and \[ass:gradInd\], whenever $\dot \phi(\zeta)^2$ has finite variance for Gaussian variable $\zeta$, $$\begin{aligned}
{\underline}{{\boldsymbol{\oldchi}}}&= (\sigma_v^2\sigma_w^2  {\mathrm{V}}\dot \phi( {\mathbf{q}}) + 1){{\boldsymbol{\oldchi}}},&
{\boldsymbol{\oldchi}}_b &= \sigma_v^2{{\boldsymbol{\oldchi}}}{\mathrm{V}}\dot \phi( {\mathbf{q}}),\\
{\boldsymbol{\oldchi}}_w &= \sigma_v^2{{\boldsymbol{\oldchi}}}{\mathrm{V}}\dot \phi( {\mathbf{q}}) {\underline}{\mathbf{p}},&
{\boldsymbol{\oldchi}}_v &= {{\boldsymbol{\oldchi}}}{\mathrm{V}}\phi( {\mathbf{q}}),&
{\boldsymbol{\oldchi}}_a &= {{\boldsymbol{\oldchi}}}\end{aligned}$$

[thm]{}[dalethExpSqrtTanhFullRes]{}\[thm:dalethExpSqrtTanhFullRes\] Assume $\phi = \tanh$ in an FRN.

-   If $\sigma_w = 0$, $\p {{\boldsymbol{\oldchi}}}m = \p {{\boldsymbol{\oldchi}}}l$ for all $l, m$.

-   If $\sigma_w > 0$, then for $l \ge m \ge 0,$ $$\log(\p {{\boldsymbol{\oldchi}}}{m} / \p {{\boldsymbol{\oldchi}}}l) = \mathcal A (\sqrt l - \sqrt m) + \mathcal B (\log l - \log m) + O(1)$$ where $$\begin{aligned}
        \mathcal A &= \f 4 3 \sqrt{\f 2 \pi} \f{\sigma_v^2 \sigma_w}{\sqrt{\sigma_v^2 + \sigma_a^2}}\\
        \mathcal B &= \f 4 {9\pi}\f{ \sigma_v^4 }{\sigma_v^2 + \sigma_a^2}\lp \f 3 {\sigma_v^2 + \sigma_a^2} - \sigma_w^2\rp
        \end{aligned}$$

\[fig:tanhasymptoticsgrid\] shows empirical verification of the asymptotic expansion of ${{\boldsymbol{\oldchi}}}$ for various values of $\sigma_\bullet$s.

![Empirical verification of the asymptotic expansion of ${{\boldsymbol{\oldchi}}}$ for various values of $\sigma_\bullet$s. Note that we have chosen all small values for $\sigma_\bullet$s. For larger values, the constant term in \[thm:dalethExpSqrtTanhFullRes\] begins to dominate (primarily because of the expansion $\log(1+x) = x + \Theta(x^2)$ has large $\Theta$ term when $x$ is large), and ${{\boldsymbol{\oldchi}}}$ behaves more like $\exp(\Theta(l))$ up to depth 1000.[]{data-label="fig:tanhasymptoticsgrid"}](graphics/tanhAsymptoticsGrid){height=".3\textheight"}

[thm]{}[dalethExpSqrtTanhFullResAllGrad]{}\[thm:dalethExpSqrtTanhFullResAllGrad\] Suppose $\phi = \tanh$ in an FRN.

-   If $\sigma_w = 0$, then $$\begin{aligned}
        \p {\boldsymbol{\oldchi}}l _b &= \sigma_v^2 \p {{\boldsymbol{\oldchi}}}L {\mathrm{V}}\dot \phi( \sigma_b^2)\\
        \p {\boldsymbol{\oldchi}}l _w &= \sigma_v^2 \p {{\boldsymbol{\oldchi}}}L {\mathrm{V}}\dot \phi( \sigma_b^2) ( (\sigma_v^2 {\mathrm{V}}\phi( \sigma_b^2) + \sigma_a^2)(l-1) +\p {\mathbf{p}}0)\\
        \p {\boldsymbol{\oldchi}}l _v &= \p {{\boldsymbol{\oldchi}}}L {\mathrm{V}}\phi( \sigma_b^2)\\
        \p {\boldsymbol{\oldchi}}l _a &= \p {{\boldsymbol{\oldchi}}}L.
        \end{aligned}$$

-   If $\sigma_w > 0$, then for $l \ge m \ge 0,$

    $$\begin{aligned}
        \log(\p{\boldsymbol{\oldchi}}m _b / \p {\boldsymbol{\oldchi}}l _b) &= \mathcal A (\sqrt l - \sqrt m) + \mathcal B_b (\log l - \log m) + O(1)\\
        \log(\p{\boldsymbol{\oldchi}}m _w / \p {\boldsymbol{\oldchi}}l _w) &= \mathcal A (\sqrt l - \sqrt m) + \mathcal B_w (\log l - \log m) + O(1)\\
        \log(\p{\boldsymbol{\oldchi}}m _a / \p {\boldsymbol{\oldchi}}l _a) &= \mathcal A (\sqrt l - \sqrt m) + \mathcal B (\log l - \log m) + O(1)\\
        \log(\p{\boldsymbol{\oldchi}}m _v / \p {\boldsymbol{\oldchi}}l _v) &= \mathcal A (\sqrt l - \sqrt m) + \mathcal B (\log l - \log m) + O(1)
        \end{aligned}$$

    where $\mathcal A = \f 4 3 \sqrt{\f 2 \pi} \f{\sigma_v^2 \sigma_w}{\sqrt{\sigma_v^2 + \sigma_a^2}}$ and $\mathcal B = \f 4 {9\pi}\f{ \sigma_v^4 }{\sigma_v^2 + \sigma_a^2}\lp \f 3 {\sigma_v^2 + \sigma_a^2} - \sigma_w^2\rp$ are as in \[thm:dalethExpSqrtTanhFullRes\] and $\mathcal B_b = \mathcal B + \f 1 2$ and $\mathcal B_w = \mathcal B - \f 1 2$.

$\alpha$-ReLU
-------------

[lemma]{}[VtPsiAlpha]{}\[lemma:VtPsiAlpha\] If $\alpha > -\f 1 2$, then ${\mathrm{V}}\psi_\alpha( q) = {\mathsf{c}}_\alpha q^{\alpha}$, where ${\mathsf{c}}_\alpha = \f 1 {\sqrt \pi} 2^{\alpha - 1}  \Gamma \left(\alpha+ \f 1 2\right)$.

Note that if $\alpha \le - \f 1 2$, then ${\mathrm{V}}\psi_\alpha( q)$ is not defined (its defining integral does not converge).

### Full Residual Network

By \[thm:fullResPQRec\] and \[lemma:VtPsiAlpha\], we have the length recurrences $$\begin{aligned}
{\mathbf{q}}&= \sigma_w^2 {\underline}{\mathbf{p}}+ \sigma_b^2\\
{\mathbf{p}}&= \sigma_v^2 {\mathsf{c}}_\alpha {\mathbf{q}}^\alpha + \sigma_a^2 + {\underline}{\mathbf{p}}\end{aligned}$$

[thm]{}[pDynamicAlphaReLU]{}\[thm:pDynamicAlphaReLU\] Suppose we have the nonlinearity $\phi = \psi_\alpha$. The in an FRN: If $\alpha = 1$, then $\p {\mathbf{p}}l = \Theta((1 + \sigma_v^2 \sigma_w^2/2)^l)$, with the hidden constant depending on the initial condition. If $0 < \alpha < 1$, then $\p {\mathbf{p}}l = \Theta(l^{\f 1 {1- \alpha}})$. More precisely, $\lim_{l \to \infty} {\mathbf{p}}/ l^{\f 1 {1-\alpha}} = [\sigma_v^2 \sigma_w^{2\alpha} {\mathsf{c}}_\alpha (1 - \alpha)]^{\f 1 {1-\alpha}}$.

\[fig:reluverifypasymptotics\] empirically verifies the asymptotics for $\alpha=1$ for various $\sigma_v$ and $\sigma_w$.

![Verification of the exponential asymptotics of $\p {\mathbf{p}}l$ when $\alpha=1$. The lines of each color correspond to different $(\sigma_w, \sigma_v)$ pairs, which are given in the legend. The solid lines are given by the recurrences \[thm:fullResPQRec\], and the dashed lines are given by our asymptotics $(1+\sigma_v^2\sigma_w^2/2)^l$ (\[thm:pDynamicAlphaReLU\]). Note that the y-axis is in log-scale.[]{data-label="fig:reluverifypasymptotics"}](graphics/relu_verify_p_asymptotics){height=".2\textheight"}

Similarly, by \[thm:full\_res\_l\_g\_recurr\], if ${\mathbf{q}}= {\mathbf{q}}'$, then $$\begin{aligned}
{\boldsymbol{\oldlambda}}&= \sigma_w^2 {\underline}{\boldsymbol{\oldgamma}}+ \sigma_b^2\\
{\boldsymbol{\oldgamma}}&= \sigma_v^2 {\mathbf{q}}^\alpha {\mathrm{W}}\psi_\alpha( 1, {\mathbf{c}}) + \sigma_a^2 + {\underline}{\boldsymbol{\oldgamma}}\end{aligned}$$

[thm]{}[ReLUSquaredConvergence]{}\[thm:ReLUSquaredConvergence\] Suppose $\phi = \psi_1$. Then in an FRN, $\p {\mathbf{e}}l \to 1$ and $1 - \p {\mathbf{e}}l \sim [\f 1 4 \sigma_v^2 \sigma_w^2 \inv B U l]^{-2}$ for $B = 1 + \sigma_v^2 \sigma_w^2/2$ and $U = \f {2\sqrt 2}{3\pi}$. As a result, $\p {\mathbf{s}}l = (1 - \p {\mathbf{e}}\l) \p {\mathbf{p}}l = \Theta(l^{-2}\exp(\Theta(l))) = \exp(\Theta(l)).$

[thm]{}[alphaReLUeConvergence]{} \[thm:alphaReLUeConvergence\] Suppose $\phi = \psi_\alpha$ for $0 < \alpha < 1$ in an FRN. Then ${\mathbf{e}}$ converges to the unique nonunit fixed point ${\mathbf{e}}^*$ of ${\mathbb{J}}_\alpha$, and $|{\mathbf{e}}^* - \p {\mathbf{e}}l|$ is ${{\check\Theta}}(l^{-\mu})$, where $\mu = (1-\dot{\mathbb{J}}_\alpha({\mathbf{e}}^*))/(1-\alpha)$. Additionally, $\p {\mathbf{s}}l = \Theta(\p {\mathbf{p}}l) = \Theta(l^{1/(1-\alpha)}).$

\[fig:6reluverifyestar\] verifies empirically that ${\mathbf{e}}^*$ is indeed the fixed point of $\p {\mathbf{e}}l$. \[fig:alphaReLUVerifyExponents\] verifies empirically the convergence rate $l^{-\mu}$. \[fig:MuPlots\] plots $\dot {\mathbb{J}}_\alpha({\mathbf{e}}^*)$ and $\mu$ versus $\alpha$. It certainly looks like $\mu = \f 1 2 (1 - \alpha)$, but we have no proof for it. Based on this conjecture, we see there is a “discontinuity” of $\mu$ at $\alpha = 1$: $\mu \to 0$ as $\alpha \to 1$, but for $\alpha = 1$, the actual convergence dynamics has exponent $-2$ by \[thm:ReLUSquaredConvergence\].

![Verification of fixed point ${\mathbf{e}}^*$ in \[thm:alphaReLUeConvergence\] for $\alpha = .6$. Different colors correspond to different initial conditions $\p {\mathbf{e}}0$, and the dashed line gives the fixed point.[]{data-label="fig:6reluverifyestar"}](graphics/6relu_verify_estar){height=".2\textheight"}

![[**(a)**]{} A plot of $\dot {\mathbb{J}}_\alpha({\mathbf{e}}^*)$ versus $\alpha$. [**(b)**]{} A plot of the exponent $\mu$ of the dynamics of $|{\mathbf{e}}^* - \p {\mathbf{e}}l|$ (see \[thm:alphaReLUeConvergence\]) []{data-label="fig:MuPlots"}](graphics/JdotPlot "fig:"){height=".2\textheight"} ![[**(a)**]{} A plot of $\dot {\mathbb{J}}_\alpha({\mathbf{e}}^*)$ versus $\alpha$. [**(b)**]{} A plot of the exponent $\mu$ of the dynamics of $|{\mathbf{e}}^* - \p {\mathbf{e}}l|$ (see \[thm:alphaReLUeConvergence\]) []{data-label="fig:MuPlots"}](graphics/MuPlot "fig:"){height=".2\textheight"}

Because of the following theorem, we cannot expect the equations of \[thm:dalethRecFull\] to hold for $\alpha \le \f 3 4$.

[thm]{}[dalethInfVarAlphaReLU]{} \[thm:dalethInfVarAlphaReLU\] Suppose we have the nonlinearity $\psi_\alpha$ in an FRN. $\Var(\dot \psi_\alpha(\zeta)^2)$ diverges for any Gaussian variable $\zeta$ with mean 0 if $\alpha \le \f 3 4$ but is finite if $\alpha > \f 3 4$.

[thm]{}[dalethDynamicsAlphaReLU]{} \[thm:dalethDynamicsAlphaReLU\] Suppose we have the nonlinearity $\psi_\alpha$ in an FRN. If $\alpha = 1$, then $\p{{\boldsymbol{\oldchi}}}{l-m} = \p {{\boldsymbol{\oldchi}}}{l} \left(\f 1 2 \sigma_v^2 \sigma_w^2 + 1\right)^m$. If $\alpha \in (\f 3 4, 1)$, then $\p{{\boldsymbol{\oldchi}}}{l-m} = \Theta(1) \p {{\boldsymbol{\oldchi}}}{l} (l/(l-m))^R$ for $R = \f{\alpha^2}{(1-\alpha)(2 \alpha - 1)}$, where the constants in $\Theta(1)$ do not depend on $l$ or $m$.

This exponent $\f{\alpha^2}{(1 - \alpha)(2\alpha - 1)}$ is minimized at $\alpha = \f 3 4$ on $\alpha \in (3/4, 1)$, where the value is $\f 9 2$ (and at $\alpha = \f 2 3$ on $\alpha \in (1/2, 1)$, where the value achieved is 4) (\[fig:backprop\_exponent\_alpha-relu\](a)).

As a corollary,

[thm]{}[alphaReLUAllGradients]{}\[thm:alphaReLUAllGradients\] If $\phi = \psi_1$ in an FRN, then for $l \ge m \ge 0,$ $$\begin{aligned}
\p {\boldsymbol{\oldchi}}{l-m} _b &= \Theta(1) \p {{\boldsymbol{\oldchi}}}l B^m,&
\p {\boldsymbol{\oldchi}}{l-m} _w &= \Theta(1) \p {{\boldsymbol{\oldchi}}}l B^l,\\
\p {\boldsymbol{\oldchi}}{l-m} _v&= \Theta(1) \p {{\boldsymbol{\oldchi}}}l B^l,&
\p {\boldsymbol{\oldchi}}{l-m} _a &= \Theta(1) \p {{\boldsymbol{\oldchi}}}{l} B^m.\end{aligned}$$ where $B = 1 + \sigma_v^2\sigma_w^2/2$.

If $\phi = \psi_\alpha$ in an FRN, for $\alpha < 1$, then for $l \ge m \ge 0,$ $$\begin{aligned}
\p {\boldsymbol{\oldchi}}{l-m} _b &= \Theta(1) \p {{\boldsymbol{\oldchi}}}l l^R (l-m)^{-R-1},&
\p {\boldsymbol{\oldchi}}{l-m} _w &= \Theta(1) \p {{\boldsymbol{\oldchi}}}l l^R (l-m)^{\f \alpha {1-\alpha} - R},\\
\p {\boldsymbol{\oldchi}}{l-m} _v &= \Theta(1) \p {{\boldsymbol{\oldchi}}}l l^R (l-m)^{\f \alpha {1-\alpha} - R},&
\p {\boldsymbol{\oldchi}}{l-m} _a &= \Theta(1) \p {{\boldsymbol{\oldchi}}}{l} (l/(l-m))^R.\end{aligned}$$

\[fig:alphaReLUVerifyExponents\] verifies the backward asymptotic dynamics empirically for different $\alpha < 1$. \[fig:backprop\_exponent\_alpha-relu\](b) graphs the exponent $\f \alpha {1-\alpha} - R$ in terms of $\alpha$. We see that on $[0.5, 1]$, the maximum of this exponent is at $\alpha = 1$.

![**(a)** The exponent of the polynomial gradient dynamics with respect to $\alpha$-ReLU versus $\alpha$. **(b)** The exponent of the dynamics of ${\boldsymbol{\oldchi}}_v$ and ${\boldsymbol{\oldchi}}_w$.[]{data-label="fig:backprop_exponent_alpha-relu"}](graphics/backprop_exponent_alpha-relu.pdf "fig:"){width=".4\textwidth"} ![**(a)** The exponent of the polynomial gradient dynamics with respect to $\alpha$-ReLU versus $\alpha$. **(b)** The exponent of the dynamics of ${\boldsymbol{\oldchi}}_v$ and ${\boldsymbol{\oldchi}}_w$.[]{data-label="fig:backprop_exponent_alpha-relu"}](graphics/alphaReLU_chi_w_chi_v_exponent.pdf "fig:"){width=".4\textwidth"}

Proofs
======

A brief note about notation: We use $\sim$ to denote both how a random variable is sampled (ex: $x \sim {\mathcal{N}}(0, 1)$ for a Gaussian $x$) and how a function behaves asymptotically, i.e. $f(x) \sim g(x)$ as $x \to a$ iff $\lim_{x \to a} f(x)/g(x) = 1$. Context should be enough to differentiate between these two cases. We in addition use $\simeq$ to denote asymptotic expansion. For example, if $\{\alpha_i\}_{i \ge 0}$ is a sequence of strictly decreasing reals and $\{\beta_i\}_{i \ge 0}$ is a sequence of nonzero reals, then $$f(x) \simeq \sum_{i \ge 0} \beta_i (x - \xi)^{\alpha_i}$$ means that as $x \to \xi$, $f(x) - \sum_{i = 0}^N \beta_i (x - \xi)^{\alpha_i} = \Theta((x- \xi)^{\alpha_{N+1}})$.

Preliminary Lemmas
------------------

\[lemma:cExpansion\] We have $$\f{\sigma_w^2  {\boldsymbol{\oldgamma}}+ \sigma_b^2}{\sigma_w^2  {\mathbf{p}}+ \sigma_b^2} =  {\mathbf{e}}(1 + O( {\boldsymbol{\oldgamma}}^{-1})).$$ regardless of whether $\p {\mathbf{e}}l = \p {\boldsymbol{\oldgamma}}l /\p {\mathbf{p}}l$ converges.

But suppose $\p {\mathbf{e}}l = \p {\boldsymbol{\oldgamma}}l /\p {\mathbf{p}}l \to {\mathbf{e}}^*$. If ${\mathbf{e}}^* < 1$, then $$\f{\sigma_w^2  {\boldsymbol{\oldgamma}}+ \sigma_b^2}{\sigma_w^2  {\mathbf{p}}+ \sigma_b^2} =  {\mathbf{e}}(1 + \Theta( {\boldsymbol{\oldgamma}}^{-1})).$$ If ${\mathbf{e}}^* = 1$, then $$\f{\sigma_w^2  {\boldsymbol{\oldgamma}}+ \sigma_b^2}{\sigma_w^2  {\mathbf{p}}+ \sigma_b^2} =  {\mathbf{e}}(1 + \Theta( {\epsilon}{\mathbf{p}}^{-1})),$$ where ${\epsilon}= 1 - {\mathbf{e}}$.

Write $M = \sigma_b^2/\sigma_w^2$. $$\begin{aligned}
\f{\sigma_w^2  {\boldsymbol{\oldgamma}}+ \sigma_b^2}{\sigma_w^2  {\mathbf{p}}+ \sigma_b^2} &=  {\mathbf{e}}(1 + \f{1 + M {\boldsymbol{\oldgamma}}^{-1}}{1 + M {\mathbf{p}}^{-1}})\\
&=  {\mathbf{e}}(1 + M(\inv {\boldsymbol{\oldgamma}}- \inv {\mathbf{p}}) + O(\inv {\mathbf{p}}(\inv {\boldsymbol{\oldgamma}}- \inv {\mathbf{p}}))).\end{aligned}$$ In any situation, $\inv {\boldsymbol{\oldgamma}}- \inv {\mathbf{p}}= O(\inv {\boldsymbol{\oldgamma}})$ because ${\boldsymbol{\oldgamma}}\le {\mathbf{p}}$, so this gives the first statement. If ${\mathbf{e}}^*$ exists and ${\mathbf{e}}^* < 1$, then $\inv {\boldsymbol{\oldgamma}}- \inv {\mathbf{p}}= \Theta(\inv {\boldsymbol{\oldgamma}})$, which yields the second statement. If ${\mathbf{e}}^*$ exists and ${\mathbf{e}}^* = 1$, then $\inv {\boldsymbol{\oldgamma}}- \inv {\mathbf{p}}= \inv {\mathbf{p}}((1-{\epsilon})^{-1} - 1) = \inv {\mathbf{p}}({\epsilon}+ O({\epsilon}^2)) = \Theta({\epsilon}\inv {\mathbf{p}})$.

For any function $f$ that is $(k+1)$-times differentiable in a neighborhood of $0$, we have the asymptotic expansion $$f(z) = \sum_{n = 0}^k \f{d^n f}{dz^n}(0) \f{z^n}{n!} + O(z^{k+1}), \text{as } z \to 0.$$ Since $$\begin{aligned}
    \left.\f{d^n}{d(1/q)^n}q^{1/2}{\mathrm{V}}\phi( q)\right\rvert_{q \to \infty} &= \f {(-1)^n} {2^n\sqrt{2\pi}}\int_{-\infty}^\infty \phi^2(z) z^{2n} \dd z\end{aligned}$$ whenever the RHS is integrable, we have

\[lemma:Vt\_asymptotic\] Suppose $\phi^2(z) z^{2n}$ is integrable over $z \in \R$ for all $0 \le n \le N + 1$. Then ${\mathrm{V}}\phi( q) = q^{-1/2} (\sum_{n = 0}^N C_n q^{-n} + O(q^{-N-1}))$ as $q \to \infty$, where $$C_n := \f {(-1)^n} {2^n n! \sqrt{2\pi}}\int_{-\infty}^\infty \phi^2(z) z^{2n} \dd z.$$

Note that ${\operatorname{sech}}^d(z) = \Theta(e^{-d|z|})$ for $z \to \infty$ as long as $d > 0$, so that $C_n$ from the above result converges when $\phi = {\operatorname{sech}}^d$. Therefore

\[lemma:Vt\_sech\_asymptotics\] Let $d > 0$. We have ${\mathrm{V}}{\operatorname{sech}}^d( q) \simeq q^{-1/2} \sum_{n \ge 0} C_n q^{-n}$, where $$C_n := \f {(-1)^n} {2^n n! \sqrt{2\pi}}\int_{-\infty}^\infty {\operatorname{sech}}^{2d}(z) z^{2n} \dd z.$$

As corollaries, we obtain the following asymptotics.

\[lemma:Vt\_dot\_tanh\_asymptotics\] ${\mathrm{V}}\dot \tanh( q) = \f 2 3 \sqrt{\f 2 \pi} q^{-1/2} + \Theta(q^{-3/2})$ as $q \to \infty$.

Use \[lemma:Vt\_sech\_asymptotics\] along with the fact that $\dot \tanh(z) = {\operatorname{sech}}^2(z)$ and $\int {\operatorname{sech}}^4 z \dd z = \f 2 3 \tanh z + \f 1 2 {\operatorname{sech}}^2 z \tanh z$.

\[lemma:vtanhSqrtConvergence\] $1 - {\mathrm{V}}\tanh( q) = \sqrt{\f 2 \pi} q^{-1/2} + \Theta(q^{-3/2})$ as $q \to \infty$.

Use \[lemma:Vt\_sech\_asymptotics\] along with the fact that $1 - \tanh^2(z) = {\operatorname{sech}}^2(z)$ and $\int {\operatorname{sech}}^2 z  \dd z = \tanh z$.

\[lemma:sech\_lower\_bound\] ${\operatorname{sech}}^2(t) \ge \exp(-t^2)$ for all $t$, with equality iff $t = 0$.

The lower bound is equivalent to $$\begin{aligned}
2 & \ge e^{t-t^2/2} + e^{-t -t^2/2}\end{aligned}$$ The RHS has derivative $(1 - t)e^{t - t^2/2} - (1 + t)e^{-t-t^2/2}$. This is 0 iff $$\begin{aligned}
\f{1 - t}{1 + t} = e^{-2t}\end{aligned}$$ which has a solution 0 and in general can only have solution $t\in (-1, 1)$ (by considering the sign of the LHS). Since each side is analytic in $t \in (-1, 1)$, we expand $$\begin{aligned}
\log \f{1 - t}{1 + t} &= \log e^{-2t}\\
\log(1-t) - \log(1+t) &= -2t\\
(-t-t^2-\cdots) - (t-t^2+\cdots) &= -2t\\
-2t-2t^3-\cdots &= -2t\end{aligned}$$ which shows that the only solution is $t = 0$. A simple plot shows that $t=0$ is a maximum, where the bound in question achieves equality.

\[lemma:V\_dot\_tanh\_lower\_bound\] Suppose $\phi = \tanh$. Then ${\mathrm{V}}\dot \phi( q) \ge \f 1 {\sqrt{4q+1}}$.

As a sanity check, \[lemma:Vt\_dot\_tanh\_asymptotics\] shows that ${\mathrm{V}}\dot \phi( q) \sim C_0 q^{1/2}$ where $C_0 \approx .5319$, which is above the .5 in this lemma.

By \[lemma:sech\_lower\_bound\], $$\begin{aligned}
    {\mathrm{V}}\dot\phi( q) &= \int \dd \mu(z) \dot \phi^2(\sqrt q z)\\
    &\ge \f 1 {\sqrt{2\pi}} \int \dd z \exp(-z^2/2 - 2qz^2)\\
    &= \f 1 {\sqrt{2\pi}} \int \dd z \exp(-(4q + 1)z^2/2)\\
    &= \f 1 {\sqrt{4q + 1}}.
    \end{aligned}$$

\[fig:V\_dot\_tanh\_lower\_bound\] demonstrates \[lemma:V\_dot\_tanh\_lower\_bound\].

![Illustration of \[lemma:V\_dot\_tanh\_lower\_bound\]: ${\mathrm{V}}\dot \phi( q)$ vs $\f 1 {\sqrt{4q+1}}$ for $\phi = \tanh$. This bound is very tight, and for most purposes, $\f 1 {\sqrt{4q+1}}$ can be taken as a good approximation of ${\mathrm{V}}\dot \phi( q)$.[]{data-label="fig:V_dot_tanh_lower_bound"}](graphics/V_dot_tanh_lower_bound.pdf){width=".4\textwidth"}

\[lemma:power\_sum\_asymptotics\] Let $d \in \R$ and $1 < M < N$ with $N - M \in \Z^{\ge 0}$. Set $\Sigma(M, N, d) := \sum_{a=M}^N a^d$. If we fix $M$ and let $N \to \infty$, $$\Sigma(M, N, d) = \begin{cases}
    \Theta(1) & \text{if $d < -1$}\\
    \log N + O(1) & \text{if $d = -1$}\\
    \f{N^{d+1}}{d+1} + O(1) & \text{if $-1 < d < 0$}\\
    N - M + 1 & \text{if $d =0$}\\
    \f 1 {d+1}N^{d+1} + \f 1 2 N^d + O(N^{\max(0, d-1)}) & \text{if $d > 0$}
\end{cases}$$

Consider the integrals $A = \int_M^{N+1} a^d \dd a$ and $B = \int^N_{M-1} a^d \dd a $. They evaluate to $A = \f 1 {d+1}((N+1)^{d+1} - M^{d+1})$ and $B = \f 1 {d+1}(N^{d+1} - (M-1)^{d+1})$ when $d \not = -1$ and to $A = \log(N+1) - \log M$ and $B = \log N - \log(M-1)$ when $d = -1$. When $d \le 0$, we have $A \le B$ and $\Sigma(M, N, d) \in [A, B]$; when $d > 0$, $B \le A$ and $\Sigma(M, N, d) \in [B, A].$ Thus, as $N \to \infty$ with $M$ fixed, when $d < -1$, $\Sigma(M, N, d) = \Theta(1)$; when $d = -1$, $\Sigma(M, N, -1) = \log N + O(1)$; and when $d > -1$, we have $\Sigma(M, N, d) = \f{N^{d+1}}{d+1} + O(N^{d})$.

Now for $a > 0$ and $d > -1$ and $d \not = 0, 1$, $$\begin{aligned}
    \int_a^{a+1} z^d - a^d \dd z &= \f 1 {d+1} ( (a+1)^{d+1} - a^d )\\
        &= (a^d + \f d 2 a^{d-1} + \cdots ) - a^d\\
        &= \f d 2 a^{d-1} + \Theta(a^{d-2}).
    \end{aligned}$$ where the hidden constants in $\Theta$ depend only on $d$ (and in fact this term vanishes if $d = 1$). Thus $$\begin{aligned}
    \Sigma(M, N, d) &= \int_M^{N+1} z^d \dd z - \sum_{a=M}^N [ \f d 2 a^{d-1} + \Theta(a^{d-2})]\\
        &= \f 1 {d+1} ((N+1)^{d+1} - M^{d+1}) - \f d 2 \Sigma(M, N, d-1) + \Theta(\Sigma(M, N, d-2))
    \end{aligned}$$ If $-1 < d < 0$, then $\Sigma(M, N, d - 1) = \Theta(1)$, so that $\Sigma(M, N, d) = \f{(N+1)^{d+1}}{d+1} + O(1) = \f{N^{d+1}}{d+1} + O(1).$ If $d > 0$ and $d \not = 1$, then $\Sigma(M, N, d - 1) = \f{N^d}d$, so that $$\begin{aligned}
    \Sigma(M, N, d) &= \f 1 {d+1} N^{d+1} + N^{d} + \Theta(N^{\max(0, d-1)}) - \f 1 2 N^d + \Theta(\Sigma(M, N, d-2)) \\
                &= \f 1 {d+1}N^{d+1} + \f 1 2 N^d + O(N^{\max(0, d-1)}).
    \end{aligned}$$

We can obtain more terms in the expansion for higher $d$ via the Euler-Maclaurin formula, but this suffices for our purposes.

Dynamics Zoo
------------

This section deduces the asymptotic behaviors of some sequences governed by recurrence equations. For the most part, the leading term of their asymptotic expansions is as one would expect from the corresponding differential equation. However, in some cases we need subleading terms for later results. They require slightly more nuanced reasoning. First we present a technical lemma.

\[lemma:simpleDynamics\] Let $F: \R \times \N \to \R$ be a function such that for a subset $U \sbe \R$, and for all $z, z' \in U, z \ge z' \implies F(z, n) \ge F(z', n)$ for every $n$. Suppose sequences $\p a l, \p b l, \p {c}l$ satisfy

-   $\p a {l+1} = F(\p a l, l)$ for all $l$;

-   $\p b {l+1} \le F(\p b l, l)$ for all $l$ above a constant $K_b$.

-   $\p {c}{l+1} \ge F(\p {c}l, l)$ for all $l$ above a constant $K_c$.

and furthermore, $\p a l, \p b l, \p {c}l$ all fall into $U$ for $l$ above a constant $K_U$. If for some $m \ge \max(K_b, K_U)$, $\p b m \le \p a m$, then $\p b l \le \p a l, \forall l \ge m$. Similarly, if for some $n \ge \max(K_c, K_U)$, $\p {c}n \ge \p a n$, then $\p {c}l \ge \p a l, \forall l \ge n$.

For the first claim: $\p b m \le \p a m \implies \p b {m+1} \le F(\p b m, m) \le F(\p a m, m) = \p a {m+1}$. Here the last inequality used the monotonicity of $F$. Induction gives the desired result.

It’s similar for the second claim, where the inductive step is $\p {c}m \ge \p a m \implies \p {c}{m+1} \ge F(\p {c}m, m) \ge F(\p a m, m) = \p a {m+1}$.

\[lemma:alphaDeltaRec\] Suppose $\p {\epsilon}l$ satisfies the recurrence $$\p {\epsilon}l = \p {\epsilon}{l-1} (1 + \f \delta {l^\beta}).$$ for some nonzero constant $\delta \in \R$ independent of $l$.

-   If $\beta > 1$, then $\p {\epsilon}l = \Theta(1)$.

-   If $\beta = 1$, then $\p {\epsilon}l = \Theta(l^{\delta})$.

-   If $0 < \beta < 1$, then $\p {\epsilon}l = \exp(\f{\delta}{1 - \beta}l^{1-\beta} + \tilde\Theta(l^{\psi_1(1-2\beta)}))$, where $\psi_1(x) = \max(0, x)$ is the ReLU function.

We have $$\begin{aligned}
    \log \p {\epsilon}l &= \log \p {\epsilon}{l-1} + \log(1 + \delta/l^\beta)\\
    &= \log \p {\epsilon}{l-1} + \delta/l^\beta + \Theta(\delta^2/l^{2\beta})
    \end{aligned}$$ for large $l$. If $\beta > 1$, then $\sum_l l^{-\beta}$ converges, and $$\begin{aligned}
    \log \p {\epsilon}l &= \log \p {\epsilon}0 - \Theta(1)\\
    \p {\epsilon}l &= \Theta(1).
    \end{aligned}$$ If $\beta = 1$, then $$\begin{aligned}
    \log \p {\epsilon}l &= \log \p {\epsilon}0 + \delta \log l + \Theta(1)\\
    \p {\epsilon}l &= \Theta(l^{\delta}).
    \end{aligned}$$ If $\beta < 1$, then $$\begin{aligned}
    \log \p {\epsilon}l &= \log \p {\epsilon}0 + \f\delta{1 - \beta} l^{1 - \beta} + \tilde\Theta(l^{1 - 2\beta})\\
    \p {\epsilon}l &= \exp( \f\delta{1 - \beta} l^{1 - \beta} + \tilde\Theta(l^{\psi_1(1 - 2\beta)})).
    \end{aligned}$$

\[lemma:alphaDeltaDynamics\] Suppose $\p {\epsilon}l = C l^{-\alpha} + \p {\epsilon}{l-1} (1 + \delta/l^{\beta})$ for $\alpha \in \R$, $C\not=0$, and $\delta \not = 0$. Then

-   If $\beta > 1$, then

    -   $\p {\epsilon}l = \Theta(l^{1 - \alpha})$ if $\alpha \in (0, 1)$;

    -   $\p {\epsilon}l = \Theta(\log l)$ if $\alpha = 1$;

    -   $\p {\epsilon}l = \Theta(1)$ if $\alpha > 1$.

-   If $\beta = 1$, then

    -   $\p {\epsilon}l = \Theta(l^{\max(\delta, 1-\alpha)})$ if $1-\delta \not = \alpha$.

    -   $\p {\epsilon}l = \Theta(l^{\delta} \log l)$ if $1-\delta = \alpha$.

Furthermore, for $\beta = -\delta = 1$, $\p {\epsilon}l \sim l^{-1}$ if $\alpha > 2$, $\p {\epsilon}l \sim l^{1-\alpha}$ if $\alpha < 2$, and $\p {\epsilon}l \sim l^{\delta} \log l$ if $\alpha = 2$.

We can unwind the recurrence to get $$\begin{aligned}
    \p {\epsilon}l &= \sum_{m=1}^l m^{-\alpha} \prod_{n=m+1}^l (1 + \f \delta {n^\beta}) + \p {\epsilon}0 \prod_{n=1}^l (1 + \f \delta {n^\beta})
    \end{aligned}$$

Suppose $\beta > 1$. By \[lemma:alphaDeltaRec\], we get $$\begin{aligned}
    \p {\epsilon}l &= \Theta(1)\sum_{m=1}^l m^{-\alpha} + \p {\epsilon}0 \Theta(1)\\
    &= \begin{cases}
    \Theta(l^{1 - \alpha}) & \text{if $\alpha \in (0, 1)$}\\
    \Theta(\log l) & \text{if $\alpha = 1$}\\
    \Theta(1) & \text{if $\alpha > 1$.}
    \end{cases}
    \end{aligned}$$

Now suppose $\beta = 1$. By \[lemma:alphaDeltaRec\], we get $$\begin{aligned}
    \p {\epsilon}l &= \sum_{m=1}^l m^{-\alpha} \Theta(m^{-\delta}l^\delta) + \p {\epsilon}0 \Theta(l^{\delta})
    \end{aligned}$$ where the constants hidden inside the $\Theta$ are the same in every term of the sum. If $\alpha > 1 - \delta$, then $m^{-\delta - \alpha} = o(m^{-1})$, so that $\sum_{m=1}^l m^{-\delta - \alpha} = \Theta(1)$, and $$\begin{aligned}
    \p {\epsilon}l &= \Theta(l^{\delta}) + \p {\epsilon}0 \Theta(l^{\delta})\\
    &= \Theta(l^{\delta}).
    \end{aligned}$$ On the other hand, if $\alpha < 1 - \delta$, then $\sum_{m=1}^l m^{-\delta - \alpha} = \Theta(l^{1 - \delta - \alpha})$. So $$\begin{aligned}
    \p {\epsilon}l &= \Theta(l^{1 - \alpha}) + \p {\epsilon}0 \Theta(l^{\delta})\\
    &= \Theta(l^{1 - \alpha}).
    \end{aligned}$$ If $\alpha = 1 - \delta$, then $\sum_{m=1}^l m^{-\delta - \alpha} = \Theta(\log l)$. So $$\begin{aligned}
    \p {\epsilon}l &= \Theta(l^{\delta}\log l) + \p {\epsilon}0 \Theta(l^{\delta})\\
    &= \Theta(l^{\delta}\log l).
    \end{aligned}$$

Finally, if $\beta \in (0, 1)$, then $$\begin{aligned}
    \p {\epsilon}l &= e^{\f\delta{1-\beta}l^{1-\beta} + \Theta(l^{1-2\beta})}\sum_{m=1}^l m^{-\alpha} e^{\f{-\delta}{1-\beta} m^{1-\beta} + \Theta(m^{1-2\beta})} + e^{\f\delta{1-\beta}l^{1-\beta} + \Theta(l^{1-2\beta})}
    \end{aligned}$$ The case of $\delta = -1$ telescopes, so that the upper and lower constants hidden in $\Theta$ can both be taken to be 1.

\[lemma:invlogDynamics\] Suppose for some $\beta > 0$, a sequence $\p {\epsilon}l$ satisfies $$\p {\epsilon}l = \p {\epsilon}{l-1} (1 - \mu (\p {\epsilon}{l-1})^\beta/l),\quad \p {\epsilon}0 \in (0, \f 1 \mu).$$ Then $\p {\epsilon}l \sim (\beta\mu \log l)^{-1/\beta}$.

Consider the differential equation $$\dot x_\mu = - \mu x^{\beta+1}_\mu/t$$ for constant $\mu$ has solution $x_\mu = [\beta(\mu \log t + C)]^{-1/\beta}$ for some constant $C$ determined by initial condition. Note that $$-\mu x_\mu(t)^{\beta+1}/t \le x_\mu(t+1) - x_\mu(t) \le - \mu x_\mu(t+1)^{\beta+1}/(t+1) = - (1 - o(t^{-1}))\mu x_\mu (t)^{\beta+1}/t.$$ For any small enough $\alpha > 0$, we apply \[lemma:simpleDynamics\] with $F({\epsilon}, l) = {\epsilon}- \mu {\epsilon}^{\beta+1}/l$ (which is monotonic in ${\epsilon}$ for small enough ${\epsilon}$), $\p {c}l = x_\mu(l)$, and $\p b l = x_{\mu-\alpha}(l)$ to obtain $$x_{\mu - \alpha}(l) \le \p {\epsilon}l \le x_\mu(l)$$ for large enough $l$ and appropriately chosen initial conditions. This shows that $\p {\epsilon}l = \Theta(\log l ^{-1/\beta})$ Taking $\alpha \to 0$, we also obtain the leading coefficient $\p {\epsilon}l \sim [\beta\mu \log l]^{-1/\beta}$.

\[lemma:polyDynamics\] Suppose a sequence $\p u l$ is governed by the equation $$\p u l - \p u {l-1} = A(\p u {l-1} + B)^\alpha,$$ where $\alpha \in [0, 1)$ and $A > 0$. Then $\p u l = K_1 l^{{\f 1 {1-\alpha}}}- K_2 l^{{\f \alpha {1-\alpha}}}\log l + o(l^{{\f \alpha {1-\alpha}}}\log l)$, where $K_1 = [A(1-\alpha)]^{{\f 1 {1-\alpha}}}$ and $K_2 = \f 1 2 A^{{{\f 1 {1-\alpha}}}} (1-\alpha)^{{{\f \alpha {1-\alpha}}}- 1} \alpha$.

**Leading term.** The differential equation $$\dot x_{A, B} = A(x_{A, B} + B)^\alpha$$ has solution $x_{A,B}(l) = [A(1-\alpha)(l + S)]^{\f 1 {1-\alpha}} - B$ for some constant $S$. Since $\dot x_{A, B}$ is monotonic, we have (writing $x = x_{A, B}$ for brevity) $$A(x_{A, B}(l) + B)^\alpha = \dot x_{A, B}(l) \le x_{A, B}(l+1) - x_{A, B}(l) \le \dot x_{A, B}(l+1) \le (A + o(1))(x_{A, B}(l) + B)^\alpha$$ for large enough $l$. We apply \[lemma:simpleDynamics\] with $F(x, l) = x + A(x + B)^\alpha$ (which is monotonic in $x$ for large $x$), $\p {c}l = x_{A, B}(l)$, and $\p b l = x_{A-{\epsilon}, B}(l)$ to obtain $$x_{A-{\epsilon},B}(l) \le \p u l \le x_{A, B}(l)$$ for large enough $l$ and appropriate initial conditions. Therefore $\lim \p u l / l^{\f 1 {1-\alpha}} \in [[(A - {\epsilon}) ( 1- \alpha)]^{\f 1 {1 - \alpha}}, [A ( 1- \alpha)]^{\f 1 {1 - \alpha}}].$ Taking ${\epsilon}\to 0$ gives the leading term.

**Subleading term.** Now let $\p v l := \p u l - \aleph l^{\f 1 {1-a}}$, where $\aleph = [A(1-\alpha)]^{\f 1 {1-\alpha}}.$ Then we have the recurrence $$\begin{aligned}
    \p v {l+1} + \aleph (l + 1)^{{{\f 1 {1-\alpha}}}} - \p v l - \aleph l^{{{\f 1 {1-\alpha}}}} &= A(\p v l + \aleph l^{{{\f 1 {1-\alpha}}}} + B)^\alpha\\
    \p v {l+1} - \p v l + \aleph({{\f 1 {1-\alpha}}}l^{{{\f \alpha {1-\alpha}}}} &+ \f 1 2 ({{\f 1 {1-\alpha}}}) ({{\f \alpha {1-\alpha}}}) l^{{{\f \alpha {1-\alpha}}}-1} + \Theta(l^{{{\f \alpha {1-\alpha}}}- 2}))  \\&= A[\aleph^\alpha l^{{{\f \alpha {1-\alpha}}}} + \alpha (\p v l+B) \aleph^{\alpha - 1} \inv l + \Theta((\p v l+B) l^{-1-{{\f 1 {1-\alpha}}}})]\\
    \p v {l+1} - \p v l &= {{\f \alpha {1-\alpha}}}\p v l \inv l - \f 1 2 \aleph({{\f 1 {1-\alpha}}})({{\f \alpha {1-\alpha}}}) l^{{{\f \alpha {1-\alpha}}}- 1} + g(l)
    \end{aligned}$$ for some $g(l) = O(l^{{{\f \alpha {1-\alpha}}}- 2} + \inv l)$ and where, to get the last equation, we have used $A \alpha^\alpha = {{\f 1 {1-\alpha}}}\aleph$ to cancel the $l^{{{\f \alpha {1-\alpha}}}}$ term and simplified $\alpha A \aleph^{\alpha-1} = {{\f \alpha {1-\alpha}}}$.

For any $J > 0$, the differential equation $\dot v_J(l) = {{\f \alpha {1-\alpha}}}v_J(l) \inv l - J l^{{{\f \alpha {1-\alpha}}}- 1}$ has solution $v_J(l) = C[l(1-\alpha)]^{{{\f \alpha {1-\alpha}}}} - J l^{{{\f \alpha {1-\alpha}}}} \log l$. Note that the functions $F_J(z, n) = z + {{\f \alpha {1-\alpha}}}z \inv n - J n^{{{\f \alpha {1-\alpha}}}- 1}$ and $G_J(z, n) = F_J(z, n) + g(n)$ is monotonic in $z$ (for positive $n$). For large $l$, we also have $\dot v_J(l)$ and $F_J(v_J(l), l) = v_J(l) + \dot v_J(l)$ decreasing in $l$. Thus for any ${\epsilon}> 0$ and $l$ large enough $$G_{J+{\epsilon}}(v_J(l), l) \le F_{J+{\epsilon}/2}(v_J(l), l) \le v_J(l) + \dot v_J(l+1) \le v_J(l+1) \le F_J(v_J(l), l) \le G_{J-{\epsilon}}(v_J(l), l).$$ Now apply \[lemma:simpleDynamics\] with $F = G_{K}$, $\p a l = \p v l, \p {c}l = v_{K-{\epsilon}}, \p b l = v_{K+{\epsilon}}$ where $K := \f 1 2 \aleph({{\f 1 {1-\alpha}}})({{\f \alpha {1-\alpha}}}) = \f 1 2 A^{{{\f 1 {1-\alpha}}}} (1-\alpha)^{{{\f \alpha {1-\alpha}}}- 1} \alpha$, with appropriately chosen initial conditions. This yields $\lim_{l\to \infty} \p v l / (l^{{{\f \alpha {1-\alpha}}}} \log l) \in [-K - {\epsilon}, -K + {\epsilon}]$ for every ${\epsilon}> 0$, and there it must be equal to $K$. We have thus obtained the asymptotic expansion $$\p u l = [A(1-\alpha) l]^{{\f 1 {1-\alpha}}}- \f 1 2 A^{{{\f 1 {1-\alpha}}}} (1-\alpha)^{{{\f \alpha {1-\alpha}}}- 1} \alpha l^{{\f \alpha {1-\alpha}}}\log l + o(l^{{\f \alpha {1-\alpha}}}\log l).$$

\[lemma:polyDynamicsPosAlpha\] Suppose a sequence $\p u l$ is governed by the equation $$\p u l - \p u {l-1} = -A(\p u {l-1} + B)^\alpha,$$ where $\alpha > 1$ and $A > 0$. Then $\p u l \sim [A(\alpha-1)l]^{\f 1 {1-\alpha}}$.

Similar to \[lemma:polyDynamics\].

\[lemma:polyDynamicsConstant\] Suppose a sequence $\p u l$ is governed by the equation $$\p u l - \p u {l-1} = A(\p u {l-1} + B)^\alpha + C,$$ where $\alpha \in (0, 1)$. Then $\p u l = K_1 l^{{\f 1 {1-\alpha}}}+ R(l)$, where the remainder $R(l)$ is $$R(l) \sim \begin{cases}
    -K_2 l^{{{\f \alpha {1-\alpha}}}} \log l &\text{if $\alpha > \f 1 2$}\\
    (C-K_2) l \log l            &\text{if $\alpha = \f 1 2$ and $K_2 \not=C$}\\
    \f{C(1-\alpha)}{1-2\alpha} l    & \text{if $ \alpha < \f 1 2$}
    \end{cases}$$ where $K_1 = [A(1-\alpha)]^{\f 1 {1-\alpha}}, K_2 = \f 1 2 A^{{{\f 1 {1-\alpha}}}} (1-\alpha)^{{{\f \alpha {1-\alpha}}}- 1} \alpha$ as in \[lemma:polyDynamics\].

$u$ is bounded below by the dynamics $\p v l - \p v {l-1} = A(\p v {l-1} + B)^\alpha$ and bounded above by the dynamics $\p w l - \p w {l-1} = (A + o(1))(\p w {l-1} + B)^\alpha$. By \[lemma:polyDynamics\], both $v$ and $w$ are asymptotic to $\p u l \sim [A(1-\alpha)l]^{\f 1 {1-\alpha}}$, which gives the result.

Now define $\p v l = \p u l - [A(1-\alpha)l]^{{\f 1 {1-\alpha}}}$, and similar to the proof of \[lemma:polyDynamics\], we find $$\p v {l+1} - \p v l = {{\f \alpha {1-\alpha}}}\p v l \inv l - K l^{{{\f \alpha {1-\alpha}}}- 1} + C + g(l)$$ where $K = \f 1 2 A^{{{\f 1 {1-\alpha}}}} (1-\alpha)^{{{\f \alpha {1-\alpha}}}- 1} \alpha$ and $g(l) = O(l^{{{\f \alpha {1-\alpha}}}- 2} + \inv l)$. If ${{\f \alpha {1-\alpha}}}> 1 \iff \alpha > \f 1 2$, then $C + g(l) = o(l^{{{\f \alpha {1-\alpha}}}- 1})$ and we can proceed as in the proof of \[lemma:polyDynamics\] to find $\p v l \sim K l^{{{\f \alpha {1-\alpha}}}} \log l.$ If ${{\f \alpha {1-\alpha}}}= 1 \iff \alpha = 1$ and $K \not = C$, then $\p v {l+1} - \p v l = {{\f \alpha {1-\alpha}}}\p v l l^{-1} - (K-C) l^{{{\f \alpha {1-\alpha}}}- 1} + g(l)$, so that the technique used in \[lemma:polyDynamics\] would obtain $\p v l \sim (K-C) l^{{{\f \alpha {1-\alpha}}}} \log l = (K-C) l \log l.$ If ${{\f \alpha {1-\alpha}}}< 1 \iff \alpha < \f 1 2$, then $\p v {l+1} - \p v l = {{\f \alpha {1-\alpha}}}\p v l l^{-1} + C + o(1)$, then by using the differential equation $\dot v_J(l) = {{\f \alpha {1-\alpha}}}v_J(l) l^{-1} + J$ to approximate the difference equation solution and applying \[lemma:simpleDynamics\] as in the proof of \[lemma:polyDynamics\], we obtain $\p v l(l) \sim \f{C(1-\alpha)}{1-2\alpha} l$.

Forward Dynamical Equations
---------------------------

Here we derive the recurrences governing the forward length and correlation quantities ${\mathbf{p}}, {\mathbf{q}}, {\boldsymbol{\oldlambda}}, {\boldsymbol{\oldgamma}}.$ We start with reduced residual networks.

We have $$\begin{aligned}
    {\mathbf{q}}&= \la h_j^2 \ra = \la \sum_i (w_{ji}{\underline}x_i + b_j)^2 \ra\\
    &= \la b_j^2 \ra + \sum_i \la w_{ji}^2 {\underline}x_i^2\ra + 2\sum_i \la w_{ji} {\underline}x_i b_j \ra + 2 \sum_{j \not= l} \la w_{ji} w_{li} x_i^2\ra 
    \end{aligned}$$ But $w_{ji}, w_{li}, {\underline}x,$ and $b_j$ form an independency, so the last two sums are 0, and the terms in the first sum split multiplicatively. Therefore $$\begin{aligned}
    {\mathbf{q}}&= \sigma_b^2 + \sum_i \la w_{ji}^2 \ra \la {\underline}x_i^2 \ra\\
    &= \sigma_b^2 + N \cdot \f {\sigma_w^2}{N} {\underline}{\mathbf{p}}\\
    &= \sigma_b^2 + \sigma_w^2 {\underline}{\mathbf{p}}. 
    \end{aligned}$$

For the recurrence of ${\mathbf{p}}$, we have $$\begin{aligned}
    {\mathbf{p}}&= \la x_i^2 \ra = \la (\phi(h_i) + {\underline}x_i)^2 \ra \\
    &= \la \phi(h_i)^2\ra + \la {\underline}x_i^2 \ra + 2 \la\phi(h_i) {\underline}x_i \ra
    \end{aligned}$$ As $N \to \infty$, the coefficient $w_{ii}$ of ${\underline}x_i$ in $h_i$ has vanishing covariance, so $h_i$ and ${\underline}x_i$ become independent. Therefore $\la \phi(h_i) {\underline}x_i \ra = \la \phi(h_i) \ra \la {\underline}x_i \ra$. Because $h_i$ is the sum of a large number of independent random variables, by CLT, $h_i$ is a Gaussian with mean $\sum_i \la w_{ji}\ra \la {\underline}x_i\ra + \la b_j\ra = 0$ since $\la w_{ji} \ra = \la b_j \ra = 0$. Our antisymmetry assumption on $\phi$ then implies $\la \phi(h_i) \ra = 0$. Therefore, $$\begin{aligned}
    {\mathbf{p}}&= \la \phi(h_i)^2\ra + \la {\underline}x_i^2 \ra\\
    &= {\mathrm{V}}\phi( {\mathbf{q}}) + {\underline}{\mathbf{p}}\end{aligned}$$ as desired.

Similar to \[lemma:p\_q\_recurrence\].

Now, for the full residual networks, the proofs are similar, but we no longer need to assume that $\phi$ is antisymmetric because of the randomization via the extra sets of weights.

$$\begin{aligned}
    {\mathbf{q}}&= \la  h_j^2 \ra
    = \la ( w_j^i {\underline}x_i +  b_j)^2\ra
    = \la ( w_j^i {\underline}x_i)^2 \ra  + \la  b_j^2\ra\\
    &= \sigma^2_w \la {\underline}x_i^2 \ra + \sigma_b^2\\
    &= \sigma^2_w {\underline}{\mathbf{p}}+ \sigma^2_b\\
    {\mathbf{p}}&= \la  x_i^2 \ra
    = \la (v^j_i\phi( h_j) +{\underline}x_i + a_i)^2 \ra\\ 
    &= \sigma_v^2\la \phi( h_i)^2 \ra + \la {\underline}x_i^2 \ra + \sigma_a^2\\
    &=\sigma_v^2 {\mathrm{V}}\phi(  {\mathbf{q}}) + \sigma_a^2 + {\underline}{\mathbf{p}}\\
    \end{aligned}$$

where in the third equality for ${\mathbf{p}}$, we are now using the independence of $v_i^j$ from all other variables to cancel out the terms, whereas before we had to rely on $\phi$ being antisymmetric.

Similar to \[thm:fullResPQRec\].

Backward Dynamical Equations
----------------------------

Here we derive the recurrences governing the gradient quantities ${{\boldsymbol{\oldchi}}}$ and ${\boldsymbol{\oldchi}}_{\bullet}$ for different $\bullet$, all under the gradient independence assumption. Write $\p {\beta}l_i = \pdf E {\p x l_i}$ for a cost function $E$.

For a reduced residual network, we have the following derivative computation: $$\begin{aligned}
    \pdf {x_i} {{\underline}x_j} &= \delta_{ji} + \dot \phi(h_i) \pdf {h_i}{{\underline}x_j},&
    \pdf {x_i}{h_j} &= \delta_{ji} \dot \phi(h_j),&
    \pdf {h_i}{{\underline}x_j} &= w_{ij},&
    \pdf {h_i}{w_{ij}} &= {\underline}x_j,&
    \pdf {h_i}{b_j} &= \delta_{ij}.
    \end{aligned}$$

Then $$\begin{aligned}
        {\underline}{\beta}_j &= {\beta}_j + \sum_i {\beta}_i \dot \phi(h_i) \pdf {h_i}{{\underline}x_j}\\
    &= {\beta}_j + \sum_i {\beta}_i\dot \phi(h_i) w_{ij}\\
    \la {\underline}{\beta}_j^2 \ra &= \la [{\beta}_j + \sum_i {\beta}_i\dot \phi(h_i) w_{ij}]^2 \ra\\
    &= \la {\beta}_j^2 \ra + \sum_i \la {\beta}_i^2 \dot \phi^2(h_i) (w_{ij})^2\ra \\
    &\phantom{={}}+ 2\sum_{i < k} \la {\beta}_i {\beta}_k\dot\phi(h_i)w_{ij}\dot\phi(h_k)w_{kj}\ra + 2 \sum_i\la {\beta}_j {\beta}_i \dot\phi(h_i)w_{ij} \ra
        \end{aligned}$$ The last two terms of the above vanish as $w_{ij}$ is independent from $w_{kj}$, $h_i, h_k$ and ${\beta}_i, {\beta}_j, {\beta}_k$ by \[ass:gradInd\], and $\la w_{ij} \ra = 0$.

Therefore, applying \[ass:symAct\], $$\begin{aligned}
    \la {\underline}{\beta}_j^2 \ra &=  \sigma_w^2\la {\beta}_j^2\ra\la\dot \phi^2(h_i)\ra + \la {\beta}_j^2\ra\\
    &= (\sigma_w^2  {\mathrm{V}}\dot \phi( {\mathbf{q}}) + 1)\la {\beta}_j^2 \ra
    \end{aligned}$$

We similarly have $$\begin{aligned}
    \pdf E {b_j} &= \sum_i \pdf E {x_i} \pdf {x_i}{h_j}
    = {\beta}_j \dot \phi(h_j), & \text{since $\pdf {x_i}{h_j} = \delta_{ji} \dot \phi(h_j)$}\\
    \la \lp\pdf E {b_j}\rp^2 \ra &= \la {\beta}_j^2 \dot \phi(h_j)^2 \ra
    = \la {\beta}_j^2 \ra {\mathrm{V}}\dot \phi( {\mathbf{q}}), & \text{by \cref{ass:gradInd}(b)};\\
    \pdf E {w_{ji}} &= \sum_i \pdf E {x_i} \pdf {x_i}{h_j}\pdf {h_j}{w_{ji}}
    = {\beta}_j \dot \phi(h_j){\underline}x_i, & \text{since $\pdf {x_i}{h_j} = \delta_{ji} \dot \phi(h_j)$}\\
    \la \lp\pdf E {w_{ji}}\rp^2 \ra &= \la {\beta}_j^2 \dot \phi^2(h_j) {\underline}x_i^2 \ra
    = \la {\beta}_j^2 \ra {\mathrm{V}}\dot \phi( {\mathbf{q}}) {\underline}{\mathbf{p}}, & \text{by \cref{ass:gradInd}(b)}
    \end{aligned}$$ In the last equation we have also used the fact that as $N \to \infty$, $h_j$ and $x_i$ become independent (they are jointly Gaussian and their correlation $\la w_{ji}^2 \ra$ goes to 0 with $N$).

For the full residual network, we have the following derivative computations: $$\begin{aligned}
    \pdf {x_i} {{\underline}x_j} &= \delta_{ji} + \sum_k v_{ik}\dot \phi(h_k) \pdf {h_k}{{\underline}x_j},&
    \pdf {x_i}{h_j} &= v_{ij} \dot \phi(h_j),&
    \pdf {h_i}{{\underline}x_j} &= w_{ij},&
    \pdf {h_i}{w_{ij}} &= {\underline}x_j,&
    \pdf {h_i}{b_i} &= 1,\\
    &&
    \pdf {x_i}{v_{ik}} &= \phi(h_k),&
    \pdf {x_i}{a_i} &= 1.
    \end{aligned}$$

Again let ${\beta}_j = \pdf E {x_j}$. Then $$\begin{aligned}
    {\underline}{\beta}_j &= \sum_i {\beta}_i (\delta_{ji} + \sum_k v_{ik}\dot \phi(h_k) \pdf {h_k}{{\underline}x_j})\\
    &= \sum_i {\beta}_i (\delta_{ji} + \sum_k v_{ik}\dot \phi(h_k) w_{kj})
    \end{aligned}$$

Thus, $$\begin{aligned}
    \la {\underline}{\beta}_j^2 \ra &= \la [\sum_i {\beta}_i (\delta_{ji} + \sum_k v_{ik}\dot \phi(h_k) w_{kj})]^2 \ra\\
    &= \la {\beta}_j^2 \ra + \sum_{i, k} \la v_{ik}^2 \ra \la w_{kj}^2\ra {\mathrm{V}}\dot \phi( {\mathbf{q}}) \la {\beta}_i^2 \ra\\
    &= \la {\beta}_j^2 \ra (1 + \sigma_v^2 \sigma_w^2 {\mathrm{V}}\dot \phi( {\mathbf{q}}))
        \end{aligned}$$ where in the second equality we applied the independence argument as in the proof of \[thm:dalethRecReduced\], leveraging \[ass:gradInd\], and in the third equality we used \[ass:symAct\] to get $\la {\beta}_i^2\ra = \la {\beta}_j^2 \ra$.

The other computations are similar to the proof of \[thm:dalethRecFull\].

Tanh: Reduced Residual Network
------------------------------

### Forward Dynamics

The case with $\sigma_w = 0$ is trivial. We assume $\sigma_w > 0$ from here on.

**${\mathbf{p}}$ and ${\mathbf{q}}$ are asymptotically linear with $l$.** We first show that, for any $\omega < 1$, $$l + \p {\mathbf{p}}0 \ge \p {\mathbf{p}}l \ge \omega l$$ and $$\sigma_w^2 (l + \p {\mathbf{p}}0) + \sigma_b^2 \ge \p {\mathbf{q}}l \ge \sigma_w^2 \omega (l-1) + \sigma_b^2,$$ so that $\p {\mathbf{p}}l \sim l$ and $\p {\mathbf{q}}l \sim \sigma_w^2 l$.

The upper bounds are trivial, given ${\mathrm{V}}\phi( {\mathbf{q}}) \le 1$ for any ${\mathbf{q}}$. We show the lower bounds for any $\omega < 1$.

For any ${\epsilon}> 0$, define $\aleph_{\epsilon}$ by $\phi^2(\aleph_{\epsilon}) = \exp(-{\epsilon})$. Then $$\begin{aligned}
{\mathrm{V}}\phi( {\mathbf{q}}) &\ge \exp(-{\epsilon}) \Pr[z \not \in [-\aleph_{\epsilon}, \aleph_{\epsilon}]: z \sim {\mathcal{N}}(0, {\mathbf{q}})]\\
    &\ge \exp(-{\epsilon})\lp 1 - \f{2 \aleph_{\epsilon}}{\sqrt{2\pi {\mathbf{q}}}} \rp\end{aligned}$$ where the second inequality follows from an overestimate of the $\Pr[z \in [-\aleph_{\epsilon}, \aleph_{\epsilon}]]$ via the mode of ${\mathcal{N}}(0, {\mathbf{q}})$.

For any ${\mathbf{q}}\ge \p {\mathbf{q}}0$, ${\mathrm{V}}\phi( {\mathbf{q}})$ is then lower bounded by $$\phi^2\lp \sqrt{\p {\mathbf{q}}0} \rp
\lp1 - \f{2\sqrt{\p {\mathbf{q}}0}}{\sqrt{2 \pi \p {\mathbf{q}}0}}\rp
= \phi^2\lp\sqrt{\p {\mathbf{q}}0}\rp
\lp1 - \sqrt{\f 2 \pi}\rp > 0.$$ Thus $\p {\mathbf{p}}l$ and $\p {\mathbf{q}}l$ are unbounded with $l$.

Furthermore, as ${\mathbf{q}}\to \infty$, the lower bound $\exp(-{\epsilon})\lp 1 - \f{2 \aleph_{\epsilon}}{\sqrt{2\pi {\mathbf{q}}}} \rp$ goes to $\exp(-{\epsilon})$, for any ${\epsilon}$. Therefore, for any $\omega < 1, \p {\mathbf{p}}l \ge \omega l$ and $\p {\mathbf{q}}l \ge \sigma_w^2 \omega (l-1) + \sigma_b^2$.

**Asymptotic expansion.** Now we repeat the following to get each successive asymptotic term of $\p {\mathbf{p}}l$ and $\p {\mathbf{q}}l$: We plug in the current asymptotic form of $\p {\mathbf{q}}l$ into ${\mathrm{V}}\tanh( {\mathbf{q}}) = 1 - C {\mathbf{q}}^{-1/2} + \Theta({\mathbf{q}}^{-3/2})$ (\[lemma:vtanhSqrtConvergence\]), where $C = \sqrt{2/\pi}$. Next we take the sum $\p {\mathbf{q}}l = \sum_{r=1}^l {\mathrm{V}}\tanh( \p {\mathbf{q}}r)$, which yields one more term in the asymptotic expansion of ${\mathbf{p}}$ than the last round. We then repeat until we get only constant terms.

The following exhibits a trace of this procedure, where in the summation step for $\p {\mathbf{q}}l$, we implicitly apply $$\begin{aligned}
{\mathbf{q}}&= \sigma_w^2 l + o(l) = \sigma_w^2 l (1 + o(1))\\
{\mathbf{q}}^{-1/2} &= \sigma_w^{-1} l^{-1/2}(1 + o(1)) = \sigma_w^{-1} l^{-1/2} + o(l^{-1/2})\\
{\mathbf{p}}&= \sum_{r=1}^l 1 - C (\p {\mathbf{q}}r)^{-1/2} + \Theta((\p {\mathbf{q}}r)^{-3/2}) \\
    &= \sum_{r=1}^l 1 - C (\sigma_w^{-1} r^{-1/2} + o(r^{-1/2})) + \Theta(r^{-3/2}) \\
    &= l - 2C \sigma_w^{-1} l^{1/2} + o(l^{1/2})\\
{\mathbf{q}}&= \sigma_w^2 l - 2 C \sigma_w l^{1/2} + o(l^{1/2}) = \sigma_w^2 l (1 - 2 C \sigma_2^{-1} l^{-1/2} + o(l^{-1/2}))\\
{\mathbf{q}}^{-1/2} &= \sigma_w^{-1} l^{-1/2} ( 1 + C \sigma_w^{-1} l^{-1/2} + o(l^{-1/2}))= \sigma_w^{-1} l^{-1/2} + C \sigma_w^{-2} l^{-1} + o(l^{-1})\\
{\mathbf{p}}&= \sum_{r=1}^l 1 - C(\sigma_w^{-1} l^{-1/2} + C \sigma_w^{-2} l^{-1} + o(l^{-1})) + \Theta(l^{-3/2})\\  
    &= l - 2 C \sigma_w^{-1} l^{1/2} - C^2 \sigma_w^{-2} \log l + o(\log l)\\
{\mathbf{q}}&= \sigma_w^2 l (1 - 2 C \sigma_w^{-1} l^{-1/2} - C^2 \sigma_w^{-2} \f{\log l}{l} + o(\f{\log l}{l}))\\
{\mathbf{q}}^{-1/2} &=\sigma_w^{-1} l^{-1/2}(1 + C \sigma_w^{-1} l ^{-1/2} + \f 1 2 C^2 \sigma_w^{-2} \f{\log l}{l} + o (\f{\log l}{l}))\\
{\mathbf{p}}&= \sum_{r=1}^l 1 - C(\sigma_w^{-1} r^{-1/2} + C \sigma_w^{-2} r ^{-1} + \f 1 2 C^2 \sigma_w^{-3} \f{\log r}{r^{3/2}} + o (\f{\log r}{r^{3/2}}))) + \Theta(r^{-3/2})\\
    &= l - 2 C \sigma_w^{-1} l^{1/2} - C^2 \sigma_w^{-2} \log l + O(1)\end{aligned}$$ which is what we want.

\[lemma:WtPhiAntisymmetric\]

Let $\phi$ is antisymmetric. Then for $\tau \in [0, \pi/2]$, $$\begin{aligned}
    {\mathrm{W}}\phi( {q}, {q}\cos \tau) &= \lim_{t \to \tau}\f 1 {\pi \sin t} \int_{w' \ge |w|}  \dd w \dd w' \Upsilon(w, w';\tau) \phi(\f {\sqrt{{q}}} {\sqrt{2}}(w + w')) \phi(\f {\sqrt {q}} {\sqrt{2}}(w' - w))\\
        &= \f 1 \pi \int_{0}^\infty r\dd r e^{-r^2/2} \int_{0}^{\pi}\dd \theta  \Sigma(\sqrt {q}r, \theta; \tau)\\
        &= \f 1 \pi \int_0^\infty s \dd s \inv {q}e^{-s^2 \inv {q}/2} \int_0^\pi \dd \theta \Sigma(s, \theta; \tau)\\
        &= \f 1 \pi \int_0^\pi \dd \theta \int_0^\infty \dd s e^{-s^2 \inv {q}/2} \pdf{}{s} \Sigma(s, \theta; \tau)
    \end{aligned}$$ where $\Upsilon(w, w';\tau) :=e^{-\f 1 2(\f{w^2}{1-c} + \f{(w')^2}{1+c})} - e^{-\f 1 2(\f{(w')^2}{1-c} + \f{w^2}{1+c})}$ with $c = \cos \tau$, and $\Sigma(s, \theta; \tau) := \phi(s \sin \theta) \phi(s \sin(\theta - \tau)).$

Of course, in the above lemma, the limit in the first equation is only necessary when $\tau = 0$ or $\tau = \pi/2$.

Let $c := \cos \tau$ and $$\Gamma := {\mathrm{W}}\phi( {q}, cq) = \f 1 {2\pi {q}\sqrt{1 - c^2}} \int \dd \mathbf z \exp(-\mathbf z^T \Sigma^{-1} \mathbf z /2) \phi(z)\phi(z'),$$ where $\Sigma = \begin{pmatrix}{q}& cq \\ cq & {q}\end{pmatrix}$.

Our proof will have two portions: Symmetrization of the $\Gamma$ integral and trigonometric change of variables for evaluation.

**Symmetrization.** $\Sigma$ is diagonalized by $\Omega = \f 1 {\sqrt{2q}}\begin{pmatrix}-1 & 1 \\ 1 & 1 \end{pmatrix}$, $$\Sigma = \Omega^T \Diag(1 - c, 1 + c) \Omega.$$ By a change of variable $\mathbf w = \Omega \mathbf z$, so that $\dd \mathbf w = \inv {q}\dd \mathbf z$, we have $$\begin{aligned}
    \Gamma &= \f 1 {2\pi \sqrt{1 - c^2}}\int \dd \mathbf w \exp(-\mathbf w^T \Diag(1 - c, 1 + c)^{-1}\mathbf w/2) \phi(\f {\sqrt {q}} {\sqrt{2}}(w' - w)) \phi(\f {\sqrt {q}} {\sqrt{2}}(w + w'))\\
    &= \f 1 {2\pi \sqrt{1 - c^2}} \int  \dd w \dd w' e^{-\f 1 2(\f{w^2}{1-c} + \f{(w')^2}{1+c})} \phi(\f {\sqrt{{q}}} {\sqrt{2}}(w' - w)) \phi(\f {\sqrt {q}} {\sqrt{2}}(w + w'))
    \end{aligned}$$

By a change of variable swapping $w$ with $w'$, we get

$$\begin{aligned}
    \Gamma &= -\f 1 {2\pi \sqrt{1 - c^2}} \int  \dd w \dd w' e^{-\f 1 2(\f{(w')^2}{1-c} + \f{w^2}{1+c})} \phi(\f {\sqrt{{q}}} {\sqrt{2}}(w + w')) \phi(\f {\sqrt {q}} {\sqrt{2}}(w' - w))
    \end{aligned}$$

Thus $$\begin{aligned}
    2 \Gamma &= \f 1 {2\pi \sqrt{1 - c^2}} \int  \dd w \dd w' \Upsilon(w, w';\tau) \phi(\f {\sqrt{{q}}} {\sqrt{2}}(w + w')) \phi(\f {\sqrt {q}} {\sqrt{2}}(w' - w))
    \end{aligned}$$ where $$\Upsilon(w, w';\tau) = e^{-\f 1 2(\f{w^2}{1-c} + \f{(w')^2}{1+c})} - e^{-\f 1 2(\f{(w')^2}{1-c} + \f{w^2}{1+c})}.$$

![The integrand of $\Gamma$ after symmetrization. Here $c = .2$ and ${q}= 100$ and $\phi = \tanh$.[]{data-label="fig:symmetrization_of_Gamma"}](graphics/symmetrized_integrand.pdf){height=".2\textheight"}

Note that, by the antisymmetry of $\phi$, the integrand $K := \Upsilon(w, w';\tau)\phi(\ldots)\phi(\ldots)$ above has the symmetries $K(w, w') = K(w', w) = K(w, -w')$, and is everywhere nonnegative. \[fig:symmetrization\_of\_Gamma\] displays a contour plot of $K$ for typical values of ${q}$ and $c$. So $$\begin{aligned}
    \Gamma &= \f 1 {\pi \sqrt{1 - c^2}} \int_{w' \ge |w|} \dd w \dd w' K(w, w').
    \end{aligned}$$ This gives the first equation in the lemma.

**Polar Coordinates.** Let $\f w {\sqrt{1 - c}} = r \cos \theta, \f {w'} {\sqrt{1 + c}} = r \sin \theta$, so that $$\begin{aligned}
    w &= r \cos \theta \sqrt{1-c} = \sqrt 2 r \cos \theta \sin \f \tau 2\\
    w' &= r \sin \theta \sqrt{1+c} = \sqrt 2 r \sin \theta \cos \f \tau 2\\
    \dd w \dd w' &= \sqrt{1-c^2}r \dd r \dd \theta = (\sin^2 \tau) r \dd r \dd \theta.
    \end{aligned}$$ Then $$\begin{aligned}
    \mathbf A &:= \int_{w' \ge |w|} e^{-(\f{w^2}{1-c} + \f{(w')^2}{1 +c})/2}\phi(\sqrt{{q}/2}(w+w'))\phi(\sqrt{{q}/2}(w'-w)) \dd w \dd w'\\
    &= \sin^2 \tau \int_{0}^\infty r\dd r e^{-r^2/2} \int_{\tau/2}^{\pi - \tau/2}\dd \theta  \phi(\sqrt {q}r \sin(\theta + \tau/2))\phi(\sqrt {q}r \sin(\theta - \tau/2)).
    \end{aligned}$$ Similarly, let $\f w {\sqrt{1 + c}} = r \cos \theta, \f {w'} {\sqrt{1 - c}} = r \sin \theta$, so that $$\begin{aligned}
    w &= r \cos \theta \sqrt{1+c} = \sqrt 2 r \cos \theta \cos \f \tau 2\\
    w' &= r \sin \theta \sqrt{1-c} = \sqrt 2 r \sin \theta \sin \f \tau 2\\
    \dd w \dd w' &= \sqrt{1-c^2}r \dd r \dd \theta = (\sin^2 \tau) r \dd r \dd \theta,
    \end{aligned}$$ and $$\begin{aligned}
    \mathbf B &= \int_{w' \ge |w|} e^{-(\f{w^2}{1+c} + \f{(w')^2}{1 - c})/2}\phi(\sqrt{{q}/2}(w+w'))\phi(\sqrt{{q}/2}(w'-w)) \dd w \dd w'\\
    &= -\sin^2 \tau \int_{0}^\infty r\dd r e^{-r^2/2} \int_{\pi/2 - \tau/2}^{\pi/2 + \tau/2}\dd \theta  \phi(\sqrt {q}r \cos(\theta + \tau/2))\phi(\sqrt {q}r \cos(\theta - \tau/2))\\
    &= -\sin^2 \tau \int_{0}^\infty r\dd r e^{-r^2/2} \int_{- \tau/2}^{\tau/2}\dd \theta  \phi(\sqrt {q}r \sin(\theta + \tau/2))\phi(\sqrt {q}r \sin(\theta - \tau/2)).
    \end{aligned}$$ Thus $$\begin{aligned}
    \Gamma &= \f 1 {\pi \sqrt{1-c^2}} (\mathbf A - \mathbf B)\\
        &= \f 1 \pi \int_{0}^\infty r\dd r e^{-r^2/2} \int_{- \tau/2}^{\pi - \tau/2}\dd \theta  \phi(\sqrt {q}r \sin(\theta + \tau/2))\phi(\sqrt {q}r \sin(\theta - \tau/2))\\
        &= \f 1 \pi \int_{0}^\infty r\dd r e^{-r^2/2} \int_{0}^{\pi}\dd \theta  \phi(\sqrt {q}r \sin(\theta))\phi(\sqrt {q}r \sin(\theta - \tau)).
    \end{aligned}$$ This gives the second equation in the lemma, and a change of variables $s = \sqrt {q}r$ gives the third.

For the fourth equality, we start from the third equality, and apply integration by parts: $$\begin{aligned}
    &\phantom{{}={}}\f 1 \pi \int_0^\infty s \dd s \inv {q}e^{-s^2 \inv {q}/2} \int_0^\pi \dd \theta \Sigma(s, \theta; \tau)\\
    &= \f 1 \pi \int_0^\pi \dd \theta \int_0^\infty \dd s s \inv {q}e^{-s^2 \inv {q}/2}  \Sigma(s, \theta; \tau)\\
    &= \f 1 \pi \int_0^\pi \dd \theta \biggl(-e^{-s^2 \inv {q}/2} \Sigma(s, \theta; \tau)\biggr\rvert_{s=0}^\infty + \int_0^\infty \dd s e^{-s^2 \inv {q}/2} \pdf{}{s} \Sigma(s, \theta; \tau) \biggr)\\
    &= \f 1 \pi \int_0^\pi \dd \theta \int_0^\infty \dd s e^{-s^2 \inv {q}/2} \pdf{}{s} \Sigma(s, \theta; \tau).
    \end{aligned}$$ where the last equality follows because $\Sigma(0, \theta; \tau) = 0$ and $e^{-s^2\inv {q}/ 2} \to 0$ as $s \to \infty$.

In the following lemmas, the “2” is not important, and can be any arbitrary finite or infinite value.

Suppose a function $f: (0, 2) \to \R$ is $C^k$ on $(0,2)$. If $\lim_{x \downarrow 0} \p f i (x)$ exists and is finite for every $i \in [0, k]$, then $f$ can be extended to $[0, 2)$ such that one sided $i$th derivatives exist at 0 for all $i \in [0, k]$.

Consider $\overline{\p f i}(0) := \p f i (1) - \int_0^1 \p f {i+1}(x) \dd x$ for $i \in [0, k-1]$, which naturally is also equal to $\p f i ({\epsilon}) - \int_0^{\epsilon}\p f {i+1}(x)\dd x$ for any ${\epsilon}> 0$. Certainly $\p f i (x) \to \overline{\p f i}(0)$ as $x \to 0$ if this limit exists — and by assumption it does, for $0 \le i \le k -1$. Therefore, we can define the extension of $\p f i$ to $x = 0$ to be $\p f i (0) := \overline{\p f i}(0)$. But we need to check that for $i \in [0, k-1]$. $$\begin{aligned}
    \lim_{{\epsilon}\to 0} \f 1 {\epsilon}(\p f i({\epsilon}) - \p f i (0)) = \p f {i+1}(0)
    \end{aligned}$$ so that all one sided $i$th derivatives exist. But $$\begin{aligned}
    \f 1 {\epsilon}(\p f i({\epsilon}) - \p f i (0)) &= \f 1 {\epsilon}\int_0^{\epsilon}\p f {i+1}(x)\dd x\\
        &=\p f {i+1}(0) + \int_0^1 (\p f {i+1}(x) - \p f i (0)) \I(x \in [0, {\epsilon}]) \dd x
    \end{aligned}$$ Since $\lim_{x \downarrow 0} \p f {i+1}(x) = \p f {i+1} (0)$, $\p f {i+1}(x) - \p f {i+1}(0)$ is bounded for small $x$, and by dominated convergence, $\int_0^1 (\p f {i+1}(x) - \p f i (0)) \I(x \in [0, {\epsilon}]) \dd x \to \int_0^1 0 \dd x = 0$ as ${\epsilon}\to 0$. Thus $$\begin{aligned}
    \lim_{{\epsilon}\to 0} \f 1 {\epsilon}(\p f i({\epsilon}) - \p f i (0)) &= \p f {i+1}(0)
    \end{aligned}$$ as desired.

If $f: [0, 2) \to \R$ is $C^k$ on $(0, 2)$ and has one sided derivatives at 0 up to order $k$, then $$\begin{aligned}
    f({\epsilon}) &= f(0) + {\epsilon}\p f 1 (0) + \cdots + \f{{\epsilon}^{i-1}}{(i-1)!} \p f {i-1} (0) + O({\epsilon}^i)
    \end{aligned}$$ for any $i \le k$.

We have $$\begin{aligned}
    f({\epsilon}) &= f(0) + \int_0^{\epsilon}\p f 1 (x) \dd x\\
        &= f(0) + {\epsilon}\p f 1(0) + \int_0^{\epsilon}\p f 1(x) - \p f 1 (0) \dd x\\
        &= f(0) + {\epsilon}\p f 1(0) + \int_0^{\epsilon}\int_0^{x_0} \p f 2 (x_2) \dd x_2 \dd x_1\\
        &= f(0) + {\epsilon}\p f 1(0) + \f{{\epsilon}^2}{2} \p f 2(0) + \int_0^{\epsilon}\int_0^{x_1} \p f 2 (x_2) - \p f 2 (0) \dd x_2 \dd x_1\\
        &\ \ \mathbin{\vdots}\ \\
    f({\epsilon}) &= f(0) + {\epsilon}\p f 1 (0) + \cdots + \f{{\epsilon}^{i-1}}{(i-1)!} \p f {i-1} (0) + \int_0^{\epsilon}\dd x_1 \int_0^{x_1} \dd x_2 \cdots \int_0^{x_{i-1}} \dd x_{i} \p f {i}(x_{i})
    \end{aligned}$$ for any $i\le k$. It suffices then to bound the size of the integral. Since $\p f i (x) \to \p f i (0)$ as $x \downarrow 0$ by assumption, $|\p f i (x_i)|$ is bounded by some constant $C$ on the integration region $\mathbb{A} := \{(x_1, \ldots, x_i): {\epsilon}\ge x_1 \ge \cdots \ge x_i\}$ for small enough ${\epsilon}$. Therefore, $$\begin{aligned}
    &\phantom{{}={}}\int_0^{\epsilon}\dd x_1 \int_0^{x_1} \dd x_2 \cdots \int_0^{x_{i-1}} \dd x_{i} \p f {i}(x_{i})\\
    &= \int \p f {i}(x_i) \I(\vec x \in \mathbb A) \dd \vec x\\
    &\le C |\mathbb A|\\
    &= \Theta({\epsilon}^i).
    \end{aligned}$$

As a corollary,

\[lemma:smoothExtensionBoundary\] If $f: (0, 2) \to \R$ is smooth on $(0, 2)$ and $\lim_{x \to 0}\p f i (x)$ exists and is finite for all $i$, then $f$ can be extended to $[0, 2)$ and be [*one-sided smooth*]{} at 0, and $$\begin{aligned}
    f({\epsilon}) &= f(0) + {\epsilon}\p f 1 (0) + \cdots + \f{{\epsilon}^{i-1}}{(i-1)!} \p f {i-1} (0) + O({\epsilon}^i)
    \end{aligned}$$ for any $i$.

\[lemma:WtExtension\]

Let $\phi = \tanh$. For any fixed $c$, ${\mathrm{W}}\phi( {q}, cq)$ is smooth (infinitely differentiable) on ${q}\in (0, \infty)$. As a function of $Q := \inv {q}$, it can be extended smoothly to the point $Q = 0$, so that $$\begin{aligned}
    {\mathrm{W}}\phi( {q}, cq ) &= \lim_{{q}' \to \infty} {\mathrm{W}}\phi( {q}', cq') + \inv {q}\lim_{{q}' \to \infty} \pd {\mathrm{W}}\phi( {q}', cq')/\pd ({q}')^{-1} + \cdots \\
        &\qquad+ \f{{q}^{-i+1}}{(i-1)!} \lim_{{q}' \to \infty} \pd^{i-1} {\mathrm{W}}\phi( {q}', cq')/\pd ({q}')^{-i+1} + O({q}^{-i})
    \end{aligned}$$ for any $i \ge 0$. Furthermore, for $c$ bounded away from $1$, the constants hidden $O$ can be taken independent of $c$.

**Smoothness on $(0, \infty)$.** By the third equation of \[lemma:WtPhiAntisymmetric\], for $Q \in (0, \infty) \iff {q}\in (0, \infty)$, $$\begin{aligned}
    &\phantom{{}={}}\f 1 \pi \int_0^\infty s \dd s \left|\pdf{^n}{Q^n}\lp Q e^{-s^2 Q/2}\rp\right| \int_0^\pi \dd \theta |\phi(s \sin \theta)\phi(s \sin(\theta - \tau))|\\
    &\le \int_0^\infty s \dd s \left|\pdf{^n}{Q^n}\lp Q e^{-s^2 Q/2}\rp\right| < \infty,
    \end{aligned}$$ so by Leibniz’s integral rule and a simple induction, all derivatives of ${\mathrm{W}}\phi( {q}, cq)$ against $Q$ exists for any $Q \in (0, \infty)$.

**Extension to $Q = 0$.** By \[lemma:smoothExtensionBoundary\], it suffices to show that the limit of $
    \pdf{^k{\mathrm{W}}\phi( {q}, cq)}{Q^k}$ exists and is finite as $Q \to 0$, for all $k$. Let $\tau = \arccos c$. By the fourth equation of \[lemma:WtPhiAntisymmetric\], we have explicitly $$\begin{aligned}
    \pdf{^k{\mathrm{W}}\phi( {q}, cq)}{Q^k}
    &= \f 1 \pi \int_0^\pi \dd \theta \int_0^\infty \dd s (\nicefrac{-s^2}{2})^k e^{-s^2 Q /2} \pdf{}{s} \Sigma(s, \theta; \tau)\\
    &= \f {(-2)^{-k}} \pi \int_0^\pi \dd \theta \int_0^\infty \dd s\ s^{2k} e^{-s^2 Q /2} \pdf{}{s} \Sigma(s, \theta; \tau)
    \end{aligned}$$ for any $Q \in (0, \infty)$. Note that for $\phi = \tanh$, $\dot \phi = {\operatorname{sech}}^2$, $$\begin{aligned}
    \pdf{}{s} \Sigma(s, \theta; \tau)
        &= \sin \theta \dot \phi(s \sin \theta)\phi(s \sin (\theta - \tau)) + \sin(\theta-\tau)\phi(s \sin \theta)\dot \phi(s \sin (\theta - \tau)).\\
    \end{aligned}$$ We split the integral of $\pdf{^k {\mathrm{W}}\phi}{Q^k}$ as follows: $$\begin{aligned}
    \pdf{^k{\mathrm{W}}\phi( {q}, cq)}{Q^k}
    &= \f {(-2)^{-k}} \pi \int_0^\pi \dd \theta \int_0^\infty \dd s\ s^{2k} e^{-s^2 Q /2} \sin \theta \dot \phi(s \sin \theta)\phi(s \sin (\theta - \tau))\\
    &\qquad + \f {(-2)^{-k}} \pi \int_0^\pi \dd \theta \int_0^\infty \dd s\ s^{2k} e^{-s^2 Q /2} \sin(\theta-\tau)\phi(s \sin \theta)\dot \phi(s \sin (\theta - \tau))
    \end{aligned}$$ We show that for each piece, the limit as $Q \to 0$ exists and is finite, for any $k$. This will prove the smooth extendability of ${\mathrm{W}}\phi$ to $Q = 0$. We will do this for the first piece; the second is similar.

For $Q > 0$, the integrand is absolutely integrable, so we may switch the integrals. $$\begin{aligned}
    &\phantom{{}={}} \int_0^\pi \dd \theta \int_0^\infty \dd s\ s^{2k} e^{-s^2 Q /2} \sin \theta \dot \phi(s \sin \theta)\phi(s \sin (\theta - \tau))\\
    &= \int_0^\infty \dd s s^{2k} e^{-s^2 Q/2}\int_0^\pi \dd \theta \sin \theta \dot \phi(s \sin\theta) \phi(s \sin (\theta - \tau))
    \end{aligned}$$ We now try to bound the inner integral by an exponentially decreasing term $e^{-s \mu}$ for some $\mu$; clearly, by monotone convergence on the outer integral as $Q \to 0$, this would show the limit of the integral exists and is finite.

Because $\phi$ is odd and $\dot \phi$ is even, the inner integrand is negative on $\theta \in [0, \tau)$ and positive on $\theta \in (\tau, \pi]$. We will break up the inner integral as follows, for some fixed ${\epsilon}> 0$ satisfying $\tau - {\epsilon}> 0$ independent of $s$ (recall $\tau \in (0, \pi/2]$). $$\begin{aligned}
    &\phantom{{}={}} \int_0^\pi \dd \theta \sin \theta \dot \phi(s \sin\theta) \phi(s \sin (\theta - \tau))\\
    &= \left(\int_0^{\epsilon}+ \int_{\pi - {\epsilon}}^\pi \right) \dd \theta \sin \theta \dot \phi(s \sin\theta) \phi(s \sin (\theta - \tau)) + \int_{\epsilon}^{\pi - {\epsilon}} \dd \theta \sin \theta \dot \phi(s \sin\theta) \phi(s \sin (\theta - \tau))
    \end{aligned}$$ Now because $\dot \phi(z) = {\operatorname{sech}}^2(z) \le 2 e^{-z}$, and $\sin \theta \ge \sin {\epsilon}$ on $\theta \in [{\epsilon}, \pi - {\epsilon}]$, $$\begin{aligned}
    &\phantom{{}={}}\left|\int_{\epsilon}^{\pi - {\epsilon}} \dd \theta \sin \theta \dot \phi(s \sin\theta) \phi(s \sin (\theta - \tau))\right|\\
    &\le 2\int_{\epsilon}^{\pi - {\epsilon}} \dd \theta \exp(-s \sin{\epsilon}) \\
    &= 2(\pi - 2{\epsilon}) \exp(-s \sin{\epsilon}).
    \end{aligned}$$

For the other part: $$\begin{aligned}
    &\phantom{{}={}}\int_{\pi - {\epsilon}}^\pi \dd \theta \sin \theta \dot \phi(s \sin\theta) \phi(s \sin (\theta - \tau))\\
        &= \int_{\epsilon}^{0} \sin(\pi - \theta) \dot \phi(s \sin \pi - \theta) \phi(s \sin (\pi - \theta - \tau)) \dd(\pi - \theta)\\
        &= \int_0^{\epsilon}\dd \theta \sin\theta \dot \phi(s \sin \theta) \phi(s \sin \theta + \tau)
    \end{aligned}$$ so that $$\begin{aligned}
    &\phantom{{}={}}\left(\int_0^{\epsilon}+ \int_{\pi - {\epsilon}}^\pi \right) \dd \theta \sin \theta \dot \phi(s \sin\theta) \phi(s \sin (\theta - \tau))\\
    &= \int_0^{\epsilon}\dd \theta \sin \theta \dot \phi(s \sin\theta) [\phi(s \sin (\tau + \theta)) - \phi(s \sin(\tau - \theta))]
    \end{aligned}$$ But by intermediate value theorem, $\phi(s \sin (\tau + \theta)) - \phi(s \sin(\tau - \theta)) = 2 \theta \pd \phi(s \sin(\tau + \theta))/\pd \theta |_{\theta = \psi} = 2 \theta \dot \phi(s \sin (\tau + \psi))s \cos(\tau + \psi)$ for some $\psi \in [- \theta, \theta]$. By the assumption on ${\epsilon}$, $\phi(s \sin (\tau + \theta)) - \phi(s \sin(\tau - \theta)) \le 2{\epsilon}\dot \phi(s \sin(\tau - {\epsilon}))s \cos(\tau - {\epsilon}).$ Then $$\begin{aligned}
    &\phantom{{}={}} \int_0^{\epsilon}\dd \theta \sin \theta \dot \phi(s \sin\theta) [\phi(s \sin (\tau + \theta)) - \phi(s \sin(\tau - \theta))]\\
    &\le \int_0^{\epsilon}\dd \theta \sin \theta \dot \phi(s \sin\theta) 2 {\epsilon}\dot \phi(s \sin(\tau - {\epsilon}))s \cos(\tau - {\epsilon})\\
    &\le 2 {\epsilon}\dot \phi(s \sin(\tau - {\epsilon}))s \cos(\tau - {\epsilon}) O(1)
    \end{aligned}$$ Because $\tau -{\epsilon}> 0$ by assumption on ${\epsilon}$, and because $\dot \phi(z) = \exp(-\Theta_+(z))$, this quantity is $\exp(-\Theta_+(z))$, as desired (here $\Theta_+$ denotes a positive quantity).

Thus $$\begin{aligned}
    &\phantom{{}={}} \int_0^\pi \dd \theta \sin \theta \dot \phi(s \sin\theta) \phi(s \sin (\theta - \tau))\\
    &= \left(\int_0^{\epsilon}+ \int_{\pi - {\epsilon}}^\pi \right) \dd \theta \sin \theta \dot \phi(s \sin\theta) \phi(s \sin (\theta - \tau)) + \int_{\epsilon}^{\pi - {\epsilon}} \dd \theta \sin \theta \dot \phi(s \sin\theta) \phi(s \sin (\theta - \tau))\\
    &= \exp(-\Theta_+(s))
    \end{aligned}$$ and similarly for the other piece of $\pdf{^k {\mathrm{W}}\phi}{Q^k}$, so that $$\begin{aligned}
    &\phantom{{}={}} \int_0^\infty \dd s s^{2k} e^{-s^2 Q/2}\int_0^\pi \dd \theta \sin \theta \dot \phi(s \sin\theta) \phi(s \sin (\theta - \tau))\\
    &= \int_0^\infty \dd s s^{2k} e^{-s^2\f Q 2 - \Theta_+(z)}\\
    &\to \int_0^\infty \dd s s^{2k} e^{-\Theta_+(z)}
    \end{aligned}$$ is finite as $Q \to 0$, by monotone convergence.

**Independence of constant hidden in $O(({q}')^{-i})$.** The constant hidden is a function of the ${\epsilon}$ chosen above, which depend on $\tau$, but only to the extent that it must satisfy $\tau - {\epsilon}> 0$. As long as we are interested in a set $\mathcal C$ of $c$ that is bounded away from 1, the corresponding set of $\tau$ is bounded away from $0$, so ${\epsilon}$ can be taken to be some number smaller than all of the corresponding $\tau$.

\[lemma:Wt\_le\_arcsin\]

Suppose $\phi$ is tanh-like. Then for $c \in [0, 1]$, $${\mathrm{W}}\phi( {q}, cq) \le \f 2 \pi \arcsin(c),$$ and weakly increases to this upper bound as ${q}\to \infty$. Furthermore,

-   If $c = 0$ or 1, then equality holds regardless of ${q}$.

-   If $c \in (0, 1)$ is held constant, $\f 2 \pi \arcsin(c) - {\mathrm{W}}\phi( {q}, cq) = \Theta(\inv {q})$, where the hidden constants in $\Theta$ depend on $c$. But the constants can be made independent of $c$ if $c \in [{\epsilon}, 1-{\epsilon}]$ for some ${\epsilon}> 0$.

The cases of $c = 0$ or 1 are obvious by the definition of ${\mathrm{W}}$. So from here on we assume $c \in (0, 1)$.

Let $\tau := \arccos c$. By the first equation of \[lemma:WtPhiAntisymmetric\] and the assumption that $\phi$ is tanh-like, it is immediate that ${\mathrm{W}}\phi({q},cq)$ is nondecreasing in ${q}$. By dominated convergence, using the second equation of \[lemma:WtPhiAntisymmetric\], we get $$\begin{aligned}
    \lim_{{q}\to \infty} {\mathrm{W}}\phi({q},cq) &= \f 1 \pi \int_{0}^\infty r\dd r e^{-r^2/2} (\pi - 2 \tau)\\
        &= \f {\pi - 2 \tau}{\pi}\\
        &= \f 2 \pi \arcsin c.
    \end{aligned}$$

Then the convergence rate is $O(\inv {q})$ by \[lemma:WtExtension\] and Taylor’s theorem. Thus to show the convergence rate is $\Theta(\inv {q})$, it suffices to show that $\mathbf D := \pdf{{\mathrm{W}}\phi( {q}, cq)}{Q} < 0$. But this is apparent from the first equation of \[lemma:WtPhiAntisymmetric\]: For $\tau \in (0, \pi/2)$, $$\begin{aligned}
    \mathbf D &= \f 1 {\pi \sin \tau} \int_{w' \ge |w|} \dd w \dd w' \Upsilon(w, w';\tau)(-\f 1 {2\sqrt 2} Q^{-3/2})\\
        &\qquad\qquad\qquad\quad \times[\dot\phi(\sqrt{{q}/2}(w + w'))\phi(\sqrt{{q}/2}(w'-w))\\
        &\qquad\qquad\qquad\qquad + \phi(\sqrt{{q}/2}(w + w'))\dot\phi(\sqrt{{q}/2}(w'-w))]\\
        & < 0
    \end{aligned}$$ since $\Upsilon$ is positive on the integration domain, and $\dot \phi$ and $\phi$ are both positive for positive arguments, by the assumption of $\phi$ being tanh-like.

**Independence of the constants in $\Theta(\inv {q})$ from $c$ when $c \in [{\epsilon}, 1-{\epsilon}]$.** By \[lemma:WtExtension\], the upper constant can be made independent from $c$. Since $\mathbf D$ is monotonically decreasing in $c$ (or monotonically increasing in $\tau$) and $|\mathbf D|$ is monotonically increasing in $c$ (or monotonically decreasing in $\tau$), we have $|\mathbf D| > |\mathbf D|\bigg]_{c = {\epsilon}}$, which can be taken to be the lower constant in $\Theta(\inv {q})$.

\[fig:verifywtanhasymptotics\] verifies empirically that the subleading term in ${\mathrm{W}}\tanh( q, cq)$ is linear in $\inv q$, for constant $c$.

![We verify empirically that the subleading term in ${\mathrm{W}}\tanh( q, cq)$ is linear in $\inv q$, for constant $c$. Indeed, observe that the curve of of ${\mathrm{W}}\tanh$ intersects the y-axis at an angle.[]{data-label="fig:verifywtanhasymptotics"}](graphics/verify_Wtanh_asymptotics){height=".15\textheight"}

We have by \[lemma:Wt\_le\_arcsin\], $$\begin{aligned}
{\boldsymbol{\oldgamma}}&= \f 2 \pi \arcsin({\boldsymbol{\oldlambda}}/{q}) - \Theta({q}^{-1}) + {\underline}{\boldsymbol{\oldgamma}}.\end{aligned}$$ Since ${q}= \sigma_w^2 {\underline}{\mathbf{p}}+ \sigma_b^2$ by \[thm:p\_q\_linear\], and ${\boldsymbol{\oldlambda}}= \sigma_w^2 {\underline}{\boldsymbol{\oldgamma}}+ \sigma_b^2$ by \[thm:lambda\_gamma\_recurrence\], $$\begin{aligned}
{\boldsymbol{\oldgamma}}&= \f 2 \pi \arcsin\lp\f{\sigma_w^2 {\underline}{\boldsymbol{\oldgamma}}+ \sigma_b^2}{\sigma_w^2 {\underline}{\mathbf{p}}+ \sigma_b^2}\rp - \Theta({q}^{-1}) + {\underline}{\boldsymbol{\oldgamma}}.\end{aligned}$$ We claim that $\p {\boldsymbol{\oldgamma}}l \to \infty$ as $l \to \infty$. Otherwise, there is some $C$ such that $\p {\boldsymbol{\oldgamma}}l \le C$ for all $l$. For large enough $l$, $\p {\mathbf{p}}l \ge \omega l$ for any $\omega < 1$ and $\arcsin \lp \f{C}{\sigma_w^2 \p {\mathbf{p}}{l-1} + \sigma_b^2} \rp = \Theta(1/l)$ by linearization of $\arcsin$. Thus $\p {\boldsymbol{\oldgamma}}l = \Theta(\log l)$, but this contradicts our assumption that ${\boldsymbol{\oldgamma}}$ is bounded. This proves our claim.

Therefore, for large enough $l$, $$\f{\sigma_w^2 {\underline}{\boldsymbol{\oldgamma}}+ \sigma_b^2}{\sigma_w^2 {\underline}{\mathbf{p}}+ \sigma_b^2} = {\underline}{\boldsymbol{\oldgamma}}/{\underline}{\mathbf{p}}+ \Theta(l^{-1}).$$

\[fig:arcsin-vs-x\] shows $\f2\pi\arcsin x$ vs $x$. One sees that 1 is an unstable fixed point; if ${\mathbf{e}}< 1 - \epsilon$, then $\f 2 \pi \arcsin {\mathbf{e}}< 1 - \epsilon - \delta$ for some $\delta$. Thus ${c}$ drops monotonically until some threshold under which the linearization of $\arcsin$, $\arcsin x = x + \Theta(x^3)$, is applicable. So for large enough $l$, $$\begin{aligned}
{\boldsymbol{\oldgamma}}- {\underline}{\boldsymbol{\oldgamma}}&= \f 2 \pi \arcsin ({\underline}{\boldsymbol{\oldgamma}}/ {\underline}{\mathbf{p}}+ \Theta(\inv l)) - \Theta(\inv l)\\
    &= \f 2 \pi  {\underline}{\boldsymbol{\oldgamma}}/ {\underline}{\mathbf{p}}+ O(\inv l)\end{aligned}$$ As $\p {\mathbf{p}}l \sim l$ by \[thm:p\_q\_linear\], this difference equation has solution ${\boldsymbol{\oldgamma}}= \Omega(l^{\f 2 \pi -{\epsilon}}), O(l^{\f 2 \pi + {\epsilon}})$ for any ${\epsilon}$ by using the dynamics of \[lemma:alphaDeltaDynamics\] to upper and lower bound this difference equation.

### Backward Dynamics

The $\sigma_w = 0$ case is obvious. We will assume $\sigma_w > 0$ from here on.

Let $\p {\mathbf{p}}l = b_0 l + b_1 l^{1/2} + b_2 \log l + O(1)$. Then for $D = \f 2 3 \sqrt{\f 2 \pi}$, we have (implicitly applying \[lemma:Vt\_dot\_tanh\_asymptotics\] and \[lemma:power\_sum\_asymptotics\]), $$\begin{aligned}
    {q}^{-1/2} &= \sigma_w^{-1} b_0^{-1/2} l^{-1/2}(1 - b_1 b_0^{-1} \inv 2 l^{-1/2} - b_2 b_0^{-1} 2^{-1} l^{-1} \log l + O(l^{-1}))\\
    {\mathrm{V}}\dot \phi( {q}) & = D {q}^{-1/2} + \Theta({q}^{-3/2})\\
        &= D\sigma_w^{-1} b_0^{-1/2} l^{-1/2}(1 - b_1 b_0^{-1} \inv 2 l^{-1/2} - b_2 b_0^{-1} 2^{-1} l^{-1} \log l + O(l^{-1}))\\
    \log(B {\mathrm{V}}\dot \phi( {q}) + 1) &= B D \sigma_w^{-1} b_0^{-1/2} l^{-1/2}\\
                            &\phantom{={}}- (BD \sigma_w^{-1} b_0^{-3/2} b_1 2^{-1} +B^2 D^2 \sigma_w^{-2} b_0^{-1}2^{-1})l^{-1} + \Theta(l^{-3/2}\log l)\\
    \sum_{r=1}^l \log(B {\mathrm{V}}\dot \phi( \p {q}r) + 1) &= 2B D \sigma_w^{-1} b_0^{-1/2} l^{1/2}
\\
        &\phantom{={}} -(BD \sigma_w^{-1} b_0^{-3/2} b_1 2^{-1} +B^2 D^2 \sigma_w^{-2} b_0^{-1}2^{-1})\log l + O(1)\\
    \end{aligned}$$ In our case, we have $b_0 = 1, b_1 = -2C \sigma_w^{-1}, b_2 = C^2 \sigma_w^{-2},  B = \sigma_w^2, C = \sqrt{\f 2 \pi}$, which gives $$\sum_{r=1}^l \log(B {\mathrm{V}}\dot \phi( \p {q}r) + 1) = \f 4 3 \sqrt{\f 2 \pi} \sigma_w l^{1/2} + (\f 4 {3\pi} - \sigma_w^2 \f 4 {9\pi}) \log l + O(1).$$ so that $$\begin{aligned}
    \p {{\boldsymbol{\oldchi}}}{m}/\p {{\boldsymbol{\oldchi}}}l &= \exp\left[\f 4 3 \sqrt{\f 2 \pi} \sigma_w (\sqrt l - \sqrt m) + (\f 4 {3\pi} - \sigma_w^2 \f 4 {9\pi}) (\log l - \log m) + O(1)\right]\\
    \end{aligned}$$

The $\sigma_w = 0$ case is obvious. We will assume $\sigma_w > 0$ from here on.

As in the proof of \[thm:dalethExpSqrtTanh\], $${\mathrm{V}}\dot \phi( {q}) = D\sigma_w^{-1} b_0^{-1/2} l^{-1/2} + \Theta(\inv l)$$ where $D = \f 2 3 \sqrt{\f 2 \pi}$. Thus by \[thm:dalethRecReduced\], $$\log(\p {{\boldsymbol{\oldchi}}}{m}/\p {{\boldsymbol{\oldchi}}}l) = \f 4 3 \sqrt{\f 2 \pi} \sigma_w (\sqrt l - \sqrt m) + (\f 4 {3\pi} - \sigma_w^2 \f 4 {9\pi}) (\log l - \log m) + O(1)$$ $$\log(\p {\boldsymbol{\oldchi}}m_b/\p {\boldsymbol{\oldchi}}l _b) =\f 4 3 \sqrt{\f 2 \pi} \sigma_w (\sqrt l - \sqrt m) + (\f 4 {3\pi} - \f 1 2 - \sigma_w^2 \f 4 {9\pi}) (\log l - \log m) + O(1)$$ Similarly, since ${\mathbf{p}}= l + \Theta(\sqrt l)$ by \[thm:p\_q\_linear\], we have $$\log(\p {\boldsymbol{\oldchi}}m_w/\p {\boldsymbol{\oldchi}}l _w) = \f 4 3 \sqrt{\f 2 \pi} \sigma_w (\sqrt l - \sqrt m) + (\f 4 {3\pi} + \f 1 2 - \sigma_w^2 \f 4 {9\pi}) (\log l - \log m) + O(1)$$

Tanh: Full Residual Network
---------------------------

### Forward Dynamics

The $\sigma_w = 0$ case is obvious. We will assume $\sigma_w > 0$ from here on.

As in \[thm:p\_q\_linear\], ${\mathbf{p}}$ will have expansion ${\mathbf{p}}= b_0 l + b_1 l^{1/2} + b_2 \log l + O(1).$ Then, for $C = \sqrt{\f 2 \pi}$, $$\begin{aligned}
    {q}^{-1/2} &= \sigma_w^{-1} b_0^{-1/2} l^{-1/2}(1 - b_1 b_0^{-1} \inv 2 l^{-1/2} - b_2 b_0^{-1} 2^{-1} l^{-1} \log l + O(l^{-1}))\\
    \sum_{r=1}^l {\mathrm{V}}\phi( \p {q}r) &= \sum_{r=1}^l 1 - C(\p {q}r)^{-1/2} + \Theta((\p {q}r)^{-3/2})\\
        &= l - 2 C \sigma_w^{-1} b_0^{-1/2} l^{1/2} + C \sigma_w^{-1} b_1 b_0^{-3/2} 2^{-1} \log l + O(1)\\
    \p {\mathbf{p}}l &= \sigma_v^2 \sum_{r=1}^l + \sigma_a^2 l\\
        &= (\sigma_v^2 + \sigma_a^2) l - 2C \sigma_v^2 \sigma_w^{-1} b_0^{-1/2} l^{1/2} + C \sigma_v^2 \sigma_w^{-1} b_1 b_0^{-3/2} 2^{-1} \log l + O(1)
    \end{aligned}$$ which yields $$\begin{aligned}
    b_0 &= \sigma_v^2 + \sigma_a^2\\
    b_1 &= -2C \sigma_v^2 \sigma_w^{-1} b_0^{-1/2} = \f{-2C \sigma_v^2 \sigma_w^{-1}}{\sqrt{\sigma_v^2 + \sigma_a^2}}\\
    b_2 &= \f{-C^2 \sigma_v^4 \sigma_w^{-2}}{(\sigma_v^2 + \sigma_a^2)^2}
    \end{aligned}$$

Suppose $\phi$ is tanh-like. Then $${\boldsymbol{\oldgamma}}\le \sigma_v^2 \f 2 \pi \arcsin\lp { {\boldsymbol{\oldlambda}}}/{{q}}\rp + \sigma_a^2 + {\underline}{\boldsymbol{\oldgamma}},$$ and $$\sigma_v^2 \f 2 \pi \arcsin\lp { {\boldsymbol{\oldlambda}}}/{{q}}\rp + \sigma_a^2 + {\underline}{\boldsymbol{\oldgamma}}- {\boldsymbol{\oldgamma}}= \Theta({q}^{-1}).$$

Similar to the proof of \[lemma:Wt\_le\_arcsin\].

\[lemma:timeDependentConvergence\] Let $u^* \in [0, 1)$. Let $f_t: [0, 1) \to [0, 1]$ be a continuous function for each $t \in \N$, to each of which we associate two numbers $0 \le a_t \le u^* \le b_t \le 1$. Suppose for each $t$, $f_t(u) > u$ for all $u \in [0, a_t)$ and $f_t(u) < u$ for all $u \in (b_t, 1)$. Assume that for each $u$, $f_t(u) - u \to 0$ as $t \to \infty$ uniformly over $u$. If $a_t \nearrow u^*$ and $b_t \searrow u^*$, then for any $u_0 \in [0, 1)$, the dynamics $u_t = f_t(u_{t-1})$ has a limiting point. Furthermore, either $u_t \to u^*$ or $u_t$ eventually converges monotonically (decreasing or increasing) to a limit point.

Fix a $u_0 \in [0, 1)$. If $u_t \to u^*$ then we are done. Otherwise, suppose there is a neighborhood $[u^* - {\epsilon}, u^* + {\epsilon}]$ such that for an infinite sequence $t_1, t_2, \ldots$, $u_{t_i} \not \in [u^* - {\epsilon}, u^* + {\epsilon}]$. WLOG assume $u_{t_i} < u^* - {\epsilon}$ for all $i$ and $(t_i)_i$ is the sequence of all $t$s that satisfy this inequality.

If $(t_i)_i$ contains $\{s: s \ge N\}$ for some $N$, then for some $M > N$, for every $t > M$, $a_t > u^* - {\epsilon}> u_t$. By assumption, $u_t$ is monotonic for all $t > M$ but is bounded above. Thus $u_t$ has a fixed point $\hat u \le u^* - {\epsilon}$ as desired.

Now assume there are infinite $i$s such that $t_i - 1 \not= t_{i-1}$ (i.e. $t_i - 1$ is not part of the sequence $(t_i)_i$). We will show that this case is contradictory. Take $T$ large enough such that $a_t > u^* - {\epsilon}/2$ and $|f_t(u) - u| < {\epsilon}/4$ for all $u$ and for all $t \ge T$ ($T$ exists by premise). Let $j$ be the smallest index such that $t_j > T$ and $t_j - 1 \not = t_{j-1}$. By the definition of $j$, $u_{t_j - 1} \ge u^* - {\epsilon}$. If $u_{t_j - 1} \ge u^* - {\epsilon}/2$, then by definition of $T$, $u^* - {\epsilon}> u_{t_j} = f_{t_j}(u_{t_j - 1}) > u_{t_j-1} - {\epsilon}/4 > u^* - 3{\epsilon}/4 > u^* - {\epsilon}$, a contradiction. If $u^* - {\epsilon}\le u_{t_j - 1} \le u^* - {\epsilon}/2$, then by the definition of $T$, $u_{t_j-1} \le a_{t_j-1}$ so that $u_{t_j } = f_{t_j}(u_{t_j - 1}) > u_{t_j - 1} \ge u^* - {\epsilon},$ a contradiction.

The “furthermore” claim is clear from our proof above.

The $\sigma_w = 0$ case is obvious. We will assume $\sigma_w > 0$ from here on.

If $\sigma_a = 0$, then ${\mathbf{e}}^*$ as defined above is 0, and ${\mathbf{e}}= \f {\boldsymbol{\oldgamma}}{\mathbf{p}}$ decreases as $\Theta(l^{\f 2 \pi - 1})$ to 0, by the same reason as before.

So from now on suppose $\sigma_a > 0$. We apply \[lemma:timeDependentConvergence\] first to show that ${\mathbf{e}}$ converges. We have $$\begin{aligned}
\sigma_v^2 {\mathrm{W}}\phi( {q}, cq) + \sigma_a^2 &= {\mathbf{e}}{\mathbf{p}}- {\underline}{\mathbf{e}}{\underline}{\mathbf{p}}\\
    &= {\mathbf{e}}{\mathbf{p}}- {\underline}{\mathbf{e}}{\mathbf{p}}+ {\underline}{\mathbf{e}}{\mathbf{p}}- {\underline}{\mathbf{e}}{\underline}{\mathbf{p}}\\
    &= ({\mathbf{e}}- {\underline}{\mathbf{e}}) {\mathbf{p}}+ {\underline}{\mathbf{e}}({\mathbf{p}}- {\underline}{\mathbf{p}})\\
    &= ({\mathbf{p}}- {\underline}{\mathbf{p}})[({\mathbf{e}}- {\underline}{\mathbf{e}})\f {\mathbf{p}}{{\mathbf{p}}- {\underline}{\mathbf{p}}} + {\underline}{\mathbf{e}}]\\
\f{\sigma_v^2 {\mathrm{W}}\phi( {q}, cq) + \sigma_a^2}{\sigma_v^2 {\mathrm{V}}\phi( {q}) + \sigma_a^2} &= ({\mathbf{e}}- {\underline}{\mathbf{e}}) \f {\mathbf{p}}{{\mathbf{p}}- {\underline}{\mathbf{p}}} + {\underline}{\mathbf{e}}\\
\f {{\mathbf{p}}- {\underline}{\mathbf{p}}}{\mathbf{p}}\left[\f{\sigma_v^2 {\mathrm{W}}\phi + \sigma_a^2}{\sigma_v^2 {\mathrm{V}}\phi + \sigma_a^2} - {\underline}{\mathbf{e}}\right] &= {\mathbf{e}}- {\underline}{\mathbf{e}}\end{aligned}$$ If we define $f_l(u) := \f{\p {\mathbf{p}}l - \p {\mathbf{p}}{l-1}}{\p {\mathbf{p}}l}\left[\f{\sigma_v^2 {\mathrm{W}}\phi( \p {q}l, \p {c}l \p {q}l) + \sigma_a^2}{\sigma_v^2 {\mathrm{V}}\phi( \p {q}l) + \sigma_a^2} - u\right] + u$ (the LHS of the above), then $f_l(u) - u = O(l^{-1})$ uniformly for all $u$ because $\p {\mathbf{p}}l = \Theta(l)$, $\p {\mathbf{p}}l - \p {\mathbf{p}}{l-1} = \Theta(1)$, and the part in the bracket is $O(1)$, with constants all (able to be taken) independent of $u$. We divide $[0, 1)$ into the following intervals $I_1 = [1, 1/2), I_2 = [1/2, 3/4), I_3 = [3/4, 7/8), \ldots$. For each $I_k$, it is clear that the trajectories of $\p {\mathbf{e}}l = f_l(\p {\mathbf{e}}{l-1})$ with $\p {\mathbf{e}}0 \in I_k$ will fall into some interval $J_k$ bounded away from 1 for all $l \ge L$, for large enough $L$ (dependent on $k$). Then we can apply \[lemma:cExpansion,lemma:vtanhSqrtConvergence,lemma:Wt\_le\_arcsin\] to get $f_l(u) = \f{\p {\mathbf{p}}l - \p {\mathbf{p}}{l-1}}{\p {\mathbf{p}}l}\left[\f{\sigma_v^2\f 2 \pi \arcsin(u) + \sigma_a^2}{\sigma_v^2 + \sigma_a^2} - u + o(1)\right] + u$ where the constants in $o(1)$ is uniform for all $\p {\mathbf{e}}0 \in I_k$. For $u < {\mathbf{e}}^*$ (as defined in the theorem statement), $\f{\sigma_v^2\f 2 \pi \arcsin(u) + \sigma_a^2}{\sigma_v^2 + \sigma_a^2} > u$ and for $u > {\mathbf{e}}^*$, $\f{\sigma_v^2\f 2 \pi \arcsin(u) + \sigma_a^2}{\sigma_v^2 + \sigma_a^2} < u$ (see \[fig:arcsin-vs-x\]). Thus as $l \to \infty$, the $o(1)$ term gets smaller and smaller, and this monotonicity holds for $f_l(u) - u = \left[\f{\sigma_v^2\f 2 \pi \arcsin(u) + \sigma_a^2}{\sigma_v^2 + \sigma_a^2} - u + o(1)\right] > 0$ (resp. $<0$) on larger and larger intervals $[0, a_l]  \cap J_k$ (resp. $[b_l, 1) \cap J_k$). This proves all the preconditions for \[lemma:timeDependentConvergence\], which yields that $I_k$ converges to a limit point. As this argument is independent of $k$, we have that for all $\p {\mathbf{e}}0 \in [0, 1)$, $\p {\mathbf{e}}l$ converges.

Now we solve for the limit point.

Suppose ${\mathbf{e}}$ has limit point ${\mathbf{e}}^\dagger$ (possibly different from ${\mathbf{e}}^*$ described in the theorem); if we express $\p {\boldsymbol{\oldgamma}}l = ({\mathbf{e}}^\dagger + \p {\epsilon}l ) \p {\mathbf{p}}l$, then $$\begin{aligned}
\sigma_v^2 {\mathrm{W}}\phi( {q}, cq) + \sigma_a^2 &= {\boldsymbol{\oldgamma}}- {\underline}{\boldsymbol{\oldgamma}}\\
                &= ({\mathbf{e}}^\dagger + {\epsilon}){\mathbf{p}}- ({\mathbf{e}}^\dagger + {\underline}{\epsilon}) {\underline}{\mathbf{p}}\\
                &= {\mathbf{e}}^\dagger({\mathbf{p}}- {\underline}{\mathbf{p}}) + {\epsilon}{\mathbf{p}}- {\underline}{\epsilon}{\underline}{\mathbf{p}}\\
\f{\sigma_v^2 {\mathrm{W}}\phi( {q}, cq) + \sigma_a^2}{\sigma_v^2 {\mathrm{V}}\phi( {q}) + \sigma_a^2} &= {\mathbf{e}}^\dagger + {\epsilon}+ ({\epsilon}- {\underline}{\epsilon})\f{{\underline}{\mathbf{p}}}{{\mathbf{p}}- {\underline}{\mathbf{p}}}\end{aligned}$$ As $l \to \infty$, ${c}\sim {\mathbf{e}}\to {\mathbf{e}}^\dagger$, and ${\mathrm{W}}\phi( {q}, {\mathbf{e}}^\dagger {q}) \to \f 2 \pi \arcsin({\mathbf{e}}^\dagger)$, and ${\mathrm{V}}\phi( {q}) \to 1$. Additionally, ${\underline}{\mathbf{p}}/({\mathbf{p}}- {\underline}{\mathbf{p}}) = \Theta(l)$ and ${\epsilon}= o(1)$ so that ${\epsilon}- {\underline}{\epsilon}= o(l^{-1})$. Then we have, taking limits $l \to \infty$, $$\begin{aligned}
\f{\sigma_v^2 \f 2 \pi \arcsin({\mathbf{e}}^\dagger)+\sigma_a^2}{\sigma_v^2 + \sigma_a^2} &= {\mathbf{e}}^\dagger.\end{aligned}$$ Since $f_l$ (as defined above) repels points away from 1, the only solution for ${\mathbf{e}}^\dagger$ when $\p {\mathbf{e}}0 < 1$ is ${\mathbf{e}}^\dagger = {\mathbf{e}}^*$ as specified in the theorem statement.

We defer the proof of the convergence rate to ${\mathbf{e}}^*$ to \[thm:TanhFullResConvergenceRate\].

![Graph of $y({\mathbf{e}}) = \f 1 {\sigma_v^2 + \sigma_a^2}[\sigma_v^2\f 2 \pi \arcsin({\mathbf{e}}) + \sigma_a^2]$ for various $\sigma_v$ and $\sigma_a$.[]{data-label="fig:arcsin-vs-x"}](graphics/arcsin_vs_x.pdf){height=".15\textheight"}

\[lemma:fixed\_point\_ineq\] Let ${\mathbf{e}}^*$ be the stable fixed point determined by $\sigma_a$ and $\sigma_v$. Then as long as $\sigma_v > 0$, $$\begin{aligned}
\f 2 \pi \f 1 {\sqrt{1 - ({\mathbf{e}}^*)^2}} \f{\sigma_v^2 }{\sigma_v^2 + \sigma_a^2} \in (\f 1 2, \f 2 \pi]\end{aligned}$$

Write $\rho := \f{\sigma_a^2}{\sigma_v^2}$. By definition of ${{{\mathbf{e}}^*}}$, we get $$\begin{aligned}
    {{{\mathbf{e}}^*}}&= (1 - \rho) \f 2 \pi \arcsin {{{\mathbf{e}}^*}}+ \rho\\
    \rho = &= \f{{{{\mathbf{e}}^*}}- \f 2 \pi \arcsin {{{\mathbf{e}}^*}}}{1 - {{{\mathbf{e}}^*}}}
    \end{aligned}$$ Substituting $\rho$ into the expression in question, it follows that we want to show $$\begin{aligned}
    \f 2 \pi (1 - {{{\mathbf{e}}^*}}^2)^{-1/2} (1 + \rho)^{-1}=  \f 2 \pi (1 - {{{\mathbf{e}}^*}}^2)^{-1/2} \left(\f{1 - \f 2 \pi \arcsin {{{\mathbf{e}}^*}}}{1 - {{{\mathbf{e}}^*}}}\right)^{-1}\in (\f 1 2, \f 2 \pi]
    \end{aligned}$$ for ${{{\mathbf{e}}^*}}\in [0, 1)$ (the endpoint at 1 is not included since $\sigma_v > 0$. But this is $$\begin{aligned}
    &\phantom{={}} \f 2 \pi (1 - {{{\mathbf{e}}^*}})^{1/2} (1 + {{{\mathbf{e}}^*}})^{-1/2} (1 - \f 2 \pi \arcsin {{{\mathbf{e}}^*}})^{-1}.
    \end{aligned}$$ Set $g({{{\mathbf{e}}^*}})$ to be this expression. We could proceed by finding critical points, but a simple plot \[fig:deltastarBound\] shows that $g$ is decreasing on $[0, 1)$, with extremal values at the end points: $$g({{{\mathbf{e}}^*}}) \in [\lim_{{{{\mathbf{e}}^*}}\to 1} g({{{\mathbf{e}}^*}}), g(0)), \quad \text{for }{{{\mathbf{e}}^*}}\in [0, 1).$$ Obviously $g(0) = \f 2 \pi$. For the limit, we note that $\arcsin {{{\mathbf{e}}^*}}$ has an asymptotic expansion $\f \pi 2 - \sqrt 2 (1 - e)^{1/2} + \Theta((1-e)^{3/2})$ at 1, so that $(1 - {{{\mathbf{e}}^*}})^{1/2}(1 - \f 2 \pi \arcsin {{{\mathbf{e}}^*}})^{-1} \to \dfrac{\pi}{2 \sqrt 2}$, and $g({{{\mathbf{e}}^*}}) \to \dfrac{1}{2}$ as ${{{\mathbf{e}}^*}}\to 1$.

![Plot of $g({{{\mathbf{e}}^*}})$ in the proof of \[lemma:fixed\_point\_ineq\][]{data-label="fig:deltastarBound"}](graphics/deltastarBound.pdf){width=".4\textwidth"}

\[thm:TanhFullResConvergenceRate\] If $\p {\mathbf{e}}0 < 1$, then $|\p {\mathbf{e}}l - {\mathbf{e}}^*|$ is $\Omega(l^{-\delta^* - \varepsilon})$ and $O(l^{-\delta^* + \varepsilon})$ for any $\varepsilon > 0$, where $$\delta^* := 1 - \f 2 \pi \f 1 {\sqrt{1 - ({\mathbf{e}}^*)^2}} \f{\sigma_v^2 }{\sigma_v^2 + \sigma_a^2} \in [1-\f 2 \pi, \f 1 2),$$ where the bounds on the right follow from \[lemma:fixed\_point\_ineq\].

Define $\omega(q, c) = \f 2 \pi \arcsin(c) - {\mathrm{W}}\tanh( q, c q)$. By \[lemma:Wt\_le\_arcsin\], for large enough $l$, ${c}$ is close to ${\mathbf{e}}^*$ bounded away from 0 or 1, so that $\omega({q}, {c}) = \Theta(\inv {q})$ with the constant hidden in $\Theta$ independent of ${c}$. Additionally, by \[lemma:vtanhSqrtConvergence\], $1 - {\mathrm{V}}\tanh( {q}) = \Theta({q}^{-1/2})$. Therefore, $$\begin{aligned}
({\mathbf{e}}^* + \epsilon){\mathbf{p}}&= \sigma_v^2 (\f 2 \pi \arcsin({\mathbf{e}}^* + {\underline}\epsilon) - \omega({q}, {c}))+ \sigma_a^2 + {\underline}{\boldsymbol{\oldgamma}}\\
&= \sigma_v^2 \f 2 \pi [\arcsin({\mathbf{e}}^*) + \f{{\underline}\epsilon}{\sqrt{1 - ({\mathbf{e}}^*)^2}} + \Theta({\underline}{\epsilon}^2)] - \Theta(\inv l) + \sigma_a^2 + {\underline}{\boldsymbol{\oldgamma}}\\
&= {\mathbf{e}}^*(\sigma_v^2 + \sigma_a^2) + ({\mathbf{e}}^* + {\underline}{\epsilon}) {\underline}{\mathbf{p}}+ \sigma_v^2 \f 2 \pi\f{{\underline}\epsilon}{\sqrt{1 - ({\mathbf{e}}^*)^2}} + \Theta({\underline}{\epsilon}^2) - \Theta(\inv l)\\
{\mathbf{e}}^*({\mathbf{p}}- {\underline}{\mathbf{p}}- \sigma_v^2 - \sigma_a^2) &= {\underline}{\epsilon}{\underline}{\mathbf{p}}- {\epsilon}{\mathbf{p}}+ \sigma_v^2 \f 2 \pi\f{{\underline}\epsilon}{\sqrt{1 - ({\mathbf{e}}^*)^2}} + \Theta({\underline}{\epsilon}^2)- \Theta(\inv l)\\
{\mathbf{e}}^*\sigma_v^2({\mathrm{V}}\phi( {q}) - 1) &=  {\underline}{\epsilon}{\underline}{\mathbf{p}}- {\epsilon}{\mathbf{p}}+ \sigma_v^2 \f 2 \pi\f{{\underline}\epsilon}{\sqrt{1 - ({\mathbf{e}}^*)^2}} + \Theta({\underline}{\epsilon}^2) - \Theta(\inv l)\\
{\epsilon}&= \f 1 {\mathbf{p}}({\mathbf{e}}^*\sigma_v^2(1 - {\mathrm{V}}\phi( {q})) + \Theta({\underline}{\epsilon}^2) - \Theta(\inv l) + {\underline}{\epsilon}({\underline}{\mathbf{p}}+ \sigma_v^2 \f 2 \pi\f{1}{\sqrt{1 - ({\mathbf{e}}^*)^2}}))\\
    &= \Theta(l^{-3/2}) + {\underline}{\epsilon}(1 - \p \delta l / l)\\
$$ where $$\begin{aligned}
\p \delta l &= \f l {\mathbf{p}}(\sigma_v^2 {\mathrm{V}}\phi( {q}) + \sigma_a^2 -\sigma_v^2 \f 2 \pi\f{1}{\sqrt{1 - ({\mathbf{e}}^*)^2}}) + \Theta({\underline}{\epsilon}/l)\\
    &= (1 + \Theta(l^{-1/2}))(\sigma_v^2 (1 - \Theta(l^{-1/2})) + \sigma_a^2 -\sigma_v^2 \f 2 \pi\f{1}{\sqrt{1 - ({\mathbf{e}}^*)^2}})/(\sigma_v^2 + \sigma_a^2) + \Theta({\underline}{\epsilon}/ l)\\
    &= \delta^* + O(l^{-1/2}),\end{aligned}$$ where $\delta^* := 1 - \f 2 \pi \f 1 {\sqrt{1 - ({\mathbf{e}}^*)^2}} \f{\sigma_v^2 }{\sigma_v^2 + \sigma_a^2}$, which is positive by \[lemma:fixed\_point\_ineq\]. By taking the $\delta$ of \[lemma:alphaDeltaDynamics\] to be $\delta^* + \varepsilon$ or $\delta^* - \varepsilon$ respectively for lower and upper bounding the dynamics of $\p {\epsilon}l$, the solution $\p {\epsilon}l$ is $\Omega(l^{-\delta^* - \varepsilon})$ and $O(l^{-\delta^* + \varepsilon})$ for any $\varepsilon > 0$ since $\f 1 2 > \delta^*$.

### Backward Dynamics

The $\sigma_w = 0$ case is obvious. We will assume $\sigma_w > 0$ from here on.

As in the proof of \[thm:dalethExpSqrtTanh\], $$\begin{aligned}
    \log(\p {{\boldsymbol{\oldchi}}}{m} / \p {{\boldsymbol{\oldchi}}}l) &= 2B D \sigma_w^{-1} b_0^{-1/2} (\sqrt l - \sqrt m)
\\
&\phantom{={}} -(BD \sigma_w^{-1} b_0^{-3/2} b_1 2^{-1} +B^2 D^2 \sigma_w^{-2} b_0^{-1}2^{-1})(\log l - \log m) + O(1)\end{aligned}$$ where $B = \sigma_v^2 \sigma_w^2, D = \f 2 3 \sqrt{\f 2 \pi},$ $$\begin{aligned}
b_0 &= \sigma_v^2 + \sigma_a^2\\
b_1 &= \f{-2C \sigma_v^2 \sigma_w^{-1}}{\sqrt{\sigma_v^2 + \sigma_a^2}}\\
b_2 &= \f{-C^2 \sigma_v^4 \sigma_w^{-2}}{(\sigma_v^2 + \sigma_a^2)^2}.\end{aligned}$$ with $C = \sqrt{\f 2 \pi}$. This simplifies to the desired form.

Similar to \[thm:dalethExpSqrtTanhAllGrad\].

$\alpha$-ReLU: Full Residual Network
------------------------------------

The following can be checked readily

Since $\dot \psi_\alpha = \alpha \psi_{\alpha - 1}$, we have as a corollary,

\[lemma:Vt\_dotpsi\_alpha\] If $\alpha > \f 1 2$, then ${\mathrm{V}}\dot \psi_\alpha( q) = \alpha^2 {\mathsf{c}}_{\alpha - 1} q^{\alpha-1}$.

As a special case, when $\alpha = 1$, ${\mathsf{c}}_\alpha = \f 1 2$.

The following is a trivial computation, but useful for many simplifications.

${\mathsf{c}}_{\alpha+1}/{\mathsf{c}}_\alpha = 2\alpha+1$.

### Forward Dynamics {#sec:AlphaReluForwardProofs}

\[thm:pDynamic1ReLU\] Suppose we have the nonlinearity $\phi = \psi_1$. Then $\p {\mathbf{p}}l = \Theta((1 + \sigma_v^2 \sigma_w^2/2)^l)$, with the hidden constant depending on the initial condition.

We have $$\begin{aligned}
{\mathbf{p}}&= \f 1 2\sigma_v^2(\sigma_w^2 {\underline}{\mathbf{p}}+ \sigma_b^2) + \sigma_a^2 + {\underline}{\mathbf{p}}\\
    &= (\f 1 2\sigma_v^2\sigma_w^2 + 1) {\underline}{\mathbf{p}}+ \f 1 2 (\sigma_v^2 \sigma_b^2 + \sigma_a^2).\end{aligned}$$ By the standard method of characteristic equation, we get that $$\p {\mathbf{p}}l = A + C B^l$$ where $A = - \f{\sigma_a^2 + \sigma_b^2 \sigma_v^2}{\sigma_v^2 \sigma_w^2}$, $B = 1 + \f{\sigma_v^2 \sigma_w^2}{2}$, and $C$ is a coefficient determined by initial conditions.

\[thm:pDynamicLT1ReLU\] Suppose $\alpha < 1$. We have the following asymptotic expansion $$\p {\mathbf{p}}l = K_1 l^{{\f 1 {1-\alpha}}}+ R(l)$$ where the remainder term $$R(l) \sim \begin{cases}
    -K_2 l^{{{\f \alpha {1-\alpha}}}} \log l &\text{if $\alpha > \f 1 2$}\\
    (C-K_2) l \log l            &\text{if $\alpha = \f 1 2$ and $K_2 \not=C$}\\
    \f{C(1-\alpha)}{1-2\alpha} l    & \text{if $ \alpha < \f 1 2$}
    \end{cases}$$ where $K_1 = [\sigma_v^2 \sigma_w^{2\alpha} {\mathsf{c}}_\alpha (1 - \alpha)]^{\f 1 {1-\alpha}}, K_2 = \f 1 2 [\sigma_v^2 {\mathsf{c}}_\alpha \sigma_w^{2\alpha}]^{{{\f 1 {1-\alpha}}}} (1-\alpha)^{{{\f \alpha {1-\alpha}}}- 1} \alpha$ and $C = \sigma_a^2$.

\[fig:55reluverifyleadingcoeffp\] verifies the leading coefficient and the exponent of the leading term.

The difference equation governing the evolution of ${\mathbf{p}}$ is $${\mathbf{p}}- {\underline}{\mathbf{p}}= A({\underline}{\mathbf{p}}+ B)^\alpha + C$$ where $A =\sigma_v^2 {\mathsf{c}}_\alpha \sigma_w^{2\alpha}$, $B = \sigma_b^2/\sigma_w^2$, and $C = \sigma_a^2$. Then \[lemma:polyDynamicsConstant\] yields the result.

![Verification of leading term of \[thm:pDynamic1ReLU\] for $\alpha = 0.55$.[]{data-label="fig:55reluverifyleadingcoeffp"}](graphics/55relu_verify_leading_coeff_p){height=".2\textheight"}

\[thm:pDynamicLT1ReLU\] combined with \[thm:pDynamic1ReLU\] gives the following result.

By [@cho_kernel_2009], we know that ${\mathrm{W}}\psi_\alpha( q, q c) = {\mathrm{V}}\psi_\alpha( q) {\mathbb{J}}_\alpha(c)$, where ${\mathbb{J}}_\alpha(c) = J_\alpha(\arccos c)$ and $$J_\alpha(\theta) := \f 1 {2\pi {\mathsf{c}}_\alpha} (\sin \theta)^{2\alpha + 1} \Gamma(\alpha + 1)\int_0^{\pi / 2} \f{\dd \eta \cos^\alpha \eta}{(1 - \cos \theta \cos \eta)^{1 + \alpha}}.\tag{$\triangle$}\label{eqn:JalphaIntegralFormula}$$

Note that ${\mathbb{J}}_\alpha(c) \in (-\infty, \infty)$ for $\alpha \in (-1, \infty)$ and any $c \in (0, 1)$, even though ${\mathrm{V}}\psi_\alpha$ is only defined for $\alpha > -1/2$.

![(a) ${\mathbb{J}}_\alpha$ for different $\alpha$s and the identity function. From this plot, it looks like ${\mathbb{J}}_\alpha(c) \ge c$ and $\dot{\mathbb{J}}_\alpha(c) \le 1$ for all $\alpha \in (\f 1 2, 1]$ with equality iff $c = 1$, but this is misleading. (b) shows $|{\mathbb{J}}_\alpha(c) - c|$ in log scale. Where the curves dip below the x-axis indicate points where ${\mathbb{J}}_\alpha(c) = c$. We see that in fact every ${\mathbb{J}}_\alpha$ has a solution ${\mathbb{J}}_\alpha(c) = c$ for a $c < 1$, when $\alpha < 1$. (c) Furthermore, at each such $c$, $\dot{\mathbb{J}}_\alpha < 1$.[]{data-label="fig:jjj_vs_id"}](graphics/jjjVsId2.pdf "fig:"){width=".3\textwidth"}   ![(a) ${\mathbb{J}}_\alpha$ for different $\alpha$s and the identity function. From this plot, it looks like ${\mathbb{J}}_\alpha(c) \ge c$ and $\dot{\mathbb{J}}_\alpha(c) \le 1$ for all $\alpha \in (\f 1 2, 1]$ with equality iff $c = 1$, but this is misleading. (b) shows $|{\mathbb{J}}_\alpha(c) - c|$ in log scale. Where the curves dip below the x-axis indicate points where ${\mathbb{J}}_\alpha(c) = c$. We see that in fact every ${\mathbb{J}}_\alpha$ has a solution ${\mathbb{J}}_\alpha(c) = c$ for a $c < 1$, when $\alpha < 1$. (c) Furthermore, at each such $c$, $\dot{\mathbb{J}}_\alpha < 1$.[]{data-label="fig:jjj_vs_id"}](graphics/jjj_vs_id_log.pdf "fig:"){width=".3\textwidth"}   ![(a) ${\mathbb{J}}_\alpha$ for different $\alpha$s and the identity function. From this plot, it looks like ${\mathbb{J}}_\alpha(c) \ge c$ and $\dot{\mathbb{J}}_\alpha(c) \le 1$ for all $\alpha \in (\f 1 2, 1]$ with equality iff $c = 1$, but this is misleading. (b) shows $|{\mathbb{J}}_\alpha(c) - c|$ in log scale. Where the curves dip below the x-axis indicate points where ${\mathbb{J}}_\alpha(c) = c$. We see that in fact every ${\mathbb{J}}_\alpha$ has a solution ${\mathbb{J}}_\alpha(c) = c$ for a $c < 1$, when $\alpha < 1$. (c) Furthermore, at each such $c$, $\dot{\mathbb{J}}_\alpha < 1$.[]{data-label="fig:jjj_vs_id"}](graphics/jjj_vs_id_close.pdf "fig:"){width=".3\textwidth"}

\[fig:jjj\_vs\_id\] shows a comparison of ${\mathbb{J}}_\alpha$ for different $\alpha$s along with the identity function. By [@daniely_toward_2016 Lemma 11], ${\mathbb{J}}_\alpha$ is an increasing and convex function as long as $\psi_\alpha^2$ is Gaussian-integrable, which is precisely when $\alpha > -1/2$. We can compute ${\mathbb{J}}_\alpha(1) = {\mathrm{W}}\psi_\alpha(q, q)/{\mathrm{V}}\psi_\alpha(q) = 1$, and ${\mathbb{J}}_\alpha(0) = {\mathrm{W}}\psi_\alpha(q, 0)/{\mathrm{V}}\psi_\alpha(q) = {\mathrm{V}}\psi_{\alpha/2}(q)^2/{\mathrm{V}}\psi_\alpha(q) = {\mathsf{c}}_{\alpha/2}^2/{\mathsf{c}}_\alpha = \f 1{2\sqrt \pi} \f{\Gamma(\f \alpha 2 + \f 1 2)^2}{\Gamma(\alpha + \f 1 2)}$. We record these observations as a lemma.

\[lemma:basicJalpha\] ${\mathbb{J}}_\alpha(c)$ is an increasing and convex function for each $\alpha > -1/2$ on $c \in [0, 1]$. ${\mathbb{J}}_\alpha(1) = 1$ and ${\mathbb{J}}_\alpha(0) = \f 1{2\sqrt \pi} \f{\Gamma(\f \alpha 2 + \f 1 2)^2}{\Gamma(\alpha + \f 1 2)}$.

![${\mathbb{J}}_1$ vs identity[]{data-label="fig:jjj1_vs_id"}](graphics/jjj1_vs_id.pdf){width=".3\textwidth"}

For $\alpha = 1$, @cho_kernel_2009 computed $${\mathbb{J}}_1(c) = \f 1 \pi (\sqrt{1 - c^2} + (\pi - \arccos(c))c).$$ \[fig:jjj1\_vs\_id\] shows a plot of ${\mathbb{J}}_1$ vs identity. It has derivative $\dot{\mathbb{J}}_1(c) = 1 - \f 1 \pi \arccos c$, which shows that $\dot{\mathbb{J}}_1(c) < 1$ with equality iff $c = 1$, and consequently ${\mathbb{J}}_1(c) \ge c$ with equality iff $c = 1$. At the same time, $\dot {\mathbb{J}}_1(c) \ge 0$ with equality iff $c = -1$, so ${\mathbb{J}}_1$ is increasing on $[-1, 1]$. It has an asymptotic expansion ${\mathbb{J}}_1(1 - \varepsilon) = 1 - \varepsilon + \f {2\sqrt 2}{3\pi} {\epsilon}^{3/2} + \Theta({\epsilon}^{5/2})$ at 1.

The zeroth Bessel function of the second kind is defined by ${{\mathcal{K}_0}}(z) = \int_1^\infty e^{-z x} (x^2 - 1)^{-1/2} \dd x$. It is one of the fundamental solutions to the homogeneous differential equation $x^2 \dot y + x \dot y - x^2 y = 0$. The following lemma shows that $J_\alpha$ can be expressed in terms of ${{\mathcal{K}_0}}$.

\[lemma:JalphaBessel\] For any $\alpha > -1$, $J_\alpha(\theta) = \f 1 {2\pi {\mathsf{c}}_\alpha}\sin^{2\alpha+1} \theta \int_0^\infty \dd x {{\mathcal{K}_0}}(x) e^{x\cos \theta} x^\alpha$

@cho_kernel_2009 gave the expression $$2\pi {\mathsf{c}}_\alpha J_\alpha(\theta) = \csc \theta \int_0^\infty \dd u \int_0^\infty \dd v e^{-(u^2+v^2 - 2uv \cos \theta)/2\sin^2\theta} u^\alpha v^\alpha.$$ Note that the integrand is symmetric in $u$ and $v$.

Thus, if ${{\mathsf V}}= \{(u, v): u, v \ge 0 \And v \ge u\}$, then $$2\pi {\mathsf{c}}_\alpha J_\alpha(\theta) = 2\csc \theta \int_{{\mathsf V}}\dd u \dd v e^{-(u^2+v^2 - 2uv \cos \theta)/2\sin^2\theta} u^\alpha v^\alpha.$$

Now make the change of variables from ${{\mathsf V}}$ to $\{({\mathbbm{p}}, {\mathbbm{q}}): {\mathbbm{q}}\ge 2 \sqrt{{\mathbbm{p}}}\}$: $$\begin{aligned}
    {\mathbbm{p}}&= uv   &   \dd {\mathbbm{p}}&= v \dd u + u \dd v\\
    {\mathbbm{q}}&= u + v&   \dd {\mathbbm{q}}&= \dd u + \dd v\\
    \dd {\mathbbm{p}}\dd {\mathbbm{q}}&= (v - u) \dd u \dd v  &   \dd u \dd v &= ({\mathbbm{q}}^2 - 4 {\mathbbm{p}})^{-1/2} \dd {\mathbbm{p}}\dd {\mathbbm{q}}\end{aligned}$$ so that we have $$2 \pi {\mathsf{c}}_\alpha J_\alpha(\theta) = 2\csc \theta \int_0^\infty \dd {\mathbbm{p}}e^{{\mathbbm{p}}(1 + \cos \theta)\csc^2\theta} {\mathbbm{p}}^\alpha \int_{2\sqrt{{\mathbbm{p}}}}^\infty \dd {\mathbbm{q}}e^{-{\mathbbm{q}}^2\csc^2 \theta}({\mathbbm{q}}^2 - 4 {\mathbbm{p}})^{-1/2}.$$ The inner integral in ${\mathbbm{q}}$ can be expressed in terms of ${{\mathcal{K}_0}}$ by a change of variable $x = {\mathbbm{q}}^2/2\sqrt {\mathbbm{p}}$: $$\begin{aligned}
    2 \pi {\mathsf{c}}_\alpha J_\alpha(\theta) &= 2\csc \theta \int_0^\infty \dd {\mathbbm{p}}e^{{\mathbbm{p}}(1 + \cos \theta)\csc^2\theta} {\mathbbm{p}}^\alpha \f 1 2 e^{-{\mathbbm{p}}\csc^2 \theta} {{\mathcal{K}_0}}({\mathbbm{p}}\csc^2 \theta)\\
        &= \csc\theta \int_0^\infty \dd {\mathbbm{p}}{{\mathcal{K}_0}}({\mathbbm{p}}\csc^2\theta) e^{{\mathbbm{p}}\cos \theta \csc^2\theta} {\mathbbm{p}}^\alpha\\
        &= \sin^{2\alpha + 1} \theta \int_0^\infty \dd x {{\mathcal{K}_0}}(x) e^{x \cos \theta} x^\alpha
    \end{aligned}$$

Define $L_\alpha(\theta) = 2 \pi {\mathsf{c}}_\alpha J_\alpha(\theta)\csc^{2\alpha+1}\theta = \int_0^\infty \dd x {{\mathcal{K}_0}}(x) e^{x \cos \theta} x^\alpha$.

\[lemma:LAlphaRec\] If $\alpha > 1$, then $$L_\alpha(\theta) = \csc^2 \theta [(2\alpha-1) \cos \theta L_{\alpha-1}(\theta) + (\alpha-1)^2 L_{\alpha - 2}(\theta)].$$

We will prove this claim for $\theta < 1$, and by continuity this also proves the case $\theta = 1$. As remarked above, ${{\mathcal{K}_0}}(z) = \ddot {{\mathcal{K}_0}}(z) + \inv z \dot {{\mathcal{K}_0}}(z).$ Thus $$\begin{aligned}
    L_\alpha(\theta) &= \int_0^\infty \dd x (\ddot {{\mathcal{K}_0}}(x) + \inv x \dot {{\mathcal{K}_0}}(x)) e^{x \cos\theta} x^\alpha\\
        &=\dot {{\mathcal{K}_0}}e^{x\cos\theta} x^\alpha \rvert_0^\infty + {{\mathcal{K}_0}}e^{x\cos \theta} x^{\alpha - 1} \rvert_0^\infty\\
            &\phantom{={}} - \int \dd x[\cos \theta e^{x \cos \theta} x^\alpha + \alpha e^{x \cos \theta} x^{\alpha - 1}] \dot {{\mathcal{K}_0}}\\
            &\phantom{={}} - \int \dd x[\cos \theta e^{x \cos \theta} x^{\alpha - 1} + (\alpha-1) e^{x\cos\theta} x^{\alpha - 2}] {{\mathcal{K}_0}}\end{aligned}$$ Asymptotically, ${{\mathcal{K}_0}}(z) \sim \sqrt{\f{\pi}{2z}} e^{-z}$ as $z \to \infty$ and ${{\mathcal{K}_0}}(z) \sim -\ln(z)$ as $z \searrow 0$, and $\dot{{\mathcal{K}_0}}(z) \sim -\sqrt{\f{\pi}{2z}} e^{-z}$ as $z \to \infty$ and $\dot{{\mathcal{K}_0}}(z) \sim -\inv z$ as $z \searrow 0$. Thus, as $\alpha > 1$, $$\begin{aligned}
    \dot {{\mathcal{K}_0}}e^{x\cos\theta} x^\alpha \rvert_0^\infty = -\lim_{x \to \infty} \sqrt{\pi/2} e^{-x(1 - \cos \theta)}x^{\alpha - 1} + \lim_{x\searrow 0} e^{x\cos\theta} x^{\alpha-1}  = 0\\
    {{\mathcal{K}_0}}e^{x\cos\theta} x^{\alpha-1} \rvert_0^\infty = -\lim_{x \to \infty} \sqrt{\pi/2} e^{-x(1 - \cos \theta)}x^{\alpha - 2} + \lim_{x\searrow 0} e^{x\cos\theta} x^{\alpha-1} \ln x  = 0
    \end{aligned}$$ So $$\begin{aligned}
    L_\alpha(\theta) &= -\cos \theta L_{\alpha-1}(\theta) - (\alpha-1)L_{\alpha-2}(\theta) - \int \dd x[\cos \theta e^{x \cos \theta} x^\alpha + \alpha e^{x \cos \theta} x^{\alpha - 1}] \dot {{\mathcal{K}_0}}\end{aligned}$$ Via another integration by parts, the integral on the right is $$\begin{aligned}
    &\phantom{={}} \cos \theta e^{x \cos \theta} x^\alpha {{\mathcal{K}_0}}\rvert_0^\infty + \alpha e^{x\cos \theta} x^{\alpha-1} {{\mathcal{K}_0}}\rvert_0^\infty\\
    &\phantom{={}} - \int \dd x[\cos^2 \theta e^{x \cos \theta} x^\alpha + 2 \alpha \cos \theta e^{x \cos \theta} x^{\alpha - 1} + \alpha (\alpha - 1) e^{x \cos \theta} x^{\alpha - 2}] {{\mathcal{K}_0}}\\
    &= - [\cos^2 \theta L_\alpha(\theta) + 2 \alpha \cos \theta L_{\alpha-1}(\theta) + \alpha(\alpha-1) L_{\alpha-2}(\theta)]
    \end{aligned}$$ where the evaluation terms vanish just like before. Altogether, we have $$\begin{aligned}
    L_\alpha(\theta) &= \cos^2 \theta L_\alpha(\theta) + (2\alpha-1)\cos\theta L_{\alpha-1}(\theta) + (\alpha-1)^2 L_{\alpha-2}(\theta)\\
                 &= \csc^2 \theta[(2\alpha-1) \cos \theta L_{\alpha-1}(\theta) + (\alpha-1)^2 L_{\alpha-2}(\theta)]
    \end{aligned}$$

As a corollary we get

\[lemma:JalphaRec\] Suppose $\alpha > 1$. Then $$\begin{aligned}
    J_\alpha(\theta) &= \cos \theta J_{\alpha-1}(\theta) + (\alpha-1)^2 (2\alpha-1)^{-1}(2\alpha-3)^{-1} \sin^2 \theta J_{\alpha-2}(\theta)\\
    {\mathbb{J}}_\alpha(c) &=  c {\mathbb{J}}_{\alpha-1}(c) + (\alpha-1)^2 (2\alpha-1)^{-1}(2\alpha-3)^{-1}(1-c^2) {\mathbb{J}}_{\alpha-2}(c)
    \end{aligned}$$

The derivative of $J_\alpha(\theta)$ turns out to be quite simple.

\[lemma:JalphaGrad\] Suppose $\alpha > 0$. Then $$\begin{aligned}
    \dot J_\alpha(\theta) &= -\alpha^2(2 \alpha-1)^{-1} J_{\alpha-1}(\theta) \sin \theta\\
    \dot {\mathbb{J}}_\alpha(c) &= \alpha^2 (2 \alpha-1)^{-1} {\mathbb{J}}_{\alpha-1}(c)
    \end{aligned}$$

We will prove the first formula. The second follows from chain rule. By \[lemma:JalphaBessel\], $$\begin{aligned}
    J_\alpha(\theta) &= \f 1 {2\pi {\mathsf{c}}_\alpha}\sin^{2\alpha+1} \theta \int \dd x {{\mathcal{K}_0}}(x) e^{x\cos \theta} x^\alpha\\
    \dot J_\alpha(\theta) &= \f 1 {2\pi {\mathsf{c}}_\alpha}[(2\alpha+1)\sin^{2\alpha} \theta \cos\theta \int \dd x {{\mathcal{K}_0}}(x) e^{x\cos \theta} x^\alpha\\
                        &\phantom{={}} -\sin^{2\alpha+2}\theta \int \dd x {{\mathcal{K}_0}}(x) e^{x\cos \theta} x^{\alpha+1}]\\
                &= (2\alpha+1) \cot \theta J_\alpha(\theta) - \f{c_{\alpha+1}}{c_\alpha}\csc\theta J_{\alpha+1}(\theta)\\
                &= (2\alpha+1) \csc \theta [\cos \theta J_\alpha(\theta) - J_{\alpha+1}(\theta)].\\
    \end{aligned}$$ As $\alpha + 1 > 1$, by \[lemma:JalphaRec\], this is $$\begin{aligned}
    &\phantom{={}}-(2\alpha+1)\csc \theta [(\alpha-1)^2 (2\alpha+1)^{-1}(2\alpha-1)^{-1} \sin^2 \theta J_{\alpha-1}(\theta)]\\
    &= -(\alpha-1)^2(2\alpha-1)^{-1} \sin \theta J_{\alpha-1}(\theta).
    \end{aligned}$$

Thus $\dot {\mathbb{J}}_\alpha(1) = \alpha^2(2\alpha-1)^{-1} {\mathbb{J}}_{\alpha-1}(1) = \alpha^2(2\alpha-1)^{-1}$ for any $\alpha > 0$ by \[lemma:basicJalpha\]. For $1/2 < \alpha \le 1$, $\dot {\mathbb{J}}_\alpha(1) \ge 1$ with equality iff $\alpha = 1$, and for $\alpha = 1/2$, $\dot {\mathbb{J}}_\alpha(1) =\infty > 1$ by continuity of $\dot {\mathbb{J}}_\alpha(c)$ in $\alpha$. Because for $\alpha > -1/2$, ${\mathbb{J}}_\alpha$ is increasing and convex on $[0, 1]$ and ${\mathbb{J}}_\alpha(0) > 0$ by \[lemma:basicJalpha\], ${\mathbb{J}}_\alpha$ intersects identity at a unique point away from 1 when $\alpha \in [1/2, 1)$. We record this as a theorem.

\[thm:stableFixedPointsJJ\] For $\alpha \in [1/2, 1)$, ${\mathbb{J}}_\alpha(c) = c$ has two solutions: an unstable solution at 1 (“unstable” meaning $\dot{\mathbb{J}}_\alpha(1) > 1$) and a stable solution in ${\mathbf{e}}^* \in (0, 1)$ (“stable” meaning $\dot {\mathbb{J}}_\alpha({\mathbf{e}}^*) < 1$).

![Left-to-right: **(a)** ${\mathbb{J}}_\alpha$ for different $\alpha$s and the identity function (black, dashed line). ${\mathbb{J}}_1$ is highlighted in red. From this plot, it looks like ${\mathbb{J}}_\alpha(c) \ge c$ and $\dot{\mathbb{J}}_\alpha(c) \le 1$ for all $\alpha \in (\f 1 2, 1]$ with equality iff $c = 1$, but this is misleading. **(b)** shows $|{\mathbb{J}}_\alpha(c) - c|$ in log scale. Where the curves dip below the x-axis indicate points where ${\mathbb{J}}_\alpha(c) = c$. We see that in fact every ${\mathbb{J}}_\alpha$ has a solution ${\mathbb{J}}_\alpha(c) = c$ for a $c < 1$, when $\alpha < 1$. **(c)** Furthermore, at each such $c$, $\dot{\mathbb{J}}_\alpha < 1$. (b) and (c) demonstrate the existence of stable fixed points away from 1 for ${\mathbb{J}}_\alpha, \alpha \in (1/2, 1)$, which is confirmed rigorously by \[thm:stableFixedPointsJJ\]. []{data-label="fig:jjj_vs_id_main"}](graphics/jjjVsId2.pdf "fig:"){height=".1\textheight"}   ![Left-to-right: **(a)** ${\mathbb{J}}_\alpha$ for different $\alpha$s and the identity function (black, dashed line). ${\mathbb{J}}_1$ is highlighted in red. From this plot, it looks like ${\mathbb{J}}_\alpha(c) \ge c$ and $\dot{\mathbb{J}}_\alpha(c) \le 1$ for all $\alpha \in (\f 1 2, 1]$ with equality iff $c = 1$, but this is misleading. **(b)** shows $|{\mathbb{J}}_\alpha(c) - c|$ in log scale. Where the curves dip below the x-axis indicate points where ${\mathbb{J}}_\alpha(c) = c$. We see that in fact every ${\mathbb{J}}_\alpha$ has a solution ${\mathbb{J}}_\alpha(c) = c$ for a $c < 1$, when $\alpha < 1$. **(c)** Furthermore, at each such $c$, $\dot{\mathbb{J}}_\alpha < 1$. (b) and (c) demonstrate the existence of stable fixed points away from 1 for ${\mathbb{J}}_\alpha, \alpha \in (1/2, 1)$, which is confirmed rigorously by \[thm:stableFixedPointsJJ\]. []{data-label="fig:jjj_vs_id_main"}](graphics/jjj_vs_id_log.pdf "fig:"){height=".1\textheight"}   ![Left-to-right: **(a)** ${\mathbb{J}}_\alpha$ for different $\alpha$s and the identity function (black, dashed line). ${\mathbb{J}}_1$ is highlighted in red. From this plot, it looks like ${\mathbb{J}}_\alpha(c) \ge c$ and $\dot{\mathbb{J}}_\alpha(c) \le 1$ for all $\alpha \in (\f 1 2, 1]$ with equality iff $c = 1$, but this is misleading. **(b)** shows $|{\mathbb{J}}_\alpha(c) - c|$ in log scale. Where the curves dip below the x-axis indicate points where ${\mathbb{J}}_\alpha(c) = c$. We see that in fact every ${\mathbb{J}}_\alpha$ has a solution ${\mathbb{J}}_\alpha(c) = c$ for a $c < 1$, when $\alpha < 1$. **(c)** Furthermore, at each such $c$, $\dot{\mathbb{J}}_\alpha < 1$. (b) and (c) demonstrate the existence of stable fixed points away from 1 for ${\mathbb{J}}_\alpha, \alpha \in (1/2, 1)$, which is confirmed rigorously by \[thm:stableFixedPointsJJ\]. []{data-label="fig:jjj_vs_id_main"}](graphics/jjj_vs_id_close.pdf "fig:"){height=".1\textheight"}

This result confirms that pictures presented in \[fig:jjj\_vs\_id\_main\]b,c are qualitatively correct, that there are indeed stable fixed points of ${\mathbb{J}}_\alpha$ away from 1.

If ${\underline}{\mathbf{e}}< 1$, then $$\begin{aligned}
{c}= \f{\sigma_w^2 {\underline}{\boldsymbol{\oldgamma}}+ \sigma_b^2}{\sigma_w^2 {\underline}{\mathbbm{p}}+ \sigma_b^2} &\ge {\underline}{\mathbf{e}}\\
{\mathbb{J}}_1({c}) &\ge {\mathbb{J}}_1({\underline}{\mathbf{e}})\\
{\mathbf{e}}= \f{\sigma_v^2 {\mathsf{c}}_\alpha {\mathbbm{q}}^\alpha{\mathbb{J}}_1({c}) + \sigma_b^2}{\sigma_v^2 {\mathsf{c}}_\alpha {\mathbbm{q}}^\alpha + \sigma_b^2} &\ge {\mathbb{J}}_1({\underline}{\mathbf{e}}) \end{aligned}$$ but ${\mathbf{e}}\ge {\mathbb{J}}_1({\underline}{\mathbf{e}}) > {\underline}{\mathbf{e}}$ as noted above. Thus by monotone convergence ${\mathbf{e}}$ converges, and ${\mathbf{e}}^* = 1$ is the only possible fixed point.

By \[lemma:cExpansion\], ${c}= {\underline}{\mathbf{e}}(1 + \Theta({\underline}{\epsilon}{\underline}{\mathbbm{p}}^{-1})) = 1 - {\underline}{\epsilon}+ \Theta({\underline}{\epsilon}\inv {\mathbbm{p}}) = 1 - u{\underline}{\epsilon}$ where $u := 1 - \Theta({\underline}{\mathbbm{p}}^{-1})$. Using the asymptotic expansion ${\mathbb{J}}_1(1 - {\epsilon}) = 1 - {\epsilon}+ U {\epsilon}^{3/2} + \Theta({\epsilon}^{5/2})$, we have $$\begin{aligned}
(1 - {\epsilon}){\mathbbm{p}}&= \sigma_v^2 \f {\mathbbm{q}}2 {\mathbb{J}}_1(1 - u{\underline}{\epsilon}) + \sigma_a^2 + (1 - {\underline}{\epsilon}) {\underline}{\mathbbm{p}}\\
-{\epsilon}{\mathbbm{p}}&= \sigma_v^2 \f {\mathbbm{q}}2 ({\mathbb{J}}_1(1 - u{\underline}{\epsilon}) - 1) - {\underline}{\epsilon}{\underline}{\mathbbm{p}}\\
    &= \sigma_v^2 \f {\mathbbm{q}}2 [- u{\underline}{\epsilon}+ U u^{3/2}{\underline}{\epsilon}^{3/2} + \Theta(u^{5/2}{\underline}{\epsilon}^{5/2})] - {\underline}{\epsilon}{\underline}{\mathbbm{p}}\\
{\epsilon}&= {\underline}{\epsilon}\f 1 {\mathbbm{p}}[{\underline}{\mathbbm{p}}+ \sigma_v^2 \f {\mathbbm{q}}2 (u - U u^{3/2} {\underline}{\epsilon}^{1/2} + \Theta(u^{5/2}{\underline}{\epsilon}^{3/2}))]\\
    &=  {\underline}{\epsilon}\f 1 {\mathbbm{p}}[{\mathbbm{p}}- \sigma_a^2 + \sigma_v^2 \f {\mathbbm{q}}2 (\Theta(\inv {{\underline}{\mathbbm{p}}}) - U u^{3/2} {\underline}{\epsilon}^{1/2} + \Theta(u^{5/2}{\underline}{\epsilon}^{3/2}))]\\
    &=  {\underline}{\epsilon}[1 + \f{-\sigma_a^2 + \sigma_v^2 \f {\mathbbm{q}}2 (\Theta(\inv {{\underline}{\mathbbm{p}}}) - U u^{3/2} {\underline}{\epsilon}^{1/2} + \Theta(u^{5/2}{\underline}{\epsilon}^{3/2}))}{{\mathbbm{p}}}]\\
    &=  {\underline}{\epsilon}[1 + \f{-\sigma_a^2\inv {\mathbbm{q}}+ \f 1 2 \sigma_v^2 (\Theta(\inv {{\underline}{\mathbbm{p}}}) - U u^{3/2} {\underline}{\epsilon}^{1/2} + \Theta(u^{5/2}{\underline}{\epsilon}^{3/2}))}{{\mathbbm{p}}\inv {\mathbbm{q}}}]
$$ Let the content of the bracket on the RHS be $\aleph$. We have ${\mathbbm{p}}\inv {\mathbbm{q}}= (1+ o(1))B/\sigma_w^2$. If ${\epsilon}= O({\underline}{\mathbbm{p}}^{-1})$, then $\aleph = 1 - O(\inv {\mathbbm{p}})$, but because ${\mathbbm{p}}$ is exponentially decreasing, this means ${\epsilon}= \Theta(1)$ and does not converge to 0 — this is a contradiction. Therefore, ${\underline}{\epsilon}= \omega({\underline}{\mathbbm{p}}^{-1})$, and $$\begin{aligned}
{\epsilon}&= {\underline}{\epsilon}[1 - \f 1 2 B^{-1}\sigma_v^2 \sigma_w^2 U{\underline}{\epsilon}^{1/2}(1 + o(1))]\\
{\epsilon}- {\underline}{\epsilon}&= - \f 1 2 B^{-1}\sigma_v^2 \sigma_w^2 U{\underline}{\epsilon}^{3/2}(1 + o(1))\end{aligned}$$ Using \[lemma:polyDynamicsPosAlpha\] to upper and lower bound our dynamics, we get that $\p {\epsilon}l \sim [\f 1 4 \sigma_v^2 \sigma_w^2 \inv B U l]^{-2}$.

\[lemma:separableDynamics\] Let $\phi$ be any nonlinearity. Suppose ${\mathrm{W}}\phi( r, rd) = {\mathrm{V}}\phi( r){\mathbb{K}}(d)$ for some twice differentiable function ${\mathbb{K}}(d)$ independent of ${\mathbbm{q}}$, where ${\mathbb{K}}(1) = 1$ naturally. Suppose further that

-   ${\mathbb{K}}(d) = d$ has a solution $d = {\mathbf{e}}^* > 0$ where $\dot{\mathbb{K}}({\mathbf{e}}^*) = \delta < 1$;

-   ${\mathbb{K}}(d) > d$ for all $d < {\mathbf{e}}^*$ and ${\mathbb{K}}(d) < d$ for all $1 > d > {\mathbf{e}}^*$; and

-   ${\mathbb{K}}$ is nondecreasing.

Let $\p {\epsilon}l := \p {\mathbf{e}}l - {\mathbf{e}}^*$ and suppose $\p {\mathbf{e}}0 < 1$. If $\p {\boldsymbol{\oldgamma}}l \to \infty$ and ${\mathrm{V}}\phi( \p {\mathbbm{q}}l) \to \infty$, then $\p {\epsilon}l \to 0$ and satisfies $${\epsilon}= {\underline}{\epsilon}\left(1 - \f{\sigma_a^2 + (1 - \delta + O({\underline}{\epsilon}))\sigma_v^2 {\mathrm{V}}\phi( {\mathbbm{q}})}{{\mathbbm{p}}}\right) + {\mathrm{V}}\phi( {\mathbbm{q}})\Theta(\inv {\boldsymbol{\oldgamma}}\inv {\mathbbm{p}}).$$

First we note that because ${\mathbf{e}}^*$ is the only stable fixed point of the dynamics $x \mapsto {\mathbb{K}}(x)$, with the basin of attraction $[0, 1)$, we can show $\p {\mathbf{e}}l \to {\mathbf{e}}^*$ as in the proof of \[thm:eDynamicsFullResTanh\] (using \[lemma:timeDependentConvergence\]).

Write $\p V l := {\mathrm{V}}\phi( \p {\mathbbm{q}}l)$. We first show that $\p {\mathbf{e}}l \to {\mathbf{e}}^*$. When $l$ is large, $$\begin{aligned}
{c}= \f{\sigma_w^2 {\underline}{\boldsymbol{\oldgamma}}+ \sigma_b^2}{\sigma_w^2 {\underline}{\mathbbm{p}}+ \sigma_b^2} & = {\underline}{\mathbf{e}}( 1 + O(\inv {\boldsymbol{\oldgamma}}))\\
{\mathbf{e}}= \f{\sigma_v^2V{\mathbb{K}}({c}) + \sigma_a^2}{\sigma_v^2 V + \sigma_a^2} &= {\mathbb{K}}({c})(1 + O(\inv V {\mathbb{K}}({c})^{-1})).\end{aligned}$$ If $\p {\boldsymbol{\oldgamma}}l$ is bounded for all $l$, then ${\mathbf{e}}\to 0$ because $\p {\mathbbm{p}}l \to \infty$. Since ${\mathbb{K}}({c}) > 0$ for ${c}\in [0, 1]$ and $\p V l\to \infty$, we have that in the limit $l\to \infty$, $\lim_{l \to \infty} {\mathbf{e}}= 0 = {\mathbb{K}}(\lim_{l \to \infty}{\mathbf{e}}) = {\mathbb{K}}(0)$ (by the continuity of ${\mathbb{K}}$), which is impossible by our assumptions. Thus $\p {\boldsymbol{\oldgamma}}l \to \infty$, and we have $\lim_{l \to \infty} {\mathbf{e}}= {\mathbb{K}}(\lim_{l \to \infty} {\mathbf{e}})$. By our assumptions, ${\mathbf{e}}^*$ is the only stable fixed point of ${\mathbb{K}}$ with basin of attraction $[0, 1)$, so this shows that ${\mathbf{e}}\to {\mathbf{e}}^*$ as desired.

Now we derive the equation in question. Note that ${c}= {\underline}{\mathbf{e}}(1 + \Theta(\inv {\boldsymbol{\oldgamma}}))$ because ${\mathbf{e}}^* < 1$. We use the Taylor expansion ${\mathbb{K}}({\mathbf{e}}^* + {\epsilon}) = {\mathbf{e}}^* + \delta {\epsilon}+ O({\epsilon}^2)$.

$$\begin{aligned}
({\mathbf{e}}^* + {\epsilon}){\mathbbm{p}}&= {\sigma_v^2}V {\mathbb{K}}\lp({\mathbf{e}}^* + {\underline}{\epsilon})(1 + \Theta(\inv {\boldsymbol{\oldgamma}}))\rp + {\sigma_a^2}+ ({\mathbf{e}}^* + {\underline}{\epsilon}) {\underline}{\mathbbm{p}}\\
    &= {\sigma_v^2}V({\mathbf{e}}^* + \delta({\underline}{\epsilon}+ \Theta(\inv {\boldsymbol{\oldgamma}})) + O({\underline}{\epsilon}^2)) + {\sigma_a^2}+ ({\mathbf{e}}^* + {\underline}{\epsilon}) {\underline}{\mathbbm{p}}\\
{\epsilon}{\mathbbm{p}}&= {\sigma_v^2}V(\delta({\underline}{\epsilon}+ \Theta(\inv {\boldsymbol{\oldgamma}})) + O({\underline}{\epsilon}^2)) + {\underline}{\epsilon}{\underline}{\mathbbm{p}}\\
{\epsilon}&= {\underline}{\epsilon}(1 - \f{{\sigma_a^2}+ (1 - \delta + O({\underline}{\epsilon})){\sigma_v^2}V}{\mathbbm{p}}) + \Theta(V \inv{\boldsymbol{\oldgamma}}\inv {\mathbbm{p}})\end{aligned}$$

We apply \[lemma:separableDynamics\]. We first check the conditions of the lemma, with ${\mathbb{K}}= {\mathbb{J}}_\alpha$. The following conditions were already verified.

-   ${\mathbb{J}}_\alpha$ has a fixed point ${\mathbf{e}}^*$ less than but very close to 1, where its slope is $\upsilon := \dot{\mathbb{J}}_\alpha({\mathbf{e}}^*) < 1$. (\[thm:stableFixedPointsJJ\])

-   ${\mathbb{J}}_\alpha(d) > d$ for all $d < {\mathbf{e}}^*$ and ${\mathbb{J}}_\alpha(d) < d$ for all $d > {\mathbf{e}}^*$. (By the convexity shown in \[lemma:basicJalpha\])

-   ${\mathbb{J}}_\alpha$ is nondecreasing (\[lemma:basicJalpha\]). Furthermore, from its integral formula (\[eqn:JalphaIntegralFormula\]), we see easily that ${\mathbb{J}}_\alpha$ is smooth at ${\mathbf{e}}^* < 1$.

We also proved the following

-   $\p {\mathbbm{p}}l \sim [\sigma_v^2 \sigma_w^{2\alpha} {\mathsf{c}}_\alpha (1 - \alpha)]^{\f 1 {1-\alpha}}l^{\f 1 {1-\alpha}}$ (\[thm:pDynamicLT1ReLU\]) and $\p {\boldsymbol{\oldgamma}}l$ is asymptotically a constant fraction of $\p {\mathbbm{p}}l$ (\[lemma:separableDynamics\]), so both go to $\infty$.

-   ${\mathrm{V}}\psi_\alpha( {\mathbbm{q}}) = {\mathsf{c}}_\alpha {\mathbbm{q}}^\alpha = {\mathsf{c}}_\alpha (\sigma_w^2 {\mathbbm{p}}+ \sigma_b^2)^\alpha = \Theta(l^{\alpha/(1-\alpha)})$, so goes to $\infty$. (\[lemma:VtPsiAlpha\])

Thus, for $\upsilon = \dot {\mathbb{J}}({\mathbf{e}}^*)$, $$\begin{aligned}
\f{\sigma_a^2 + (1 - \upsilon + O({\underline}{\epsilon}))\sigma_v^2 {\mathrm{V}}\phi( {\mathbbm{q}})}{{\mathbbm{p}}} &\sim \f{(1 - \upsilon)\sigma_v^2\sigma_w^{2\alpha}{\mathsf{c}}_\alpha}{{\mathbbm{p}}^{1-\alpha}}\\
    &= \inv l (1-\upsilon)/(1 - \alpha).\end{aligned}$$

Now, ${\mathrm{V}}\phi( {\mathbbm{q}})\inv {\boldsymbol{\oldgamma}}\inv {\mathbbm{p}}= \Theta(l^{-\f1{1-\alpha} - 1})$. By using the dynamics of \[lemma:alphaDeltaDynamics\] to upper and lower bound our dynamics, we have $\p {\epsilon}l = \Omega(l^{-\mu - {\epsilon}}), O(l^{-\mu + {\epsilon}})$ for any ${\epsilon}> 0$, where $\mu = \min((1-\upsilon)/(1-\alpha), 1/(1-\alpha)) = (1-\upsilon)/(1-\alpha).$

### Backward Dynamics

\[lemma:infVarAlphaLe75\] Suppose random variable $X \sim {\mathcal{N}}(0, \sigma^2)$, and $Y = \psi_{-\beta}(X)$ for some $\beta > 0 $, where $\psi_\alpha$ is $\alpha$-ReLU. Then for $\xi > 0$, $Y$ has density $$\Pr[Y \in [\xi, \xi+\dd \xi]] = \f 1 {\beta\sqrt{2\pi \sigma^2}} \xi^{-\f 1 \beta - 1} e^{-\xi^{-2/\beta}/2\sigma^2}.$$ At $\xi = 0$, $Y$ has density given by a Dirac delta of mass $\f 1 2$.

Furthermore, $Y$ has finite second moment iff $\beta < \f 1 2$.

We have $$\begin{aligned}
    \Pr[Y \in [\xi, \infty)] &= \Pr[X \in [0, \xi^{-1/\beta}]]\\
    &= \f 1 {\sqrt{2\pi \sigma^2}}\int_0^{\xi^{-1/\beta}} e^{-x^2/2\sigma^2}\dd x.
    \end{aligned}$$ Differentiating the RHS against $\xi$ using Leibniz’s rule, we get $$\begin{aligned}
    d\Pr[Y \in [\xi, \infty)]/d\xi &= \f 1 {\sqrt{2\pi \sigma^2}}  e^{-\xi^{-2/\beta}/2\sigma^2} \f d {d\xi}\xi^{-1/\beta}\\
    &= \f {-1} {\beta\sqrt{2\pi \sigma^2}} \xi^{-\f 1 \beta - 1} e^{-\xi^{-2/\beta}/2\sigma^2}.
    \end{aligned}$$ Negating both sides gives the density $f_Y$ of $Y$ for $\xi > 0$. For $\xi = 0$, observe that $\lim_{\xi \to 0} f_Y(\xi) = 0$ because, while $\xi^{-\f 1 \beta -1}$ blows up polynomially, $e^{-\xi^{-2/\beta}/2\sigma^2}$ blows up exponentially. Thus the contribution of $Y$’s mass at $Y = 0$ from $X > 0$ is 0. On the other hand, all $X < 0$ gets mapped to $Y = 0$, so $f_Y(0) = \f 1 2 \delta_0$, where $\delta_0$ is the Dirac delta.

For the second assertion, observe that $$\begin{aligned}
    f_Y(\xi) \sim \f 1 {\beta \sqrt{2\pi \sigma^2}} \xi^{-\f 1 \beta - 1} & \text{as ${\xi \to \infty}$}.
    \end{aligned}$$ Thus, $\xi^2 f_Y(\xi)$ is integrable iff $2 -\f 1 \beta - 1 < -1 \iff \beta < \f 1 2$.

Note that $\dot \psi_\alpha \propto \psi_{\alpha - 1}$, so it suffices to show that $\Var(\psi_{\alpha-1}(\zeta)^2) = \Var(\psi_{2\alpha - 2}(\zeta))$ is infinite for $\zeta \sim {\mathcal{N}}(0, \sigma^2)$. By \[lemma:infVarAlphaLe75\] with $\beta = 2 - 2\alpha$, $\psi_{2\alpha - 2}(\zeta)$ has finite variance iff $\beta < \f 1 2 \iff \alpha > \f 3 4$.

If $\alpha = 1$, then $${\underline}{{\boldsymbol{\oldchi}}}= {{\boldsymbol{\oldchi}}}(1 + \f 1 2\sigma_v^2 \sigma_w^2).$$ So $\p {{\boldsymbol{\oldchi}}}{l-m} / \p {{\boldsymbol{\oldchi}}}l = \Theta(1) B^m$ for $B = 1 + \f 1 2\sigma_v^2 \sigma_w^2$.

If $\f 1 2 < \alpha < 1$, then ${\underline}{{\boldsymbol{\oldchi}}}/{{\boldsymbol{\oldchi}}}- 1$ is $$\begin{aligned}
 &\phantom{=} \sigma_v^2 \sigma_w^2 {\mathrm{V}}\dot \phi( {\mathbbm{q}})\\
    &= \sigma_v^2 \sigma_w^2 \alpha^2 {\mathsf{c}}_{\alpha - 1} {\mathbbm{q}}^{\alpha - 1}\\
    &= \sigma_v^2 \sigma_w^2 \alpha^2 {\mathsf{c}}_{\alpha - 1} (\sigma_w^2 {\mathbbm{p}})^{\alpha - 1} + \Theta({\mathbbm{p}}^{\alpha - 2})\\
    &= \sigma_v^2 \sigma_w^{2\alpha} \alpha^2 {\mathsf{c}}_{\alpha - 1}(K_1 l^{{\f 1 {1-\alpha}}}- K_2 l^{{\f \alpha {1-\alpha}}}\log l + o(l^{{\f \alpha {1-\alpha}}}\log l))^{\alpha - 1} + \Theta(l^{\f{\alpha - 2}{1 - \alpha}}) & \text{by \cref{thm:pDynamicLT1ReLU}}\\
    &= \sigma_v^2 \sigma_w^{2\alpha} \alpha^2 {\mathsf{c}}_{\alpha - 1} [K_1^{\alpha - 1} \inv l + \Theta(l^{-2} \log l)] + O(l^{-3})\\
    &= \sigma_v^2 \sigma_w^{2\alpha} \alpha^2 {\mathsf{c}}_{\alpha - 1} K_1^{\alpha - 1} \inv l + \Theta(l^{-2} \log l)\\
    &= R \inv l + \Theta(l^{-2} \log l)\end{aligned}$$ where $R = \sigma_v^2 \sigma_w^{2\alpha} \alpha^2 {\mathsf{c}}_{\alpha - 1} K_1^{\alpha - 1} = \f{\alpha^2}{(1-\alpha)(2 \alpha - 1)} $ and $K_1 = [\sigma_v^2 \sigma_w^{2\alpha} {\mathsf{c}}_\alpha (1 - \alpha)]^{\f 1 {1-\alpha}}$. So $$\begin{aligned}
{\underline}{{\boldsymbol{\oldchi}}}&= {{\boldsymbol{\oldchi}}}\exp(R \inv l + \Theta(l^{-2} \log l))\\
\p {{\boldsymbol{\oldchi}}}{l -m} &= \Theta(1)\p {{\boldsymbol{\oldchi}}}{l} \lp\f l {l-m} \rp^R\end{aligned}$$ as desired.

The proof is similar to that of \[thm:dalethExpSqrtTanhAllGrad\].

[^1]: Work done while at Harvard University

[^2]: Under simplified conditions, @daniely_toward_2016 showed that there exists a fixed point for any “well-behaved” activation function in a feedforward net. However, this result does not apply to architectures with residual connections.

[^3]: Note that in practice, to avoid the diverging gradient $\dot \psi_\alpha(x) \to \infty$ as $x \to 0$, we can use a tempered version $\Psi_\alpha(x)$ of $\alpha$-ReLU, defined by $\Psi_\alpha(x) = (x + {\epsilon})^\alpha - {\epsilon}^\alpha$ on $x > 0$ and 0 otherwise, for some small ${\epsilon}> 0$. The conclusions of this paper on $\psi_\alpha$ should hold similarly for $\Psi_\alpha$ as well.

[^4]: @daniely_toward_2016 called the version of ${\mathrm{W}}\phi$ with fixed $\rho = 1$ the “dual function” of $\phi$.

[^5]: [A more natural visualization is to graph $\p {\mathbf{e}}l - {\mathbf{e}}^*$ versus $l^{-\delta^*}$, but because of floating point precision, $\p {\mathbf{e}}l - {\mathbf{e}}^*$ doesn’t converge to 0, but a small number close to 0, so that the log-log plot wouldn’t look like what is expected.\[footnote:plotDelta\]]{}

[^6]: [A more natural visualization is to graph $\p {\mathbf{e}}l - {\mathbf{e}}^*$ versus $l^{-\delta^*}$, but because of floating point precision, $\p {\mathbf{e}}l - {\mathbf{e}}^*$ doesn’t converge to 0, but a small number close to 0, so that the log-log plot wouldn’t look like what is expected.\[footnote:plotDelta\]]{}

[^7]: Our derivations actually apply to all $\alpha \in (\f 1 2, 1]$, where at $\alpha = \f 1 2$, the expected norm of the gradient diverges within our mean field formalism. However, at $\alpha \le \f 3 4$, the variance of the gradient already diverges (\[thm:dalethInfVarAlphaReLU\]), so we cannot expect the empirical values to agree with our theoretical predictions. But in fact, empirically our theoretical predictions seem to form an upper bound on the gradient norms (see \[fig:alphaReLUTheoryVsEmpirics\]).

[^8]: the contour for $\p {\mathbf{p}}l$ is similar, but its slopes are slightly off from the heatmap contours.
