Mean Field Residual Networks: On the Edge of Chaos
Greg Yang Samuel S. Schoenholz
cs.NE cond-mat.dis-nn cs.LG math.DS nlin.CD


Impact of energy dissipation on interface shapes and on rates for dewetting from liquid substrates

    
Dirk PeschkaE-mail: dirk.peschka@wias-berlin.de ^1,

Stefan Bommer^2, 

Sebastian Jachalski^1, 


Ralf Seemann^2, Barbara Wagner^1


^1 Weierstraß Institute, Mohrenstrasse 39, 10117 Berlin, Germany 

^2 Experimental Physics, Saarland University, 66123 Saarbrücken, Germany


    December 30, 2023
=======================================================================================================================================================================================================================================================================================


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.




§ INTRODUCTION


Previous works <cit.> 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 <cit.> 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 <cit.> (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 α-ReLU, the αth-power of ReLU, for α < 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 (α-)ReLU resnets is determined by expressivity (how far can two different input vectors be pulled apart) (see <ref>).
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 <cit.>), so that common initialization schemes like Xavier <cit.> or He <cit.> 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 <cit.>.

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



Consider a vanilla feedforward neural network of L layers, with each layer l having N l neurons; here layer 0 is the input layer.
For the ease of presentation we assume all hidden layer widths are the same N l = N for all l > 0.
Let x 0 = ( x 0_1, …,  x 0_ N 0) be the input vector to the network, and let x l for l > 0 be the activation of layer l.
Then a neural network is given by the equations
	
    x l _i    = ϕ( h l _i),     h l _i    = ∑_j=1^N  w l _ij x l-1_j +  b l_i

	where
	[label=(*)]
	    
  * h l is the pre-activation at layer l,
	    
  * w l is the weight matrix,
	    
  * b l is the bias vector, and
	    
  * ϕ 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
	
    x_i    = ϕ(h_i),     h_i    = ∑_j w_ij x_j + b_i

where ∙ implicitly denotes ∙ l, and ∙ denotes ∙l-1 (and analogously, ∙ denotes ∙l+1).

A series of papers <cit.> investigated the “average behavior” of random neural networks sampled via w l_ij∼(0, σ_w^2/N),  b l_i ∼(0, σ_b^2), for fixed parameters σ_w and σ_b, independent of l.
Consider the expectation of 1 N ∑_i=1^N x^2_i, the normalized squared length of x, over the sampling of w and b.
<cit.> showed that this quantity converges to a fixed point exponentially fast for sigmoid nonlinearities.
Now suppose we propagate two different vectors x 0 and ( x 0)' through the network.
<cit.> also showed that the expectation of the normalized dot product 1 N ∑_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, <cit.>, 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 σ_w and σ_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, <cit.> 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 σ_∙, 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 <cit.> and later refined by <cit.>, 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
	
    x_i    = ϕ(h_i) +  x,     h_i    = ∑_j w_ij x_j + b_i.

A (full) residual network (FRN) in addition has an affine connection given by weights v and biases a from the nonlinearity ϕ(h) to the next layer:
	
    x_i    = ∑_j v_ijϕ( h_j) +  x_i + a_i,      h_i    = ∑_j w_ij x_j + b_i


We are interested in the “average behavior” of these network when the weights and biases, w l_ij,  b l_i,  v l_ij, and a l_i are sampled i.i.d. from Gaussian distributions resp. with standard deviations σ_w, σ_b, σ_v, and σ_a, independent from l.
Here we take the variance of w l_ij to be σ_w^2/N so that the variance of each h_i is  σ_w^2, assuming each x_j is fixed (similarity for v l_ij).
Such an initialization scheme is standard in practice.

We make several key “physical assumptions” to make theoretical computations tractable:
[Symmetry of activations and gradients]
	(a) We assume ( h l _i)^2  =  ( h l _j)^2 and ( x 0 _i)^2  =  ( x 0 _j)^2 for any i, j, l.
	(b) We also assume that the gradient E/ x l _i with respect to the loss function E satisfies ( E/ x l _i)^2  =  ( E/  x l _j)^2 for any i, j, l.

One can see that <ref>(a) is satisfied if the input x 0 ∈{± 1}^N and <ref>(b) is satisfied if <ref> below is true and the gradient at the last layer E/ x L ∈{± 1 }^N.
But in general it is justified both empirically and theoretically as an approximation, because ( h l _i)^2 - ( h l _j)^2 stays about constant with l, but ( h l _i)^2 and ( 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 ( x l _i)^2 and ( E/ h l _i)^2.
[Gradient independence]
	(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, E/ x l _i, is independent from all activations h l _j and x l-1 _j from the previous layer.

<ref>(a) was first made in <cit.> for computing the mean field theory of gradients for feedforward tanh networks.
This is similar to the practice of feedback alignment <cit.>.
Even though we are the first to explicitly formulate <ref>(b), in fact it was already applied implicitly in the gradient calculations of <cit.>.
Note that a priori <ref>(b) is not true, as E/ x l _i depends on ϕ̇( h l+1 _k) for every k, which depend on h l _j for each j, and which depends on 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 (<ref>) to ameliorate notation confusion.

Fix an input x 0.
Define the length quantities l :=  ( h l_1)^2 and l :=  ( x l_1)^2 for l > 0 and 0 =  x 0^2/N.
Here the expectations ∙ are taken over all random initialization of weights and biases for all layers l, as N →∞ (large width limit).

Note that in our definition, the index 1 does not matter by <ref>.

Fix two inputs x 0 and x 0'.
We write ∙' to denote a quantity ∙ with respect to the input x 0'.
Then define the correlation quantities γ l:=  h l_1  h l_1' and λ l:=  x l_1  x l_1' for l > 0 and γ 0 =  x 0 · x 0 ' / N, where the expectations ∙ are taken over all random initialization of weights and biases for all layers l, as N →∞ (large width limit).
Again, here the index 1 does not matter by <ref>.
By metric expressivity, we mean l :=  1 2N x l -  x l '^2 =  1 2N ( x l^2 +  x l '^2  - 2  x l · x l ') =  1 2 ( l +  l ') - γ l.
Additionally, define the cosine distance quantities l := γ l / √( l  l ') and l := λ l / √( l  l'), and we will also call l angular expressivity.

In this paper, for the ease of presentation, we assume 0 =  0 '.
Then, as we will see, l =  l',  l =  l' for all l, and as a result, l = γ l /  l and l =  l - γ l = (1 -  l)  l.


Fix an input x 0 and a gradient vector ( E/  x L_i)_i of some loss function E with respect to the last layer x L.

Then define the gradient quantities ℸ l:=  ( E/ x l _1)^2 , χ l _∙ :=  ( E/∙ l _1)^2 for ∙ = a, b, and χ l _∙ :=  ( E/∙ l _11)^2 for ∙ = w, v.
Here the expectations are taken with <ref> in mind, over both random initialization of forward and backward weights and biases, as N →∞ (large width limit).
Again, the index 1 or 11 does not matter by <ref>.





  
Asymptotic notations.
The expressions f = O(g)  g = Ω(f) have their typical meanings, and f = Θ(g) iff f = O(g), g = O(f).
We take f(x) = Õ(g(x))  g(x) = Ω̃(f(x)) to mean f(x) = O(glog^k x) for some k ∈ (this is slightly different from the standard usage of Õ), and f = Θ̃(g)  f = Õ(g)  g = Õ(f).
We introduce a new notation: f = (g) if f(x) = O(g(x) · x^) and f(x) = Ω(g(x) · x^-), as x →∞, for any > 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 , , γ, λ, ℸ for different architectures and different activation functions.
We use these equations to investigate the dynamics of  and , the key quantities in the forward pass, and the dynamics of ℸ, the key quantity in the backward pass.

The cosine distance  in some sense measures the angular geometry of two vectors.
If = 1, then the vectors are parallel; if = 0, then they are orthogonal.
Just as in <cit.> and <cit.>, we will show that in all of the architectures and activations we consider in this paper, l converges to a fixed point ^* as l →∞ [Under simplified conditions, <cit.> 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.].
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 <cit.> and <cit.> showed that the convergence rate to ^* is exponential in a vanilla network, the convergence rate is rather only polynomial in residual networks, for tanh and α-ReLU (<ref>) nonlinearities; see <ref>, <ref>, <ref>, and <ref>.

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 l converges to 0, corresponding to “chaos”, (resp. 1, corresponding to “stability”), for the number of layers usually seen in practice, l will reside well away from 0 (resp. 1).

Similarly, the quantity  measures the metric geometry of two vectors.
The evolution of 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 α-ReLU (α < 1) nonlinearities,  varies only polynomially with l.

On the other hand, ℸ 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 ℸ l as l → 1; see <ref> and <ref>.
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 ℸ 0 ≈exp(Θ(√(l))) ℸ l (<ref>), α-ReLU residual networks for α < 1 have ℸ 0 ≈(l) ℸ l (<ref>).
Instead of E/ x_i, we may also consider the size of gradients of actual trainable parameters.
For tanh and α-ReLU with α < 1, they are still subexponential and polynomial (<ref>).
On the other hand, while ℸ 0 = exp(Θ(l))ℸ l for a ReLU resnet, its weight gradients have size independent of layer, within O(1) (<ref>)!
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 ℸ 0 / ℸ L; whereas for ReLU and α-ReLU, the optimal initialization is obtained by maximizing  without running into numerical issues (with floating point computation).

See <ref> for details.

Over the course of our investigation of α-ReLU, we derived several new identities involving the associated kernel functions, first defined in <cit.>, which relate them to the zeroth Bessel functions (<ref>).



§ THEORETICAL RESULTS


In what follows in the main text, we assume σ_∙ > 0 for all ∙ = 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 ϕ is tanh-like if ϕ is antisymmetric (ϕ(-x) = -ϕ(x)), |ϕ(x)| ≤ 1 for all x, ϕ(x) ≥ 0, ∀ x ≥ 0, and ϕ(x) monotonically increases to 1 as x →∞.


Define the α-ReLU ψ_α(x) = x^α if x > 0 and 0 otherwise. [
	Note that in practice, to avoid the diverging gradient ψ̇_α(x) →∞ as x → 0, we can use a tempered version Ψ_α(x) of α-ReLU, defined by Ψ_α(x) = (x + )^α - ^α on x > 0 and 0 otherwise, for some small > 0.
	The conclusions of this paper on ψ_α should hold similarly for Ψ_α as well.]




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 <ref> (see appendix for proofs).
They leverage certain integral transforms [<cit.> called the version of ϕ with fixed ρ = 1 the “dual function” of ϕ.] as in the following

Define the transforms  and  by ϕ(q) := [ϕ(z)^2: z ∼(0, q)] and ϕ(ρ, ν) := [ϕ(z)ϕ(z'): (z, z') ∼(0, [ ρ ν; ν ρ ])].


These recurrences are able to track the corresponding quantities in practice very well.
For example, <ref> 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 σ_v = 1 and σ_a = 0) and α-ReLU/FRN as well (see <ref>).

As mentioned in previous sections, we seek to characterize the long term/high depth behavior of all of the quantities defined in <ref>.
To do so, we solve for the asymptotics of the recurrences in <ref>, where ϕ is instantiated with tanh or α-ReLU.
Our main dynamics results are summarized in <ref>.






 §.§ Tanh









  
Forward dynamics.
When ϕ = tanh, l and l increase as Θ(l) in either RRN or FRN (<ref>), as one might expect by observing that tanh( ) → 1 as →∞ so that, for example in the RRN case, the recurrence = tanh( ) + becomes = 1 +.
This is confirmed graphically by the black lines of the leftmost chart of <ref>.
We carefully verify that this intuition is correct in its proof in the appendix, and find that in fact l ∼ l in the RRN case and l ∼ (σ_v^2 + σ_a^2)l in the FRN case.

What about γ l?
The middle chart of <ref> shows that over time, l = γ l /  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 ^* of  or not.  
We prove that, in fact, all trajectories of  not starting at 1 do converge to a single fixed point, but only at a polynomial rate, in both the RRN and FRN cases (<ref> and <ref>); we can even explicitly compute the fixed point and the rate of convergence:
For FRN, there is a unique stable fixed point ^* < 1 determined by the equation

    ^* =  1 σ_v^2 + σ_a^2[σ_v^2  2 πarcsin^*  + σ_a^2],

and |^* -  l| decreases like l^-δ^*, where

    δ^* := 1 -  2 π 1 √(1 - (^*)^2)σ_v^2 σ_v^2 + σ_a^2.

Since ^* < 1, = (1 - )  = Θ() = Θ(l).
The case of RRN can be viewed as a special case of the above, setting σ_v^2 = 1 and σ_a^2 = 0, which yields ^* = 0 and δ^* = 1 -  2 π.
We observe that both ^* and δ^* only depend on the ratio ρ := σ_a/σ_v, so in <ref> we graph these two quantities as a function of ρ.
^* and δ^* both increase with ρ and asymptotically approach 1 and 1/2 respectively from below.
When ρ = σ_a = 0, ^* = 0 and δ^* = 1 -  2 π.
Thus the rate of convergence at its slowest for tanh/FRN is δ^* = 1 -  2 π≈ 0.36338, where asymptotically the network tends toward a chaotic regime ^* = 0, corresponding to a large weight variance and a small bias variance; it at its fastest is δ^* = 1/2, where asymptotically the network tends toward a stable regime ^* = 1, corresponding to a large bias variance and small weight variance.
We verify δ^* by comparing l - l-1 to l^-δ^* - 1 in log-log scale.
If l = Θ(l^-δ^*), then l - l-1 = Θ(l^-δ^* - 1) and should obtain the same slope as l^-δ^* - 1 as l →∞.
The middle figure of <ref> ascertains that this is indeed the case, starting around layer number 400. 





  
Backward dynamics. Finally, we show that the gradient is approximated by 

    ℸm   = exp(𝒜(√(l) - √(m)) + O(log l - log m))ℸl⋆

where 𝒜 =  4 3 √( 2 π)σ_w in the RRN case and 𝒜 =  4 3 √( 2 π)σ_v^2 σ_w√(σ_v^2 + σ_a^2) in the FRN case (<ref> and <ref>).
The rightmost plot of <ref> verifies that indeed, for large l ≥ 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 𝒜 can be decomposed into 𝒜 =  4 3 √( 2 π)·σ_v ·σ_w · (1 + σ_a^2/σ_v^2)^-1/2.
Consider the ratio ρ := σ_a/σ_v.
If ρ≫ 1, then ^* ≈ 1 (<ref>), meaning that the typical network essentially computes a constant function, and thus unexpressive; at the same time, large ρ makes 𝒜 small, and thus ameliorating the gradient explosion problem, making the network more trainable.
On the other hand, if ρ≪ 1, then ^* ≈ 0 (<ref>), 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 <cit.>; at the same time, small ρ increases 𝒜, worsening the gradient explosion problem, making the network less trainable.
This is the same expressivity-trainability tradeoff discussed in <cit.>.



 §.§ α-ReLU




  
Forward dynamics. As with the tanh case, to deduce the asymptotic behavior of random α-ReLU resnets, we need to understand the transforms ψ_α and ψ_α.
Fortunately, ψ_α has a closed form, and ψ_α has been studied before <cit.>.
In particular, if α > - 1 2, then
ψ_α( ) = _α^α, where _α is a constant with a closed form given by <ref>.
In addition, by <cit.>, we know that ψ_α( , ) = ψ_α( ) _α() for _α given in <ref>.
<ref> shows a comparison of _α for different αs along with the identity function.

Substituting in _α^α for ψ_α, we get a difference equation -  = σ_v^2 _α (σ_w^2  + σ_b^2)^α + σ_a^2 governing the evolution of .
This should be reminiscent of the differential equation Ṗ(l) = C P(l)^α, which has solution ∝ l^1/(1-α) for α < 1, and ∝exp(Cl) when α = 1.
And indeed, the solutions l to these difference equations behave asymptotically exactly like so (<ref>).
Thus ReLU behaves very explosively compared to α-ReLU with α<1.
In fact, in simulations, for σ_w^2 = 1.69 and σ_v^2 = 1.5, the ReLU resnets overflows into s after around 100 layers, while there's no problem from any other kind of networks we consider.

Regardless, α-ReLU for all α massages l toward a fixed point ^* that depends on α.
 When ϕ = ψ_1, the standard ReLU, l converges to 1 asymptotically as C l^-2 for an explicit constant C depending on σ_v and σ_w only (<ref>), so that = (1 - ) = Θ(l^-2exp(Θ(l))) = exp(Θ(l)).
When ϕ = ψ_α for α < 1, then l converges to the nonunit fixed point ^* of _α at a rate of (l^-μ), where μ = (1-_α(^*))/(1-α) is independent of the variances (<ref>), so that = Θ().
These rates are verified in <ref>.



  
Backward dynamics. Finally, we have also characterized the rate of gradient growth for any α∈ ( 3 4, 1].
[Our derivations actually apply to all α∈ ( 1 2, 1], where at α =  1 2, the expected norm of the gradient diverges within our mean field formalism.
However, at α≤ 3 4, the variance of the gradient already diverges (<ref>), 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 <ref>).
]
In the case of α = 1, the dynamics of ℸ is exponential, the same as that of , ℸl-m = ℸl B^m where B = 1 2 σ_v^2 σ_w^2 + 1.
For α∈ ( 3 4, 1), the dynamics is polynomial, but with different exponent in general from that of the forward pass: ℸl-m = Θ(1) ℸl (l/(l-m))^R for R = α^2(1-α)(2 α - 1), where the constants in Θ(1) do not depend on l or m.
This exponent R is minimized on α∈ [ 3 4, 1) at α = 3/4, where R = 9/2 (but on α∈ ( 1 2, 1) it is minimized at α = 2/3, where R = 4); see <ref>.
These exponents are verified empirically in <ref>.

Looking only at ℸ and the gradients against the biases, it seems that ReLU suffers from a dramatic case of exploding gradients.
But in fact, because ℸ gains a factor of B moving backwards while  loses a factor of B, the gradient norm χ_w^(l-m) (and similarly for χ_v^(l-m)) is independent of how far, m, the gradient has been propagated (<ref>) — 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 α-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 α-ReLU is useful to the extent that smaller α avoids numerical issues with exponentiating forward and backward dynamics, but the best performance is given by the largest α that avoids them (<ref>(c, d)); in fact, the metric expressivity , determines performance, not gradient explosion (see α-ReLU experiments).



§ EXPERIMENTAL RESULTS






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 ℸ 0 / ℸ L); for (α-)ReLU resnets, it is determined by how expressive the random network is (measured by the metric expressivity L).
We hypothesize this is because in tanh resnets, the gradient dynamics is much more explosive than the expressivity dynamics (exp(Θ(√(l))) vs Θ(l)), whereas for ReLU it's somewhat the opposite (χ_w, χ_v = Θ(1) vs = exp(Θ(l))).



  
Tanh, vary σ_w.
We train a grid of reduced and full tanh resnets on MNIST, varying the variance σ_w^2 and the number of layers (for FRN we fix σ_v = 1).
The results are indicated in <ref>(a, b).
We see that in either model, deeper resnets favor much smaller σ_w than shallower ones.

The white dotted lines in <ref>(a, b) confirm our theory: according to <ref>, for the same gradient ratio R = ℸ 0 / ℸ L, we want log R ≈σ_w √(L).


Indeed, the white dotted lines in <ref>(a, b) trace out such a level curve and it remarkably pinpoints the largest σ_w that gives the optimal test set accuracy for each depth L.
Why isn't the best initialization given by R = 1 σ_w = 0?
We believe that when L and/or σ_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 σ_a^2/σ_v^2.
As suggested in the analysis of <ref>, the ratio ρ^2 = σ_a^2/σ_v^2 determines the fixed point ^* 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 σ_v with ρ such that σ_v/√(1 + ρ^2) is kept constant, so that the leading term of the log gradient ratio is kept approximately equal for each L and ρ.
<ref>(c) shows the test accuracies of a grid of tanh FRN initialized with such an ensemble of σ_∙s.
What stands out the most is that performance is maximized essentially around a fixed value of L regardless of ρ, which shows that indeed gradient dynamics determines the initialization quality in tanh resnets.
There is also a minor increase in performance with increasing ρ regardless of L; this is counterintuitive as increasing ρ means “decreasing expressivity.”
It is currently not clear what accounts for this effect.




  
ReLU, vary σ_w
We train a grid of ReLU FRN on MNIST, varying σ_w^2 ∈ [0, 1.5] while fixing σ_v^2 = 1, σ_a^2 = σ_b^2 =  1 2.
The resulting test set accuracies are shown in <ref>(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 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 χ 0 / χ L ≥ 1, then we would be wrong — in fact it is exactly the opposite.
In this case, the dynamics of l,  l, and χ 0 / χ l are approximately the same (all exp(Θ(l)) with the same hidden constants), and optimal performance corresponds to the highest L, L, and χ 0 / χ L without running into s.






  
α-ReLU, vary α.
We similarly trained a grid of α-ReLU FRN on MNIST, varying only α and the depth, fixing all σ_∙.
<ref>(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 χ 0 / χ L (white dashed line) now gives no information at all on the test set accuracy.
In contrast, the contours for l succeeds remarkably well at this prediction (yellow/green lines).
[the contour for l is similar, but its slopes are slightly off from the heatmap contours.]
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  dynamics to be predictive of neural network performance.




§ CONCLUSION

In this paper, we have extended the mean field formalism developed by <cit.> 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 <cit.> or He <cit.> 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 <ref> 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


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.
plainnat











§ ADDITIONAL FIGURES

In figures appearing in the appendix,  means χ (due to legacy reasons).












§ A LISTING OF MAIN THEOREMS



 §.§ Tanh




  §.§.§ Reduced Residual Network

lemmapqrecurrence
Suppose ϕ is antisymmetric.
Then in an RRN,  and  satisfy the recurrence

    = σ_w^2  + σ_b^2
       = ϕ( ) + .




thmpqlinear
Suppose ϕ is tanh-like.
Assume RRN architecture.

	
  * If σ_w = 0, then l = l ϕ( σ_b^2) +  0 and l = σ_b^2.
	
  * If σ_w > 0, lim_l →∞ l/ l = 1 and lim_l →∞ l /(σ_w^2 l) = 1.
	If ϕ = tanh, then we can obtain more terms of the asymptotic expansions:
	
    l    = l - 2 C σ_w^-1 l^1/2 - C^2 σ_w^-2log l + O(1)
     l    = σ_w^2 l - 2 C σ_w l^1/2 - C^2 log l + O(1)

	as l →∞, where C = √(2 / π).





thmlambdagammarecurrence
Suppose ϕ is antisymmetric.
Then in an RRN, λ and γ satisfy the recurrence

    λ   = σ_w^2 γ + σ_b^2
    γ   = ϕ(, λ) + γ.



thmedynamics
Suppose ϕ is a tanh-like nonlinearity in an RRN.
Assume 0 < 1.

	
  * If σ_w = 0, then γ l = l ϕ( σ_b^2, σ_b^2) + γ 0 = l ϕ( σ_b^2) + γ 0 and λ l = σ_b^2, so that l → 1 and 1 -  l = Θ(l^-1).
	As a result, l =  l (1 -  l) = Θ(1).
	
  * If σ_w > 0, then γ l = (l^ 2π), and l → 0 like (l^ 2 π - 1).
	Thus l = Θ( l) = Θ(l).




thmdalethRecReduced
For any nonlinearity ϕ in an RRN, under assumptions <ref> and <ref>, whenever ϕ̇^2(ζ) has finite variance for Gaussian variable ζ,

    ℸ   = (σ_w^2  ϕ̇( ) + 1)ℸ,   χ_b    = ℸϕ̇( ),   χ_w    = ℸϕ̇( ) .




thmdalethExpSqrtTanh
For ϕ = tanh in an RRN,

	
  * If σ_w = 0, ℸ m = ℸ l for all l, m.
	
  * If σ_w > 0, 
    log(ℸm/ℸ l)    = 𝒜 (√(l) - √(m)) + ℬ (log l - log m) + O(1)

	where 𝒜 =  4 3 √( 2 π)σ_w and ℬ =  4 3π - σ_w^2  4 9π.





thmdalethExpSqrtTanhAllGrad
Suppose ϕ = tanh.
Then in an RRN

	
  * If σ_w = 0, χ l _b = ℸ L ϕ̇( σ_b^2) and χ l _w = ℸ L ϕ̇( σ_b^2) ((l-1) ϕ( σ_b^2) +  0), where L is the last layer.
	
  * If σ_w > 0, 
    log(χ m _b / χ l _b)    = 𝒜 (√(l) - √(m)) + ℬ_b (log l - log m) + O(1)
    log(χ m _w / χ l _w)    = 𝒜 (√(l) - √(m)) + ℬ_w (log l - log m) + O(1)

	where 𝒜 =  4 3 √( 2 π)σ_w (same as 𝒜 in <ref>) and ℬ_b = ℬ +  1 2, ℬ_w = ℬ -  1 2, with ℬ =  4 3π - σ_w^2  4 9π (same as ℬ in <ref>).




  §.§.§ Full Residual Network


thmfullResPQRec 
For any nonlinearity ϕ in an FRN,

    = σ_w^2  + σ_b^2
       = σ_v^2 ϕ( ) + σ_a^2 +



thmpIsLinearTanh
Suppose ϕ is tanh-like.
Assume the FRN architecture.

	
  * If σ_w = 0, then l = (σ_v^2 ϕ( σ_b^2) + σ_a^2)l + 0, and l = σ_b^2.
	
  * If σ_w > 0, then l = b_0 l + b_1 l^1/2 + b_2 log l + O(1), where
	
    b_0    = σ_v^2 + σ_a^2
    
    	b_1    = -2C σ_v^2 σ_w^-1√(σ_v^2 + σ_a^2)
    
    	b_2    = -C^2 σ_v^4 σ_w^-2(σ_v^2 + σ_a^2)^2

	and C = √( 2 π).
	Additionally, l = σ_w^2 b_0 l + σ_w^2 b_1 l^1/2 + σ_w^2 b_2 log l + O(1).






thmLGRecFullRes
For any nonlinearity ϕ, in an FRN

    λ   = σ_w^2 γ + σ_b^2
    γ   = σ_v^2 ϕ(, λ) + σ_a^2 + γ



thmeDynamicsFullResTanh 
	Assume ϕ = tanh in an FRN.
	Suppose 0 < 1.
	
		
  * If σ_w = 0, then λ l = σ_b^2 and γ l = l (σ_v^2 ϕ( σ_b^2, σ_b^2) + σ_a^2) + γ 0 = l (σ_v^2 ϕ( σ_b^2) + σ_a^2) + γ 0.
		Thus l → 1 and 1 -  l = Θ(l^-1).
		As a result, l =  l (1 -  l) = Θ(1).
		
  * If σ_w > 0, then l converges to the unique fixed point ^* ≠ 1 determined by the equation
		
    ^* =  1 σ_v^2 + σ_a^2[σ_v^2  2 πarcsin^*  + σ_a^2].

		Furthermore, l converges to ^* polynomially:
		| l - ^*| is (l^-δ^*), where
		
    δ^* := 1 -  2 π 1 √(1 - (^*)^2)σ_v^2 σ_v^2 + σ_a^2∈ [ 2 π - 1,  1 2)

		Since ^* < 1, l = Θ( l) = Θ(l).
	


thmdalethRecFull 
For any nonlinearity ϕ in an FRN, under assumptions <ref> and <ref>, whenever ϕ̇(ζ)^2 has finite variance for Gaussian variable ζ,

    ℸ   = (σ_v^2σ_w^2  ϕ̇( ) + 1)ℸ,   χ_b    = σ_v^2ℸϕ̇( ),
    χ_w    = σ_v^2ℸϕ̇( ) ,   χ_v    = ℸϕ( ),   χ_a    = ℸ



thmdalethExpSqrtTanhFullRes
Assume ϕ = tanh in an FRN.

	
  * If σ_w = 0, ℸ m = ℸ l for all l, m.
	
  * If σ_w > 0, then for l ≥ m ≥ 0,
	
    log(ℸm / ℸ l) = 𝒜 (√(l) - √(m)) + ℬ (log l - log m) + O(1)

	where
	
    𝒜   =  4 3 √( 2 π)σ_v^2 σ_w√(σ_v^2 + σ_a^2)
    ℬ   =  4 9πσ_v^4 σ_v^2 + σ_a^2 3 σ_v^2 + σ_a^2 - σ_w^2




<ref> shows empirical verification of the asymptotic expansion of ℸ for various values of σ_∙s.



thmdalethExpSqrtTanhFullResAllGrad
Suppose ϕ = tanh in an FRN.

	
  * If σ_w = 0, then
	
    χ l _b    = σ_v^2 ℸ L ϕ̇( σ_b^2)
    χ l _w    = σ_v^2 ℸ L ϕ̇( σ_b^2) ( (σ_v^2 ϕ( σ_b^2) + σ_a^2)(l-1) + 0)
    χ l _v    = ℸ L ϕ( σ_b^2)
    χ l _a    = ℸ L.

	
  * If σ_w > 0, then for l ≥ m ≥ 0,
	
	
    log(χ m _b / χ l _b)    = 𝒜 (√(l) - √(m)) + ℬ_b (log l - log m) + O(1)
    log(χ m _w / χ l _w)    = 𝒜 (√(l) - √(m)) + ℬ_w (log l - log m) + O(1)
    log(χ m _a / χ l _a)    = 𝒜 (√(l) - √(m)) + ℬ (log l - log m) + O(1)
    log(χ m _v / χ l _v)    = 𝒜 (√(l) - √(m)) + ℬ (log l - log m) + O(1)

	where 𝒜 =  4 3 √( 2 π)σ_v^2 σ_w√(σ_v^2 + σ_a^2) and ℬ =  4 9πσ_v^4 σ_v^2 + σ_a^2 3 σ_v^2 + σ_a^2 - σ_w^2 are as in <ref> and ℬ_b = ℬ +  1 2 and ℬ_w = ℬ -  1 2.
	






 §.§ α-ReLU

lemmaVtPsiAlpha
If α > - 1 2, then
ψ_α( q) = _α q^α, where _α =  1 √(π) 2^α - 1Γ(α+  1 2).


Note that if α≤ -  1 2, then ψ_α( q) is not defined (its defining integral does not converge).



  §.§.§ Full Residual Network

By <ref> and <ref>, we have the length recurrences

    = σ_w^2  + σ_b^2
       = σ_v^2 _α^α + σ_a^2 +


thmpDynamicAlphaReLU
Suppose we have the nonlinearity ϕ = ψ_α.
The in an FRN:
If α = 1, then l = Θ((1 + σ_v^2 σ_w^2/2)^l), with the hidden constant depending on the initial condition.
If 0 < α < 1, then l = Θ(l^ 1 1- α).
More precisely,
lim_l →∞/ l^ 1 1-α = [σ_v^2 σ_w^2α_α (1 - α)]^ 1 1-α.

<ref> empirically verifies the asymptotics for α=1 for various σ_v and σ_w.


Similarly, by <ref>, if = ', then

    λ   = σ_w^2 γ + σ_b^2
    γ   = σ_v^2 ^αψ_α( 1, ) + σ_a^2 + γ



thmReLUSquaredConvergence
Suppose ϕ = ψ_1.
Then in an FRN, l → 1 and 1 -  l ∼ [ 1 4 σ_v^2 σ_w^2  B U l]^-2 for B = 1 + σ_v^2 σ_w^2/2 and U = 2√(2)3π.
As a result, l = (1 - ł)  l = Θ(l^-2exp(Θ(l))) = exp(Θ(l)).


thmalphaReLUeConvergence 
Suppose ϕ = ψ_α for 0 < α < 1 in an FRN.

Then  converges to the unique nonunit fixed point ^* of _α, and |^* -  l| is (l^-μ), where μ = (1-_α(^*))/(1-α).
Additionally, l = Θ( l) = Θ(l^1/(1-α)).


<ref> verifies empirically that ^* is indeed the fixed point of l.
<ref> verifies empirically the convergence rate l^-μ.
<ref> plots _α(^*) and μ versus α.
It certainly looks like μ =  1 2 (1 - α), but we have no proof for it.
Based on this conjecture, we see there is a “discontinuity” of μ at α = 1: μ→ 0 as α→ 1, but for α = 1, the actual convergence dynamics has exponent -2 by <ref>.





Because of the following theorem, we cannot expect the equations of <ref> to hold for α≤ 3 4.
thmdalethInfVarAlphaReLU 
	Suppose we have the nonlinearity ψ_α in an FRN.
	(ψ̇_α(ζ)^2) diverges for any Gaussian variable ζ with mean 0 if α≤ 3 4 but is finite if α >  3 4.


thmdalethDynamicsAlphaReLU 
Suppose we have the nonlinearity ψ_α in an FRN.
If α = 1, then ℸl-m = ℸl( 1 2 σ_v^2 σ_w^2 + 1)^m.
If α∈ ( 3 4, 1), then ℸl-m = Θ(1) ℸl (l/(l-m))^R for R = α^2(1-α)(2 α - 1), where the constants in Θ(1) do not depend on l or m.

This exponent α^2(1 - α)(2α - 1) is minimized at α =  3 4 on α∈ (3/4, 1), where the value is 9 2 (and at α =  2 3 on α∈ (1/2, 1), where the value achieved is 4) (<ref>(a)).

As a corollary,
thmalphaReLUAllGradients
If ϕ = ψ_1 in an FRN, then for l ≥ m ≥ 0,

    χl-m _b    = Θ(1) ℸ l B^m,   χl-m _w    = Θ(1) ℸ l B^l,
    χl-m _v   = Θ(1) ℸ l B^l,   χl-m _a    = Θ(1) ℸl B^m.

where B = 1 + σ_v^2σ_w^2/2.

If ϕ = ψ_α in an FRN, for α < 1, then for l ≥ m ≥ 0,

    χl-m _b    = Θ(1) ℸ l l^R (l-m)^-R-1,   χl-m _w    = Θ(1) ℸ l l^R (l-m)^α1-α - R,
    χl-m _v    = Θ(1) ℸ l l^R (l-m)^α1-α - R,   χl-m _a    = Θ(1) ℸl (l/(l-m))^R.



<ref> verifies the backward asymptotic dynamics empirically for different α < 1.
<ref>(b) graphs the exponent α1-α - R in terms of α.
We see that on [0.5, 1], the maximum of this exponent is at α = 1.






§ PROOFS


A brief note about notation: We use ∼ to denote both how a random variable is sampled (ex: x ∼(0, 1) for a Gaussian x) and how a function behaves asymptotically, i.e. f(x) ∼ g(x) as x → a iff lim_x → a f(x)/g(x) = 1.
Context should be enough to differentiate between these two cases.
We in addition use ≃ to denote asymptotic expansion.
For example, if {α_i}_i ≥ 0 is a sequence of strictly decreasing reals and {β_i}_i ≥ 0 is a sequence of nonzero reals, then

    f(x) ≃∑_i ≥ 0β_i (x - ξ)^α_i

means that as x →ξ, f(x) - ∑_i = 0^N β_i (x - ξ)^α_i = Θ((x- ξ)^α_N+1).


 §.§ Preliminary Lemmas



We have

    σ_w^2  γ + σ_b^2σ_w^2   + σ_b^2 =   (1 + O( γ^-1)).

regardless of whether l = γ l / l converges.

But suppose l = γ l / l →^*.
If ^* < 1, then

    σ_w^2  γ + σ_b^2σ_w^2   + σ_b^2 =   (1 + Θ( γ^-1)).

If ^* = 1, then

    σ_w^2  γ + σ_b^2σ_w^2   + σ_b^2 =   (1 + Θ( ^-1)),

where = 1 -.


Write M = σ_b^2/σ_w^2.

    σ_w^2  γ + σ_b^2σ_w^2   + σ_b^2   =   (1 + 1 + M γ^-11 + M ^-1)
       =   (1 + M(γ - ) + O( (γ - ))).

In any situation, γ -  = O(γ) because γ≤, so this gives the first statement.
If ^* exists and ^* < 1, then γ -  = Θ(γ), which yields the second statement.
If ^* exists and ^* = 1, then γ -  = ((1-)^-1 - 1) =  ( + O(^2)) = Θ().


For any function f that is (k+1)-times differentiable in a neighborhood of 0, we have the asymptotic expansion

    f(z) = ∑_n = 0^k d^n fdz^n(0) z^nn! + O(z^k+1), as  z → 0.

Since

    .d^nd(1/q)^nq^1/2ϕ( q)|_q →∞   = (-1)^n2^n√(2π)∫_-∞^∞ϕ^2(z) z^2n z

whenever the RHS is integrable, we have

	Suppose ϕ^2(z) z^2n is integrable over z ∈ for all 0 ≤ n ≤ N + 1.
	Then ϕ( q) = q^-1/2 (∑_n = 0^N C_n q^-n + O(q^-N-1)) as q →∞, where
	
    C_n := (-1)^n2^n n! √(2π)∫_-∞^∞ϕ^2(z) z^2n z.



Note that ^d(z) = Θ(e^-d|z|) for z →∞ as long as d > 0, so that C_n from the above result converges when ϕ = ^d.
Therefore 

	Let d > 0.
	We have ^d( q) ≃ q^-1/2∑_n ≥ 0 C_n q^-n, where
	
    C_n := (-1)^n2^n n! √(2π)∫_-∞^∞^2d(z) z^2n z.


As corollaries, we obtain the following asymptotics.

	ṫȧṅḣ( q) =  2 3 √( 2 π) q^-1/2 + Θ(q^-3/2) as q →∞.


	Use <ref> along with the fact that ṫȧṅḣ(z) = ^2(z) and ∫^4 z  z =  2 3 tanh z +  1 2 ^2 z tanh z.


	1 - tanh( q) = √( 2 π) q^-1/2 + Θ(q^-3/2) as q →∞.


	Use <ref> along with the fact that 1 - tanh^2(z) = ^2(z) and ∫^2 z   z = tanh z.




^2(t) ≥exp(-t^2) for all t, with equality iff t = 0.


The lower bound is equivalent to

    2    ≥ e^t-t^2/2 + e^-t -t^2/2

The RHS has derivative (1 - t)e^t - t^2/2 - (1 + t)e^-t-t^2/2.
This is 0 iff

    1 - t1 + t = e^-2t

which has a solution 0 and in general can only have solution t∈ (-1, 1) (by considering the sign of the LHS).
Since each side is analytic in t ∈ (-1, 1), we expand

    log1 - t1 + t   = log e^-2t
    log(1-t) - log(1+t)    = -2t
    
    (-t-t^2-⋯) - (t-t^2+⋯)    = -2t
    
    -2t-2t^3-⋯   = -2t

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. 




 
Suppose ϕ = tanh.
Then ϕ̇( q) ≥ 1 √(4q+1).

As a sanity check, <ref> shows that ϕ̇( q) ∼ C_0 q^1/2 where C_0 ≈ .5319, which is above the .5 in this lemma.

	By <ref>,
	
    ϕ̇( q)    = ∫μ(z) ϕ̇^2(√(q) z)
       ≥ 1 √(2π)∫ z exp(-z^2/2 - 2qz^2)
       =  1 √(2π)∫ z exp(-(4q + 1)z^2/2)
       =  1 √(4q + 1).



<ref> demonstrates <ref>.




Let d ∈ and 1 < M < N with N - M ∈^≥ 0.
Set Σ(M, N, d) := ∑_a=M^N a^d.
If we fix M and let N →∞,

    Σ(M, N, d) = Θ(1)    if d < -1
    log N + O(1)    if d = -1
    N^d+1d+1 + O(1)    if -1 < d < 0
    
    	N - M + 1    if d =0
     1 d+1N^d+1 +  1 2 N^d + O(N^max(0, d-1))    if d > 0



	Consider the integrals A = ∫_M^N+1 a^d  a and B = ∫^N_M-1 a^d  a.
	They evaluate to A =  1 d+1((N+1)^d+1 - M^d+1) and B =  1 d+1(N^d+1 - (M-1)^d+1) when d ≠ -1 and to A = log(N+1) - log M and B = log N - log(M-1) when d = -1.
	When d ≤ 0, we have A ≤ B and Σ(M, N, d) ∈ [A, B]; when d > 0, B ≤ A and Σ(M, N, d) ∈ [B, A].
	Thus, as N →∞ with M fixed, when d < -1, Σ(M, N, d) = Θ(1);
		when d = -1, Σ(M, N, -1) = log N + O(1);
		and when d > -1, we have Σ(M, N, d) = N^d+1d+1 + O(N^d).
	
	Now for a > 0 and d > -1 and d ≠ 0, 1,
	
    ∫_a^a+1 z^d - a^d  z    =  1 d+1 ( (a+1)^d+1 - a^d )
       = (a^d +  d 2 a^d-1 + ⋯ ) - a^d
       =  d 2 a^d-1 + Θ(a^d-2).

	where the hidden constants in Θ depend only on d (and in fact this term vanishes if d = 1).
	Thus
	
    Σ(M, N, d)    = ∫_M^N+1 z^d  z - ∑_a=M^N [  d 2 a^d-1 + Θ(a^d-2)]
       =  1 d+1 ((N+1)^d+1 - M^d+1) -  d 2 Σ(M, N, d-1) + Θ(Σ(M, N, d-2))

	If -1 < d < 0, then Σ(M, N, d - 1) = Θ(1), so that Σ(M, N, d) = (N+1)^d+1d+1 + O(1) = N^d+1d+1 + O(1).
	If d > 0 and d ≠ 1, then Σ(M, N, d - 1) = N^dd, so that 
	
    Σ(M, N, d)    =  1 d+1 N^d+1 + N^d + Θ(N^max(0, d-1)) -  1 2 N^d + Θ(Σ(M, N, d-2)) 
       =  1 d+1N^d+1 +  1 2 N^d + O(N^max(0, d-1)).


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.

	Let F: ×→ be a function such that for a subset U, and for all z, z' ∈ U, z ≥ z'  F(z, n) ≥ F(z', n) for every n.
	Suppose sequences a l,  b l,  l satisfy 
	
		
  * a l+1 = F( a l, l) for all l;
		
  * b l+1≤ F( b l, l) for all l above a constant K_b.
		
  * l+1≥ F( l, l) for all l above a constant K_c.
	
	and furthermore, a l,  b l,  l all fall into U for l above a constant K_U.


	If for some m ≥max(K_b, K_U), b m ≤ a m, then b l ≤ a l, ∀ l ≥ m.
	Similarly, if for some n ≥max(K_c, K_U), n ≥ a n, then l ≥ a l, ∀ l ≥ n.


	For the first claim: b m ≤ a m  b m+1≤ F( b m, m) ≤ F( a m, m) =  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 m ≥ a m m+1≥ F( m, m) ≥ F( a m, m) =  a m+1.



 
	Suppose l satisfies the recurrence
	
    l = l-1 (1 + δl^β).

	for some nonzero constant δ∈ independent of l.
	
		
  * If β > 1, then l = Θ(1).
		
  * If β = 1, then l = Θ(l^δ).
		
  * If 0 < β < 1, then l = exp(δ1 - βl^1-β + Θ̃(l^ψ_1(1-2β))), where ψ_1(x) = max(0, x) is the ReLU function.
	


	We have
	
    log l    = logl-1 + log(1 + δ/l^β)
       = logl-1 + δ/l^β + Θ(δ^2/l^2β)

	for large l.
	If β > 1, then ∑_l l^-β converges, and
	
    log l    = log 0 - Θ(1)
     l    = Θ(1).

	If β = 1, then
	
    log l    = log 0 + δlog l + Θ(1)
     l    = Θ(l^δ).

	If β < 1, then
	
    log l    = log 0 + δ1 - β l^1 - β + Θ̃(l^1 - 2β)
     l    = exp( δ1 - β l^1 - β + Θ̃(l^ψ_1(1 - 2β))).

	


 
	Suppose l = C l^-α + l-1 (1 + δ/l^β) for α∈, C≠0, and δ≠ 0.
	Then
	
		
  * If β > 1, then
		
			
  * l = Θ(l^1 - α) if α∈ (0, 1);
			
  * l = Θ(log l) if α = 1;
			
  * l = Θ(1) if α > 1.
		
		
  * If β = 1, then
		
			
  * l = Θ(l^max(δ, 1-α)) if 1-δ≠α.
			
  * l = Θ(l^δlog l) if 1-δ = α.
		
		
		
		
		
	
	Furthermore, for β = -δ = 1, l ∼ l^-1 if α > 2, l ∼ l^1-α if α < 2, and l ∼ l^δlog l if α = 2.


	We can unwind the recurrence to get
	
    l    = ∑_m=1^l m^-α∏_n=m+1^l (1 + δn^β) +  0 ∏_n=1^l (1 + δn^β)

	
	Suppose β > 1.
	By <ref>, we get
	
    l    = Θ(1)∑_m=1^l m^-α +  0 Θ(1)
       = Θ(l^1 - α)    if α∈ (0, 1)
    Θ(log l)    if α = 1
    Θ(1)    if α > 1.

	
	Now suppose β = 1.
	By <ref>, we get
	
    l    = ∑_m=1^l m^-αΘ(m^-δl^δ) +  0 Θ(l^δ)

	where the constants hidden inside the Θ are the same in every term of the sum.
	If α > 1 - δ, then m^-δ - α = o(m^-1), so that ∑_m=1^l m^-δ - α = Θ(1), and
	
    l    = Θ(l^δ) +  0 Θ(l^δ)
       = Θ(l^δ).

	On the other hand, if α < 1 - δ, then ∑_m=1^l m^-δ - α = Θ(l^1 - δ - α).
	So
	
    l    = Θ(l^1 - α) +  0 Θ(l^δ)
       = Θ(l^1 - α).

	If α = 1 - δ, then ∑_m=1^l m^-δ - α = Θ(log l).
	So
	
    l    = Θ(l^δlog l) +  0 Θ(l^δ)
       = Θ(l^δlog l).

	
	Finally, if β∈ (0, 1), then 
	
    l    = e^δ1-βl^1-β + Θ(l^1-2β)∑_m=1^l m^-α e^-δ1-β m^1-β + Θ(m^1-2β) + e^δ1-βl^1-β + Θ(l^1-2β)

	The case of δ = -1 telescopes, so that the upper and lower constants hidden in Θ can both be taken to be 1.



 
	Suppose for some β > 0, a sequence l satisfies
	
    l = l-1 (1 - μ (l-1)^β/l),   0 ∈ (0,  1 μ).

	Then l ∼ (βμlog l)^-1/β.


	
	
	
	
	
	Consider the differential equation
	
    ẋ_μ = - μ x^β+1_μ/t

	for constant μ has solution x_μ = [β(μlog t + C)]^-1/β for some constant C determined by initial condition.
	Note that
	
    -μ x_μ(t)^β+1/t ≤ x_μ(t+1) - x_μ(t) ≤ - μ x_μ(t+1)^β+1/(t+1) = - (1 - o(t^-1))μ x_μ (t)^β+1/t.

	For any small enough α > 0, we apply <ref> with F(, l) =  - μ^β+1/l (which is monotonic in  for small enough ), l = x_μ(l), and b l = x_μ-α(l) to obtain
	
	
    x_μ - α(l) ≤ l ≤ x_μ(l)

	for large enough l and appropriately chosen initial conditions.
	This shows that l = Θ(log l ^-1/β)

	Taking α→ 0, we also obtain the leading coefficient l ∼ [βμlog l]^-1/β.
	




	Suppose a sequence u l is governed by the equation
	
    u l -  u l-1 = A( u l-1 + B)^α,

	where α∈ [0, 1) and A > 0.
	Then u l = K_1 l^ - K_2 l^log l + o(l^log l), where K_1 = [A(1-α)]^ and K_2 =  1 2 A^ (1-α)^ - 1α.


	Leading term.
	The differential equation
	
    ẋ_A, B = A(x_A, B + B)^α

	has solution x_A,B(l) = [A(1-α)(l + S)]^ 1 1-α - B for some constant S.
	Since ẋ_A, B is monotonic, we have (writing x = x_A, B for brevity)
	
    A(x_A, B(l) + B)^α = ẋ_A, B(l) ≤ x_A, B(l+1) - x_A, B(l) ≤ẋ_A, B(l+1) ≤ (A + o(1))(x_A, B(l) + B)^α

	for large enough l.
	We apply <ref> with F(x, l) = x + A(x + B)^α (which is monotonic in x for large x), l = x_A, B(l), and b l = x_A-, B(l) to obtain
	
    x_A-,B(l) ≤ u l ≤ x_A, B(l)

	for large enough l and appropriate initial conditions.
	Therefore lim u l / l^ 1 1-α∈ [[(A - ) ( 1- α)]^ 1 1 - α, [A ( 1- α)]^ 1 1 - α].
	Taking → 0 gives the leading term.
	
	Subleading term.
	Now let v l :=  u l - ℵ l^ 1 1-a, where ℵ = [A(1-α)]^ 1 1-α.
	Then we have the recurrence
	
    v l+1 + ℵ (l + 1)^ -  v l - ℵ l^   = A( v l + ℵ l^ + B)^α
     v l+1 -  v l + ℵ( l^   +  1 2 () () l^-1 + Θ(l^ - 2))  
       = A[ℵ^α l^ + α ( v l+B) ℵ^α - 1 l + Θ(( v l+B) l^-1-)]
     v l+1 -  v l    =  v l  l -  1 2 ℵ()() l^ - 1 + g(l)

	for some g(l) = O(l^ - 2 +  l) and where, to get the last equation, we have used A α^α = ℵ to cancel the l^ term and simplified α A ℵ^α-1 =.
	
	For any J > 0, the differential equation v̇_J(l) =  v_J(l)  l - J l^ - 1 has solution v_J(l) = C[l(1-α)]^ - J l^log l.
	Note that the functions F_J(z, n) = z +  z  n - J n^ - 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 v̇_J(l) and F_J(v_J(l), l) = v_J(l) + v̇_J(l) decreasing in l.
	Thus for any > 0 and l large enough
	
    G_J+(v_J(l), l) ≤ F_J+/2(v_J(l), l) ≤ v_J(l) + v̇_J(l+1) ≤ v_J(l+1) ≤ F_J(v_J(l), l) ≤ G_J-(v_J(l), l).

	Now apply <ref> with F = G_K, a l =  v l,  l = v_K-,  b l = v_K+ where K :=  1 2 ℵ()() =  1 2 A^ (1-α)^ - 1α, with appropriately chosen initial conditions.
	This yields lim_l→∞ v l / (l^log l) ∈ [-K - , -K + ] for every > 0, and there it must be equal to K.
	We have thus obtained the asymptotic expansion
	
    u l = [A(1-α) l]^ -  1 2 A^ (1-α)^ - 1α l^log l + o(l^log l).

	



	Suppose a sequence u l is governed by the equation
	
    u l -  u l-1 = -A( u l-1 + B)^α,

	where α > 1 and A > 0.
	Then u l ∼ [A(α-1)l]^ 1 1-α.


	Similar to <ref>.



	Suppose a sequence u l is governed by the equation
	
    u l -  u l-1 = A( u l-1 + B)^α + C,

	where α∈ (0, 1).
	Then u l = K_1 l^ + R(l), where the remainder R(l) is
	
    R(l) ∼
    	-K_2 l^log l	   if α >  1 2
    
    	(C-K_2) l log l			   if α =  1 2 and K_2 ≠C
    C(1-α)1-2α l	   if α <  1 2

	where K_1 = [A(1-α)]^ 1 1-α, K_2 =  1 2 A^ (1-α)^ - 1α as in <ref>.


	u is bounded below by the dynamics v l -  v l-1 = A( v l-1 + B)^α and bounded above by the dynamics w l -  w l-1 = (A + o(1))( w l-1 + B)^α.
	By <ref>, both v and w are asymptotic to u l ∼ [A(1-α)l]^ 1 1-α, which gives the result.
	
	Now define v l =  u l - [A(1-α)l]^, and similar to the proof of <ref>, we find
	
    v l+1 -  v l =  v l  l - K l^ - 1 + C + g(l)

	where K =  1 2 A^ (1-α)^ - 1α and g(l) = O(l^ - 2 +  l).
	If > 1 α >  1 2, then C + g(l) = o(l^ - 1) and we can proceed as in the proof of <ref> to find v l ∼ K l^log l.
	If = 1 α = 1 and K ≠ C, then v l+1 -  v l =  v l l^-1 - (K-C) l^ - 1 + g(l), so that the technique used in <ref> would obtain v l ∼ (K-C) l^log l = (K-C) l log l.
	If < 1 α <  1 2, then v l+1 -  v l =  v l l^-1 + C + o(1), then by using the differential equation v̇_J(l) =  v_J(l) l^-1 + J to approximate the difference equation solution and applying <ref> as in the proof of <ref>, we obtain v l(l) ∼C(1-α)1-2α l.




 §.§ Forward Dynamical Equations

Here we derive the recurrences governing the forward length and correlation quantities , , λ, γ.
We start with reduced residual networks.
*

	We have
	
    =  h_j^2  = ∑_i (w_ji x_i + b_j)^2 
       =  b_j^2  + ∑_i  w_ji^2  x_i^2 + 2∑_i  w_ji x_i b_j  + 2 ∑_j ≠ l w_ji w_li x_i^2

	But w_ji, w_li,  x, and b_j form an independency, so the last two sums are 0, and the terms in the first sum split multiplicatively.
	Therefore
	
    = σ_b^2 + ∑_i  w_ji^2  x_i^2 
       = σ_b^2 + N ·σ_w^2N
       = σ_b^2 + σ_w^2 .

	
	For the recurrence of , we have
	
    =  x_i^2  =  (ϕ(h_i) +  x_i)^2 
       = ϕ(h_i)^2 +  x_i^2  + 2 ϕ(h_i)  x_i

	As N →∞, the coefficient w_ii of x_i in h_i has vanishing covariance, so h_i and x_i become independent.
	Therefore ϕ(h_i)  x_i  = ϕ(h_i)  x_i.
	Because h_i is the sum of a large number of independent random variables, by CLT, h_i is a Gaussian with mean ∑_i  w_ji x_i +  b_j = 0 since w_ji =  b_j  = 0.
	Our antisymmetry assumption on ϕ then implies ϕ(h_i)  = 0.
	Therefore,
	
    = ϕ(h_i)^2 +  x_i^2 
       = ϕ( ) +

	as desired.


*

	Similar to <ref>.


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


	
    =   h_j^2 
    	=  ( w_j^i  x_i +  b_j)^2
    	=  ( w_j^i  x_i)^2   +   b_j^2
       = σ^2_w  x_i^2  + σ_b^2
       = σ^2_w  + σ^2_b
       =   x_i^2 
    	=  (v^j_iϕ( h_j) + x_i + a_i)^2 
       = σ_v^2ϕ( h_i)^2  +  x_i^2  + σ_a^2
       =σ_v^2 ϕ(  ) + σ_a^2 +

	where in the third equality for , 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 ϕ being antisymmetric.


*

	Similar to <ref>.




 §.§ Backward Dynamical Equations

Here we derive the recurrences governing the gradient quantities ℸ and χ_∙ for different ∙, all under the gradient independence assumption.
Write ℶ l_i =  E  x l_i for a cost function E.
*

	For a reduced residual network, we have the following derivative computation:
	
    x_i x_j   = δ_ji + ϕ̇(h_i) h_i x_j,   x_ih_j   = δ_jiϕ̇(h_j),   h_i x_j   = w_ij,   h_iw_ij   =  x_j,   h_ib_j   = δ_ij.

	
	Then
	
    ℶ_j    = ℶ_j + ∑_i ℶ_i ϕ̇(h_i) h_i x_j
       = ℶ_j + ∑_i ℶ_iϕ̇(h_i) w_ij
    ℶ_j^2    =  [ℶ_j + ∑_i ℶ_iϕ̇(h_i) w_ij]^2 
       = ℶ_j^2  + ∑_i ℶ_i^2 ϕ̇^2(h_i) (w_ij)^2
       =+ 2∑_i < kℶ_i ℶ_kϕ̇(h_i)w_ijϕ̇(h_k)w_kj + 2 ∑_iℶ_j ℶ_i ϕ̇(h_i)w_ij

	The last two terms of the above vanish as w_ij is independent from w_kj, h_i, h_k and ℶ_i, ℶ_j, ℶ_k by <ref>, and w_ij = 0.
	
	Therefore, applying <ref>,
	
    ℶ_j^2    =  σ_w^2ℶ_j^2ϕ̇^2(h_i) + ℶ_j^2
       = (σ_w^2  ϕ̇( ) + 1)ℶ_j^2

	
	We similarly have
	
    E b_j   = ∑_i  E x_ix_ih_j
    	= ℶ_j ϕ̇(h_j),    since x_ih_j = δ_jiϕ̇(h_j)
     E b_j^2    = ℶ_j^2 ϕ̇(h_j)^2 
    	= ℶ_j^2 ϕ̇( ),    by <ref>(b);
     E w_ji   = ∑_i  E x_ix_ih_jh_jw_ji
    	= ℶ_j ϕ̇(h_j) x_i,    since x_ih_j = δ_jiϕ̇(h_j)
     E w_ji^2    = ℶ_j^2 ϕ̇^2(h_j)  x_i^2 
    	= ℶ_j^2 ϕ̇( )  ,    by <ref>(b)

	In the last equation we have also used the fact that as N →∞, h_j and x_i become independent (they are jointly Gaussian and their correlation w_ji^2 goes to 0 with N).



*

	For the full residual network, we have the following derivative computations:
	
    x_i x_j   = δ_ji + ∑_k v_ikϕ̇(h_k) h_k x_j,   x_ih_j   = v_ijϕ̇(h_j),   h_i x_j   = w_ij,   h_iw_ij   =  x_j,   h_ib_i   = 1,
          x_iv_ik   = ϕ(h_k),   x_ia_i   = 1.

	
	Again let ℶ_j =  E x_j.
	Then
	
    ℶ_j    = ∑_i ℶ_i (δ_ji + ∑_k v_ikϕ̇(h_k) h_k x_j)
       = ∑_i ℶ_i (δ_ji + ∑_k v_ikϕ̇(h_k) w_kj)

	
	Thus,
	
    ℶ_j^2    =  [∑_i ℶ_i (δ_ji + ∑_k v_ikϕ̇(h_k) w_kj)]^2 
       = ℶ_j^2  + ∑_i, k v_ik^2  w_kj^2ϕ̇( ) ℶ_i^2 
       = ℶ_j^2  (1 + σ_v^2 σ_w^2 ϕ̇( ))

	where in the second equality we applied the independence argument as in the proof of <ref>, leveraging <ref>, and in the third equality we used <ref> to get ℶ_i^2 = ℶ_j^2.
	
	The other computations are similar to the proof of <ref>.





 §.§ Tanh: Reduced Residual Network



  §.§.§ Forward Dynamics

*

The case with σ_w = 0 is trivial.
We assume σ_w > 0 from here on.

 and  are asymptotically linear with l.
We first show that, for any ω < 1, 

    l +  0 ≥ l ≥ω l

and

    σ_w^2 (l +  0) + σ_b^2 ≥ l ≥σ_w^2 ω (l-1) + σ_b^2,

so that l ∼ l and l ∼σ_w^2 l.

The upper bounds are trivial, given ϕ( ) ≤ 1 for any .
We show the lower bounds for any ω < 1.

For any > 0, define ℵ_ by ϕ^2(ℵ_) = exp(-).
Then

    ϕ( )    ≥exp(-) [z ∉[-ℵ_, ℵ_]: z ∼(0, )]
       ≥exp(-) 1 - 2 ℵ_√(2π)

where the second inequality follows from an overestimate of the [z ∈ [-ℵ_, ℵ_]] via the mode of (0, ).

For any ≥ 0, ϕ( ) is then lower bounded by 

    ϕ^2√( 0)1 - 2√( 0)√(2 π 0)
    = ϕ^2√( 0)1 - √( 2 π) > 0.

Thus l and l are unbounded with l.

Furthermore, as →∞, the lower bound exp(-) 1 - 2 ℵ_√(2π) goes to exp(-), for any .
Therefore, for any ω < 1,  l ≥ω l and l ≥σ_w^2 ω (l-1) + σ_b^2.

Asymptotic expansion.
Now we repeat the following to get each successive asymptotic term of l and l:
We plug in the current asymptotic form of l into tanh( ) = 1 - C ^-1/2 + Θ(^-3/2) (<ref>), where C = √(2/π).
Next we take the sum l = ∑_r=1^l tanh(  r), which yields one more term in the asymptotic expansion of  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 l, we implicitly apply 

    = σ_w^2 l + o(l) = σ_w^2 l (1 + o(1))
    ^-1/2   = σ_w^-1 l^-1/2(1 + o(1)) = σ_w^-1 l^-1/2 + o(l^-1/2)
       = ∑_r=1^l 1 - C ( r)^-1/2 + Θ(( r)^-3/2) 
       = ∑_r=1^l 1 - C (σ_w^-1 r^-1/2 + o(r^-1/2)) + Θ(r^-3/2) 
       = l - 2C σ_w^-1 l^1/2 + o(l^1/2)
       = σ_w^2 l - 2 C σ_w l^1/2 + o(l^1/2) = σ_w^2 l (1 - 2 C σ_2^-1 l^-1/2 + o(l^-1/2))
    ^-1/2   = σ_w^-1 l^-1/2 ( 1 + C σ_w^-1 l^-1/2 + o(l^-1/2))= σ_w^-1 l^-1/2 + C σ_w^-2 l^-1 + o(l^-1)
       = ∑_r=1^l 1 - C(σ_w^-1 l^-1/2 + C σ_w^-2 l^-1 + o(l^-1)) + Θ(l^-3/2)
       = l - 2 C σ_w^-1 l^1/2 - C^2 σ_w^-2log l + o(log l)
       = σ_w^2 l (1 - 2 C σ_w^-1 l^-1/2 - C^2 σ_w^-2log ll + o(log ll))
    ^-1/2   =σ_w^-1 l^-1/2(1 + C σ_w^-1 l ^-1/2 +  1 2 C^2 σ_w^-2log ll + o (log ll))
       = ∑_r=1^l 1 - C(σ_w^-1 r^-1/2 + C σ_w^-2 r ^-1 +  1 2 C^2 σ_w^-3log rr^3/2 + o (log rr^3/2))) + Θ(r^-3/2)
       = l - 2 C σ_w^-1 l^1/2 - C^2 σ_w^-2log l + O(1)

which is what we want.





	
	Let ϕ is antisymmetric.
	Then for τ∈ [0, π/2],
	
    ϕ( , cosτ)    = lim_t →τ 1 πsin t∫_w' ≥ |w| w  w' Υ(w, w';τ) ϕ(√()√(2)(w + w')) ϕ(√()√(2)(w' - w))
       =  1 π∫_0^∞ r r e^-r^2/2∫_0^πθΣ(√()r, θ; τ)
       =  1 π∫_0^∞ s  s  e^-s^2  /2∫_0^πθΣ(s, θ; τ)
       =  1 π∫_0^πθ∫_0^∞ s e^-s^2  /2sΣ(s, θ; τ)

	where Υ(w, w';τ) :=e^- 1 2(w^21-c + (w')^21+c) - e^- 1 2((w')^21-c + w^21+c) with c = cosτ, and Σ(s, θ; τ) := ϕ(s sinθ) ϕ(s sin(θ - τ)).

Of course, in the above lemma, the limit in the first equation is only necessary when τ = 0 or τ = π/2.

	
	Let c := cosτ and
	
    Γ := ϕ( , cq) =  1 2π√(1 - c^2)∫𝐳exp(-𝐳^T Σ^-1𝐳 /2) ϕ(z)ϕ(z'),

	where Σ = [    cq; cq    ].
	
	Our proof will have two portions: Symmetrization of the Γ integral and trigonometric change of variables for evaluation.
	
	Symmetrization. Σ is diagonalized by Ω =  1 √(2q)[ -1  1;  1  1 ],
	
    Σ = Ω^T (1 - c, 1 + c) Ω.

	By a change of variable 𝐰 = Ω𝐳, so that 𝐰 = 𝐳, we have
	
    Γ   =  1 2π√(1 - c^2)∫𝐰exp(-𝐰^T (1 - c, 1 + c)^-1𝐰/2) ϕ(√()√(2)(w' - w)) ϕ(√()√(2)(w + w'))
       =  1 2π√(1 - c^2)∫ w  w' e^- 1 2(w^21-c + (w')^21+c)ϕ(√()√(2)(w' - w)) ϕ(√()√(2)(w + w'))

	
	By a change of variable swapping w with w', we get
	
	
    Γ   = - 1 2π√(1 - c^2)∫ w  w' e^- 1 2((w')^21-c + w^21+c)ϕ(√()√(2)(w + w')) ϕ(√()√(2)(w' - w))

	
	Thus
	
    2 Γ   =  1 2π√(1 - c^2)∫ w  w' Υ(w, w';τ) ϕ(√()√(2)(w + w')) ϕ(√()√(2)(w' - w))

	where
	
    Υ(w, w';τ) = e^- 1 2(w^21-c + (w')^21+c) - e^- 1 2((w')^21-c + w^21+c).

	
	
	Note that, by the antisymmetry of ϕ, the integrand K := Υ(w, w';τ)ϕ(…)ϕ(…) above has the symmetries K(w, w') = K(w', w) = K(w, -w'), and is everywhere nonnegative.
	<ref> displays a contour plot of K for typical values of  and c.
	So
	
    Γ   =  1 π√(1 - c^2)∫_w' ≥ |w| w  w' K(w, w').

	This gives the first equation in the lemma.
	
	Polar Coordinates.
	Let w √(1 - c) = r cosθ, w'√(1 + c) = r sinθ, so that
	
    w    = r cosθ√(1-c) = √(2) r cosθsinτ 2
    
    	w'    = r sinθ√(1+c) = √(2) r sinθcosτ 2
     w  w'    = √(1-c^2)r  r θ = (sin^2 τ) r  r θ.

	Then
	
    𝐀   := ∫_w' ≥ |w| e^-(w^21-c + (w')^21 +c)/2ϕ(√(/2)(w+w'))ϕ(√(/2)(w'-w))  w  w'
       = sin^2 τ∫_0^∞ r r e^-r^2/2∫_τ/2^π - τ/2θϕ(√()r sin(θ + τ/2))ϕ(√()r sin(θ - τ/2)).

	Similarly, let w √(1 + c) = r cosθ, w'√(1 - c) = r sinθ, so that
	
    w    = r cosθ√(1+c) = √(2) r cosθcosτ 2
    
    	w'    = r sinθ√(1-c) = √(2) r sinθsinτ 2
     w  w'    = √(1-c^2)r  r θ = (sin^2 τ) r  r θ,

	and
	
    𝐁   = ∫_w' ≥ |w| e^-(w^21+c + (w')^21 - c)/2ϕ(√(/2)(w+w'))ϕ(√(/2)(w'-w))  w  w'
       = -sin^2 τ∫_0^∞ r r e^-r^2/2∫_π/2 - τ/2^π/2 + τ/2θϕ(√()r cos(θ + τ/2))ϕ(√()r cos(θ - τ/2))
       = -sin^2 τ∫_0^∞ r r e^-r^2/2∫_- τ/2^τ/2θϕ(√()r sin(θ + τ/2))ϕ(√()r sin(θ - τ/2)).

	Thus
	
    Γ   =  1 π√(1-c^2) (𝐀 - 𝐁)
       =  1 π∫_0^∞ r r e^-r^2/2∫_- τ/2^π - τ/2θϕ(√()r sin(θ + τ/2))ϕ(√()r sin(θ - τ/2))
       =  1 π∫_0^∞ r r e^-r^2/2∫_0^πθϕ(√()r sin(θ))ϕ(√()r sin(θ - τ)).

	This gives the second equation in the lemma, and a change of variables s = √()r gives the third.
	
	For the fourth equality, we start from the third equality, and apply integration by parts:
	
    = 1 π∫_0^∞ s  s  e^-s^2  /2∫_0^πθΣ(s, θ; τ)
       =  1 π∫_0^πθ∫_0^∞ s s  e^-s^2  /2Σ(s, θ; τ)
       =  1 π∫_0^πθ(-e^-s^2  /2Σ(s, θ; τ)|_s=0^∞ + ∫_0^∞ s e^-s^2  /2sΣ(s, θ; τ) )
       =  1 π∫_0^πθ∫_0^∞ s e^-s^2  /2sΣ(s, θ; τ).

	where the last equality follows because Σ(0, θ; τ) = 0 and e^-s^2 / 2→ 0 as s →∞.


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

	Suppose a function f: (0, 2) → is C^k on (0,2).
	If lim_x ↓ 0 f i (x) exists and is finite for every i ∈ [0, k], then f can be extended to [0, 2) such that one sided ith derivatives exist at 0 for all i ∈ [0, k].


	Consider f i(0) :=  f i (1) - ∫_0^1  f i+1(x)  x for i ∈ [0, k-1], which naturally is also equal to f i () - ∫_0^ f i+1(x) x for any > 0.
	Certainly f i (x) → f i(0) as x → 0 if this limit exists — and by assumption it does, for 0 ≤ i ≤ k -1.
	Therefore, we can define the extension of f i to x = 0 to be f i (0) :=  f i(0).
	But we need to check that for i ∈ [0, k-1].
	
    lim_→ 0 1 ( f i() -  f i (0)) =  f i+1(0)

	so that all one sided ith derivatives exist.
	But
	
    1 ( f i() -  f i (0))    =  1 ∫_0^ f i+1(x) x
       = f i+1(0) + ∫_0^1 ( f i+1(x) -  f i (0)) (x ∈ [0, ])  x

	Since lim_x ↓ 0 f i+1(x) =  f i+1 (0), f i+1(x) -  f i+1(0) is bounded for small x, and by dominated convergence, ∫_0^1 ( f i+1(x) -  f i (0)) (x ∈ [0, ])  x →∫_0^1 0  x = 0 as → 0.
	Thus
	
    lim_→ 0 1 ( f i() -  f i (0))    =  f i+1(0)

	as desired.



	If f: [0, 2) → is C^k on (0, 2) and has one sided derivatives at 0 up to order k, then
	
    f()    = f(0) +  f 1 (0) + ⋯ + ^i-1(i-1)! f i-1 (0) + O(^i)

	for any i ≤ k.


	We have
	
    f()    = f(0) + ∫_0^ f 1 (x)  x
       = f(0) +  f 1(0) + ∫_0^ f 1(x) -  f 1 (0)  x
       = f(0) +  f 1(0) + ∫_0^∫_0^x_0 f 2 (x_2)  x_2  x_1
       = f(0) +  f 1(0) + ^22 f 2(0) + ∫_0^∫_0^x_1 f 2 (x_2) -  f 2 (0)  x_2  x_1
         ⋮ 
    
    	f()    = f(0) +  f 1 (0) + ⋯ + ^i-1(i-1)! f i-1 (0) + ∫_0^ x_1 ∫_0^x_1 x_2 ⋯∫_0^x_i-1 x_i f i(x_i)

	for any i≤ k.
	It suffices then to bound the size of the integral.
	Since f i (x) → f i (0) as x ↓ 0 by assumption, | f i (x_i)| is bounded by some constant C on the integration region 𝔸 := {(x_1, …, x_i): ≥ x_1 ≥⋯≥ x_i} for small enough .
	Therefore,
	
    =∫_0^ x_1 ∫_0^x_1 x_2 ⋯∫_0^x_i-1 x_i f i(x_i)
       = ∫ f i(x_i) (x⃗∈𝔸) x⃗
       ≤ C |𝔸|
       = Θ(^i).

		


As a corollary,

	If f: (0, 2) → is smooth on (0, 2) and lim_x → 0 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 
	
    f()    = f(0) +  f 1 (0) + ⋯ + ^i-1(i-1)! f i-1 (0) + O(^i)

	for any i.



	
	Let ϕ = tanh.
	For any fixed c, ϕ( , cq) is smooth (infinitely differentiable) on ∈ (0, ∞).
	As a function of Q :=, it can be extended smoothly to the point Q = 0, so that
	
    ϕ( , cq )    = lim_' →∞ϕ( ', cq') + lim_' →∞ϕ( ', cq')/ (')^-1 + ⋯
           + ^-i+1(i-1)!lim_' →∞^i-1ϕ( ', cq')/ (')^-i+1 + O(^-i)

	for any i ≥ 0.
	Furthermore, for c bounded away from 1, the constants hidden O can be taken independent of c.


	
	Smoothness on (0, ∞).
	By the third equation of <ref>, for Q ∈ (0, ∞) ∈ (0, ∞),
	
    = 1 π∫_0^∞ s  s |^nQ^n Q e^-s^2 Q/2| ∫_0^πθ |ϕ(s sinθ)ϕ(s sin(θ - τ))|
       ≤∫_0^∞ s  s |^nQ^n Q e^-s^2 Q/2| < ∞,

	so by Leibniz's integral rule and a simple induction, all derivatives of ϕ( , cq) against Q exists for any Q ∈ (0, ∞).
	
	Extension to Q = 0.
	By <ref>, it suffices to show that the limit of ^kϕ( , cq)Q^k exists and is finite as Q → 0, for all k.
	Let τ = arccos c.
	By the fourth equation of <ref>, we have explicitly
	
    ^kϕ( , cq)Q^k   =  1 π∫_0^πθ∫_0^∞ s (-s^2/2)^k e^-s^2 Q /2sΣ(s, θ; τ)
       = (-2)^-kπ∫_0^πθ∫_0^∞ s s^2k e^-s^2 Q /2sΣ(s, θ; τ)

	for any Q ∈ (0, ∞).
	Note that for ϕ = tanh, ϕ̇= ^2,
	
    sΣ(s, θ; τ)
    		   = sinθϕ̇(s sinθ)ϕ(s sin (θ - τ)) + sin(θ-τ)ϕ(s sinθ)ϕ̇(s sin (θ - τ)).

	We split the integral of ^k ϕQ^k as follows:
	
    ^kϕ( , cq)Q^k   = (-2)^-kπ∫_0^πθ∫_0^∞ s s^2k e^-s^2 Q /2sinθϕ̇(s sinθ)ϕ(s sin (θ - τ))
            + (-2)^-kπ∫_0^πθ∫_0^∞ s s^2k e^-s^2 Q /2sin(θ-τ)ϕ(s sinθ)ϕ̇(s sin (θ - τ))

	We show that for each piece, the limit as Q → 0 exists and is finite, for any k.
	This will prove the smooth extendability of ϕ 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.
	
    =∫_0^πθ∫_0^∞ s s^2k e^-s^2 Q /2sinθϕ̇(s sinθ)ϕ(s sin (θ - τ))
       = ∫_0^∞ s s^2k e^-s^2 Q/2∫_0^πθsinθϕ̇(s sinθ) ϕ(s sin (θ - τ))

	We now try to bound the inner integral by an exponentially decreasing term e^-s μ for some μ; clearly, by monotone convergence on the outer integral as Q → 0, this would show the limit of the integral exists and is finite.
	
	Because ϕ is odd and ϕ̇ is even, the inner integrand is negative on θ∈ [0, τ) and positive on θ∈ (τ, π].
	We will break up the inner integral as follows, for some fixed > 0 satisfying τ -  > 0 independent of s (recall τ∈ (0, π/2]).
	
    =∫_0^πθsinθϕ̇(s sinθ) ϕ(s sin (θ - τ))
       = (∫_0^ + ∫_π - ^π) θsinθϕ̇(s sinθ) ϕ(s sin (θ - τ)) + ∫_ ^π - θsinθϕ̇(s sinθ) ϕ(s sin (θ - τ))

	Now because ϕ̇(z) = ^2(z) ≤ 2 e^-z, and sinθ≥sin on θ∈ [, π - ],
	
    =|∫_ ^π - θsinθϕ̇(s sinθ) ϕ(s sin (θ - τ))|
       ≤ 2∫_ ^π - θexp(-s sin) 
       = 2(π - 2) exp(-s sin).

	
	For the other part:
	
    =∫_π - ^πθsinθϕ̇(s sinθ) ϕ(s sin (θ - τ))
       = ∫_^0sin(π - θ) ϕ̇(s sinπ - θ) ϕ(s sin (π - θ - τ)) (π - θ)
       = ∫_0^θsinθϕ̇(s sinθ) ϕ(s sinθ + τ)

	so that
	
    =(∫_0^ + ∫_π - ^π) θsinθϕ̇(s sinθ) ϕ(s sin (θ - τ))
       = ∫_0^θsinθϕ̇(s sinθ) [ϕ(s sin (τ + θ)) - ϕ(s sin(τ - θ))]

	But by intermediate value theorem, ϕ(s sin (τ + θ)) - ϕ(s sin(τ - θ)) = 2 θϕ(s sin(τ + θ))/θ |_θ = ψ = 2 θϕ̇(s sin (τ + ψ))s cos(τ + ψ) for some ψ∈ [- θ, θ].
	By the assumption on , ϕ(s sin (τ + θ)) - ϕ(s sin(τ - θ)) ≤ 2ϕ̇(s sin(τ - ))s cos(τ - ).
	Then
	
    =∫_0^θsinθϕ̇(s sinθ) [ϕ(s sin (τ + θ)) - ϕ(s sin(τ - θ))]
       ≤∫_0^θsinθϕ̇(s sinθ) 2 ϕ̇(s sin(τ - ))s cos(τ - )
       ≤ 2 ϕ̇(s sin(τ - ))s cos(τ - ) O(1)

	Because τ - > 0 by assumption on , and because ϕ̇(z) = exp(-Θ_+(z)), this quantity is exp(-Θ_+(z)), as desired (here Θ_+ denotes a positive quantity).
	
	Thus
	
    =∫_0^πθsinθϕ̇(s sinθ) ϕ(s sin (θ - τ))
       = (∫_0^ + ∫_π - ^π) θsinθϕ̇(s sinθ) ϕ(s sin (θ - τ)) + ∫_ ^π - θsinθϕ̇(s sinθ) ϕ(s sin (θ - τ))
       = exp(-Θ_+(s))

	and similarly for the other piece of ^k ϕQ^k, so that
	
    =∫_0^∞ s s^2k e^-s^2 Q/2∫_0^πθsinθϕ̇(s sinθ) ϕ(s sin (θ - τ))
       = ∫_0^∞ s s^2k e^-s^2 Q 2 - Θ_+(z)
       →∫_0^∞ s s^2k e^-Θ_+(z)

	is finite as Q → 0, by monotone convergence.
	
	Independence of constant hidden in O((')^-i).
	The constant hidden is a function of the  chosen above, which depend on τ, but only to the extent that it must satisfy τ -  > 0.
	As long as we are interested in a set 𝒞 of c that is bounded away from 1, the corresponding set of τ is bounded away from 0, so  can be taken to be some number smaller than all of the corresponding τ.



 
	
	
	Suppose ϕ is tanh-like.
	Then for c ∈ [0, 1],
	
    ϕ( , cq) ≤ 2 πarcsin(c),

	and weakly increases to this upper bound as →∞.
	Furthermore, 
	
		
  * If c = 0 or 1, then equality holds regardless of .
		
  * If c ∈ (0, 1) is held constant, 2 πarcsin(c) - ϕ( , cq) = Θ(), where the hidden constants in Θ depend on c.
		But the constants can be made independent of c if c ∈ [, 1-] for some > 0.
	



	
	The cases of c = 0 or 1 are obvious by the definition of .
	So from here on we assume c ∈ (0, 1).
	
	Let τ := arccos c.
	By the first equation of <ref> and the assumption that ϕ is tanh-like, it is immediate that ϕ( ,cq) is nondecreasing in .
	By dominated convergence, using the second equation of <ref>, we get
	
    lim_→∞ϕ( ,cq)    =  1 π∫_0^∞ r r e^-r^2/2 (π - 2 τ)
       = π - 2 τπ
       =  2 πarcsin c.

	
	Then the convergence rate is O() by <ref> and Taylor's theorem.
	Thus to show the convergence rate is Θ(), it suffices to show that 𝐃 := ϕ( , cq)Q < 0.
	But this is apparent from the first equation of <ref>:
	For τ∈ (0, π/2),
	
    𝐃   =  1 πsinτ∫_w' ≥ |w| w  w' Υ(w, w';τ)(- 1 2√(2) Q^-3/2)
                     ×[ϕ̇(√(/2)(w + w'))ϕ(√(/2)(w'-w))
                        + ϕ(√(/2)(w + w'))ϕ̇(√(/2)(w'-w))]
        < 0

	since Υ is positive on the integration domain, and ϕ̇ and ϕ are both positive for positive arguments, by the assumption of ϕ being tanh-like.
	
	Independence of the constants in Θ() from c when c ∈ [, 1-].
	By <ref>, the upper constant can be made independent from c.
	Since 𝐃 is monotonically decreasing in c (or monotonically increasing in τ) and |𝐃| is monotonically increasing in c (or monotonically decreasing in τ), we have |𝐃| > |𝐃|]_c =, which can be taken to be the lower constant in Θ().

<ref> verifies empirically that the subleading term in tanh( q, cq) is linear in q, for constant c.




*


We have by <ref>,

    γ   =  2 πarcsin(λ /) - Θ(^-1) + γ.


Since = σ_w^2  + σ_b^2 by <ref>, and λ = σ_w^2 γ + σ_b^2 by <ref>,

    γ   =  2 πarcsinσ_w^2 γ + σ_b^2σ_w^2  + σ_b^2 - Θ(^-1) + γ.



We claim that γ l →∞ as l →∞.
Otherwise, there is some C such that γ l ≤ C for all l.
For large enough l, l ≥ω l for any ω < 1 and arcsinCσ_w^2 l-1 + σ_b^2 = Θ(1/l) by linearization of arcsin.
Thus γ l = Θ(log l), but this contradicts our assumption that γ is bounded.
This proves our claim.

Therefore, for large enough l,

    σ_w^2 γ + σ_b^2σ_w^2  + σ_b^2 = γ/ + Θ(l^-1).






<ref> shows 2πarcsin x vs x.
One sees that 1 is an unstable fixed point; if < 1 - ϵ, then 2 πarcsin < 1 - ϵ - δ for some δ.
Thus  drops monotonically until some threshold under which the linearization of arcsin, arcsin x = x + Θ(x^3), is applicable.
So for large enough l,

    γ - γ   =  2 πarcsin (γ /  + Θ( l)) - Θ( l)
       =  2 πγ /  + O( l)

As l ∼ l by <ref>, this difference equation has solution γ = Ω(l^ 2 π -), O(l^ 2 π + ) for any  by using the dynamics of <ref> to upper and lower bound this difference equation.






  §.§.§ Backward Dynamics

*

	The σ_w = 0 case is obvious.
	We will assume σ_w > 0 from here on.
	
	Let l = b_0 l + b_1 l^1/2 + b_2 log l + O(1).
	Then for D =  2 3 √( 2 π), we have (implicitly applying <ref> and <ref>),
	
    ^-1/2   = σ_w^-1 b_0^-1/2 l^-1/2(1 - b_1 b_0^-1 2 l^-1/2 - b_2 b_0^-1 2^-1 l^-1log l + O(l^-1))
    ϕ̇( )     = D ^-1/2 + Θ(^-3/2)
       = Dσ_w^-1 b_0^-1/2 l^-1/2(1 - b_1 b_0^-1 2 l^-1/2 - b_2 b_0^-1 2^-1 l^-1log l + O(l^-1))
    log(B ϕ̇( ) + 1)    = B D σ_w^-1 b_0^-1/2 l^-1/2
       =- (BD σ_w^-1 b_0^-3/2 b_1 2^-1 +B^2 D^2 σ_w^-2 b_0^-12^-1)l^-1 + Θ(l^-3/2log l)
    ∑_r=1^l log(B ϕ̇(  r) + 1)    = 2B D σ_w^-1 b_0^-1/2 l^1/2
       = -(BD σ_w^-1 b_0^-3/2 b_1 2^-1 +B^2 D^2 σ_w^-2 b_0^-12^-1)log l + O(1)

	In our case, we have b_0 = 1, b_1 = -2C σ_w^-1, b_2 = C^2 σ_w^-2,  B = σ_w^2, C = √( 2 π), which gives
	
    ∑_r=1^l log(B ϕ̇(  r) + 1) =  4 3 √( 2 π)σ_w l^1/2 + ( 4 3π - σ_w^2  4 9π) log l + O(1).

	so that
	
    ℸm/ℸ l    = exp[ 4 3 √( 2 π)σ_w (√(l) - √(m)) + ( 4 3π - σ_w^2  4 9π) (log l - log m) + O(1)]



*

	The σ_w = 0 case is obvious.
	We will assume σ_w > 0 from here on.
	
	As in the proof of <ref>,
	
    ϕ̇( ) = Dσ_w^-1 b_0^-1/2 l^-1/2 + Θ( l)

	where D =  2 3 √( 2 π).
	Thus by <ref>,
	
    log(ℸm/ℸ l) =  4 3 √( 2 π)σ_w (√(l) - √(m)) + ( 4 3π - σ_w^2  4 9π) (log l - log m) + O(1)

	
    log(χ m_b/χ l _b) = 4 3 √( 2 π)σ_w (√(l) - √(m)) + ( 4 3π -  1 2 - σ_w^2  4 9π) (log l - log m) + O(1)

	Similarly, since = l + Θ(√(l)) by <ref>, we have
	
    log(χ m_w/χ l _w) =  4 3 √( 2 π)σ_w (√(l) - √(m)) + ( 4 3π +  1 2 - σ_w^2  4 9π) (log l - log m) + O(1)





 §.§ Tanh: Full Residual Network



  §.§.§ Forward Dynamics



*

	The σ_w = 0 case is obvious.
	We will assume σ_w > 0 from here on.
	
	As in <ref>,  will have expansion = b_0 l + b_1 l^1/2 + b_2 log l + O(1).
	Then, for C = √( 2 π),
	
    ^-1/2   = σ_w^-1 b_0^-1/2 l^-1/2(1 - b_1 b_0^-1 2 l^-1/2 - b_2 b_0^-1 2^-1 l^-1log l + O(l^-1))
    ∑_r=1^l ϕ(  r)    = ∑_r=1^l 1 - C( r)^-1/2 + Θ(( r)^-3/2)
       = l - 2 C σ_w^-1 b_0^-1/2 l^1/2 + C σ_w^-1 b_1 b_0^-3/2 2^-1log l + O(1)
     l    = σ_v^2 ∑_r=1^l + σ_a^2 l
       = (σ_v^2 + σ_a^2) l - 2C σ_v^2 σ_w^-1 b_0^-1/2 l^1/2 + C σ_v^2 σ_w^-1 b_1 b_0^-3/2 2^-1log l + O(1)

	which yields
	
    b_0    = σ_v^2 + σ_a^2
    
    	b_1    = -2C σ_v^2 σ_w^-1 b_0^-1/2 = -2C σ_v^2 σ_w^-1√(σ_v^2 + σ_a^2)
    
    	b_2    = -C^2 σ_v^4 σ_w^-2(σ_v^2 + σ_a^2)^2






Suppose ϕ is tanh-like.
Then

    γ≤σ_v^2  2 πarcsinλ/ + σ_a^2 + γ,

and

    σ_v^2  2 πarcsinλ/ + σ_a^2 + γ - γ = Θ(^-1).




Similar to the proof of <ref>.



	Let u^* ∈ [0, 1).
	Let f_t: [0, 1) → [0, 1] be a continuous function for each t ∈, to each of which we associate two numbers 0 ≤ a_t ≤ u^* ≤ b_t ≤ 1.
	Suppose for each t, f_t(u) > u for all u ∈ [0, a_t) and f_t(u) < u for all u ∈ (b_t, 1).
	Assume that for each u, f_t(u) - u → 0 as t →∞ uniformly over u.
	If a_t ↗ u^* and b_t ↘ u^*, then for any u_0 ∈ [0, 1), the dynamics u_t = f_t(u_t-1) has a limiting point.
	Furthermore, either u_t → u^* or u_t eventually converges monotonically (decreasing or increasing) to a limit point.



	Fix a u_0 ∈ [0, 1).
	If u_t → u^* then we are done.
	Otherwise, suppose there is a neighborhood [u^* - , u^* + ] such that for an infinite sequence t_1, t_2, …, u_t_i∉[u^* - , u^* + ].
	WLOG assume u_t_i < u^* - for all i and (t_i)_i is the sequence of all ts that satisfy this inequality.
	
	If (t_i)_i contains {s: s ≥ N} for some N, then for some M > N, for every t > M, a_t > u^* -  > u_t.
	By assumption, u_t is monotonic for all t > M but is bounded above.
	Thus u_t has a fixed point û≤ u^* - as desired.
	
	Now assume there are infinite is such that t_i - 1 ≠ 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^* - /2 and |f_t(u) - u| < /4 for all u and for all t ≥ T (T exists by premise).
	Let j be the smallest index such that t_j > T and t_j - 1 ≠ t_j-1.
	By the definition of j, u_t_j - 1≥ u^* -.
	If u_t_j - 1≥ u^* - /2, then by definition of T, u^* -  > u_t_j = f_t_j(u_t_j - 1) > u_t_j-1 - /4 > u^* - 3/4 > u^* -, a contradiction.
	If u^* - ≤ u_t_j - 1≤ u^* - /2, then by the definition of T, u_t_j-1≤ a_t_j-1 so that u_t_j  = f_t_j(u_t_j - 1) > u_t_j - 1≥ u^* - , a contradiction.
	
	The “furthermore” claim is clear from our proof above.


*

	The σ_w = 0 case is obvious.
	We will assume σ_w > 0 from here on.
	

If σ_a = 0, then ^* as defined above is 0, and = γ decreases as Θ(l^ 2 π - 1) to 0, by the same reason as before.

So from now on suppose σ_a > 0.
We apply <ref> first to show that  converges.
We have

    σ_v^2 ϕ( , cq) + σ_a^2    =  - 
       =  -  +  - 
       = ( - )  + ( - )
       = ( - )[( - ) -  + ]
    σ_v^2 ϕ( , cq) + σ_a^2σ_v^2 ϕ( ) + σ_a^2   = ( - )  -  + 
    - [σ_v^2 ϕ + σ_a^2σ_v^2 ϕ + σ_a^2 - ]    =  -

If we define f_l(u) :=  l - l-1 l[σ_v^2 ϕ(  l,  l  l) + σ_a^2σ_v^2 ϕ(  l) + σ_a^2 - u] + u (the LHS of the above), then f_l(u) - u = O(l^-1) uniformly for all u because l = Θ(l), l - l-1 = Θ(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), ….
For each I_k, it is clear that the trajectories of l = f_l(l-1) with 0 ∈ I_k will fall into some interval J_k bounded away from 1 for all l ≥ L, for large enough L (dependent on k).
Then we can apply <ref> to get f_l(u) =  l - l-1 l[σ_v^2 2 πarcsin(u) + σ_a^2σ_v^2 + σ_a^2 - u + o(1)] + u where the constants in o(1) is uniform for all 0 ∈ I_k.
For u < ^* (as defined in the theorem statement), σ_v^2 2 πarcsin(u) + σ_a^2σ_v^2 + σ_a^2 > u and for u > ^*, σ_v^2 2 πarcsin(u) + σ_a^2σ_v^2 + σ_a^2 < u (see <ref>).
Thus as l →∞, the o(1) term gets smaller and smaller, and this monotonicity holds for f_l(u) - u = [σ_v^2 2 πarcsin(u) + σ_a^2σ_v^2 + σ_a^2 - u + o(1)] > 0 (resp. <0) on larger and larger intervals [0, a_l]  ∩ J_k (resp. [b_l, 1) ∩ J_k).
This proves all the preconditions for <ref>, which yields that I_k converges to a limit point.
As this argument is independent of k, we have that for all 0 ∈ [0, 1), l converges.

Now we solve for the limit point.







Suppose  has limit point ^† (possibly different from ^* described in the theorem); if we express γ l = (^† +  l )  l, then

    σ_v^2 ϕ( , cq) + σ_a^2    = γ - γ
       = (^† + ) - (^† + ) 
       = ^†( - ) +  - 
    σ_v^2 ϕ( , cq) + σ_a^2σ_v^2 ϕ( ) + σ_a^2   = ^† +  + ( - ) -

As l →∞, ∼→^†, and ϕ( , ^†) → 2 πarcsin(^†), and ϕ( ) → 1.
Additionally, /( - ) = Θ(l) and = o(1) so that -  = o(l^-1).
Then we have, taking limits l →∞,

    σ_v^2  2 πarcsin(^†)+σ_a^2σ_v^2 + σ_a^2   = ^†.

Since f_l (as defined above) repels points away from 1, the only solution for ^† when 0 < 1 is ^† = ^* as specified in the theorem statement.

We defer the proof of the convergence rate to ^* to <ref>.





 
Let ^* be the stable fixed point determined by σ_a and σ_v.
Then as long as σ_v > 0,

    2 π 1 √(1 - (^*)^2)σ_v^2 σ_v^2 + σ_a^2∈ ( 1 2,  2 π]



	
	Write ρ := σ_a^2σ_v^2.
	By definition of , we get
	
    = (1 - ρ)  2 πarcsin + ρ
    ρ =    =  -  2 πarcsin1 -

	Substituting ρ into the expression in question, it follows that we want to show
	
    2 π (1 - ^2)^-1/2 (1 + ρ)^-1=   2 π (1 - ^2)^-1/2(1 -  2 πarcsin1 - )^-1∈ ( 1 2,  2 π]

	for ∈ [0, 1) (the endpoint at 1 is not included since σ_v > 0.
	But this is
	
    = 2 π (1 - )^1/2 (1 + )^-1/2 (1 -  2 πarcsin)^-1.

	Set g() to be this expression.
	We could proceed by finding critical points, but a simple plot <ref> shows that g is decreasing on [0, 1), with extremal values at the end points:
	
    g() ∈ [lim_→ 1 g(), g(0)),   for ∈ [0, 1).

	Obviously g(0) =  2 π.
	For the limit, we note that arcsin has an asymptotic expansion π 2 - √(2) (1 - e)^1/2 + Θ((1-e)^3/2) at 1, so that (1 - )^1/2(1 -  2 πarcsin)^-1→π2 √(2), and g() →12 as → 1.
	
	
	





If 0 < 1, then | l - ^*| is Ω(l^-δ^* - ε) and O(l^-δ^* + ε) for any ε > 0, where

    δ^* := 1 -  2 π 1 √(1 - (^*)^2)σ_v^2 σ_v^2 + σ_a^2∈ [1- 2 π,  1 2),

where the bounds on the right follow from <ref>.



Define ω(q, c) =  2 πarcsin(c) - tanh( q, c q).
By <ref>, for large enough l,  is close to ^* bounded away from 0 or 1, so that ω(, ) = Θ() with the constant hidden in Θ independent of .
Additionally, by <ref>, 1 - tanh( ) = Θ(^-1/2).
Therefore,

    (^* + ϵ)   = σ_v^2 ( 2 πarcsin(^* + ϵ) - ω(, ))+ σ_a^2 + γ
       = σ_v^2  2 π [arcsin(^*) + ϵ√(1 - (^*)^2) + Θ(^2)] - Θ( l) + σ_a^2 + γ
       = ^*(σ_v^2 + σ_a^2) + (^* + )  + σ_v^2  2 πϵ√(1 - (^*)^2) + Θ(^2) - Θ( l)
    ^*( -  - σ_v^2 - σ_a^2)    =  -  + σ_v^2  2 πϵ√(1 - (^*)^2) + Θ(^2)- Θ( l)
    ^*σ_v^2(ϕ( ) - 1)    =   -  + σ_v^2  2 πϵ√(1 - (^*)^2) + Θ(^2) - Θ( l)
       =  1  (^*σ_v^2(1 - ϕ( )) + Θ(^2) - Θ( l) +  (+ σ_v^2  2 π1√(1 - (^*)^2)))
       = Θ(l^-3/2) + (1 - δ l / l)

where

    δ l    =  l  (σ_v^2 ϕ( ) + σ_a^2 -σ_v^2  2 π1√(1 - (^*)^2)) + Θ(/l)
       = (1 + Θ(l^-1/2))(σ_v^2 (1 - Θ(l^-1/2)) + σ_a^2 -σ_v^2  2 π1√(1 - (^*)^2))/(σ_v^2 + σ_a^2) + Θ( / l)
       = δ^* + O(l^-1/2),

where δ^* := 1 -  2 π 1 √(1 - (^*)^2)σ_v^2 σ_v^2 + σ_a^2, which is positive by <ref>.
By taking the δ of <ref> to be δ^* + ε or δ^* - ε respectively for lower and upper bounding the dynamics of l, the solution l is Ω(l^-δ^* - ε) and O(l^-δ^* + ε) for any ε > 0 since 1 2 > δ^*.





  §.§.§ Backward Dynamics


*

	The σ_w = 0 case is obvious.
	We will assume σ_w > 0 from here on.
	
As in the proof of <ref>,

    log(ℸm / ℸ l)    = 2B D σ_w^-1 b_0^-1/2 (√(l) - √(m))
    
       = -(BD σ_w^-1 b_0^-3/2 b_1 2^-1 +B^2 D^2 σ_w^-2 b_0^-12^-1)(log l - log m) + O(1)

where B = σ_v^2 σ_w^2, D =  2 3 √( 2 π),

    b_0    = σ_v^2 + σ_a^2
    
    b_1    = -2C σ_v^2 σ_w^-1√(σ_v^2 + σ_a^2)
    
    b_2    = -C^2 σ_v^4 σ_w^-2(σ_v^2 + σ_a^2)^2.

with C = √( 2 π).
This simplifies to the desired form.


*

	Similar to <ref>.




 §.§ α-ReLU: Full Residual Network


The following can be checked readily
*

Since ψ̇_α = αψ_α - 1, we have as a corollary,

If α >  1 2, then
ψ̇_α( q) = α^2 _α - 1 q^α-1.

As a special case, when α = 1, _α =  1 2.

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

	_α+1/_α = 2α+1.




  §.§.§ Forward Dynamics




Suppose we have the nonlinearity ϕ = ψ_1.
Then l = Θ((1 + σ_v^2 σ_w^2/2)^l), with the hidden constant depending on the initial condition.


We have

    =  1 2σ_v^2(σ_w^2  + σ_b^2) + σ_a^2 + 
       = ( 1 2σ_v^2σ_w^2 + 1)  +  1 2 (σ_v^2 σ_b^2 + σ_a^2).

By the standard method of characteristic equation, we get that

    l = A + C B^l

where A = - σ_a^2 + σ_b^2 σ_v^2σ_v^2 σ_w^2, B = 1 + σ_v^2 σ_w^22, and C is a coefficient determined by initial conditions.




	Suppose α < 1.
	We have the following asymptotic expansion
	
    l = K_1 l^ + R(l)

	where the remainder term
	
    R(l) ∼
    	-K_2 l^log l	   if α >  1 2
    
    	(C-K_2) l log l			   if α =  1 2 and K_2 ≠C
    C(1-α)1-2α l	   if α <  1 2

	where K_1 = [σ_v^2 σ_w^2α_α (1 - α)]^ 1 1-α, K_2 =  1 2 [σ_v^2 _ασ_w^2α]^ (1-α)^ - 1α and C = σ_a^2.

<ref> verifies the leading coefficient and the exponent of the leading term.

The difference equation governing the evolution of  is

    -  = A( + B)^α + C

where A =σ_v^2 _ασ_w^2α, B = σ_b^2/σ_w^2, and C = σ_a^2.
Then <ref> yields the result.




<ref> combined with <ref> gives the following result.
*


By <cit.>, we know that ψ_α( q, q c) = ψ_α( q) _α(c), where _α(c) = J_α(arccos c) and

    J_α(θ) :=  1 2π_α (sinθ)^2α + 1Γ(α + 1)∫_0^π / 2ηcos^αη(1 - cosθcosη)^1 + α.


Note that _α(c) ∈ (-∞, ∞) for α∈ (-1, ∞) and any c ∈ (0, 1), even though ψ_α is only defined for α > -1/2.


<ref> shows a comparison of _α for different αs along with the identity function.
By <cit.>, _α is an increasing and convex function as long as ψ_α^2 is Gaussian-integrable, which is precisely when α > -1/2.
We can compute _α(1) = ψ_α(q, q)/ψ_α(q) = 1, and _α(0) = ψ_α(q, 0)/ψ_α(q) = ψ_α/2(q)^2/ψ_α(q) = _α/2^2/_α =  12√(π)Γ(α 2 +  1 2)^2Γ(α +  1 2).
We record these observations as a lemma.

	_α(c) is an increasing and convex function for each α > -1/2 on c ∈ [0, 1].
	_α(1) = 1 and _α(0) =  12√(π)Γ(α 2 +  1 2)^2Γ(α +  1 2).



For α = 1, <cit.> computed

    _1(c) =  1 π (√(1 - c^2) + (π - arccos(c))c).

<ref> shows a plot of _1 vs identity.
It has derivative _1(c) = 1 -  1 πarccos c, which shows that _1(c) < 1 with equality iff c = 1, and consequently _1(c) ≥ c with equality iff c = 1.
At the same time, _1(c) ≥ 0 with equality iff c = -1, so _1 is increasing on [-1, 1].
It has an asymptotic expansion _1(1 - ε) = 1 - ε + 2√(2)3π^3/2 + Θ(^5/2) at 1.

The zeroth Bessel function of the second kind is defined by (z) = ∫_1^∞ e^-z x (x^2 - 1)^-1/2 x.
It is one of the fundamental solutions to the homogeneous differential equation x^2 ẏ + x ẏ - x^2 y = 0.
The following lemma shows that J_α can be expressed in terms of .

	For any α > -1,
	J_α(θ) =  1 2π_αsin^2α+1θ∫_0^∞ x (x) e^xcosθ x^α


	<cit.> gave the expression
	
    2π_α J_α(θ) = θ∫_0^∞ u ∫_0^∞ v e^-(u^2+v^2 - 2uv cosθ)/2sin^2θ u^α v^α.

	Note that the integrand is symmetric in u and v.
	
	Thus, if = {(u, v): u, v ≥ 0  v ≥ u}, then
	
    2π_α J_α(θ) = 2θ∫_ u  v e^-(u^2+v^2 - 2uv cosθ)/2sin^2θ u^α v^α.

	
	
	Now make the change of variables from  to {(, ): ≥ 2 √()}:
	
    = uv	      = v  u + u  v
       = u + v      =  u +  v
       = (v - u)  u  v	    u  v    = (^2 - 4 )^-1/2

	so that we have
	
    2 π_α J_α(θ) = 2θ∫_0^∞ e^(1 + cosθ)^2θ^α∫_2√()^∞ e^-^2^2 θ(^2 - 4 )^-1/2.

	The inner integral in  can be expressed in terms of  by a change of variable x = ^2/2√():
	
    2 π_α J_α(θ)    = 2θ∫_0^∞ e^(1 + cosθ)^2θ^α 1 2 e^-^2 θ(^2 θ)
       = θ∫_0^∞(^2θ) e^cosθ^2θ^α
       = sin^2α + 1θ∫_0^∞ x (x) e^x cosθ x^α



Define L_α(θ) = 2 π_α J_α(θ)^2α+1θ = ∫_0^∞ x (x) e^x cosθ x^α.

	If α > 1, then
	
    L_α(θ) = ^2 θ [(2α-1) cosθ L_α-1(θ) + (α-1)^2 L_α - 2(θ)].



	We will prove this claim for θ < 1, and by continuity this also proves the case θ = 1.
	As remarked above, (z) = (z) +  z (z).
	Thus
	
    L_α(θ)    = ∫_0^∞ x ((x) +  x (x)) e^x cosθ x^α
       =e^xcosθ x^α|_0^∞ +  e^xcosθ x^α - 1|_0^∞
       = - ∫ x[cosθ e^x cosθ x^α + α e^x cosθ x^α - 1] 
       = - ∫ x[cosθ e^x cosθ x^α - 1 + (α-1) e^xcosθ x^α - 2]

	Asymptotically, (z) ∼√(π2z) e^-z as z →∞ and (z) ∼ -ln(z) as z ↘ 0, and (z) ∼ -√(π2z) e^-z as z →∞ and (z) ∼ - z as z ↘ 0.
	Thus, as α > 1,
	
    e^xcosθ x^α|_0^∞ = -lim_x →∞√(π/2) e^-x(1 - cosθ)x^α - 1 + lim_x↘ 0 e^xcosθ x^α-1  = 0
     e^xcosθ x^α-1|_0^∞ = -lim_x →∞√(π/2) e^-x(1 - cosθ)x^α - 2 + lim_x↘ 0 e^xcosθ x^α-1ln x  = 0

	So
	
    L_α(θ)    = -cosθ L_α-1(θ) - (α-1)L_α-2(θ) - ∫ x[cosθ e^x cosθ x^α + α e^x cosθ x^α - 1]

	Via another integration by parts, the integral on the right is
	
    =cosθ e^x cosθ x^α|_0^∞ + α e^xcosθ x^α-1|_0^∞
       = - ∫ x[cos^2 θ e^x cosθ x^α + 2 αcosθ e^x cosθ x^α - 1 + α (α - 1) e^x cosθ x^α - 2] 
       = - [cos^2 θ L_α(θ) + 2 αcosθ L_α-1(θ) + α(α-1) L_α-2(θ)]

	where the evaluation terms vanish just like before.
	Altogether, we have
	
    L_α(θ)    = cos^2 θ L_α(θ) + (2α-1)cosθ L_α-1(θ) + (α-1)^2 L_α-2(θ)
       = ^2 θ[(2α-1) cosθ L_α-1(θ) + (α-1)^2 L_α-2(θ)]


As a corollary we get

	Suppose α > 1.
	Then
	
    J_α(θ)    = cosθ J_α-1(θ) + (α-1)^2 (2α-1)^-1(2α-3)^-1sin^2 θ J_α-2(θ)
    _α(c)    =  c _α-1(c) + (α-1)^2 (2α-1)^-1(2α-3)^-1(1-c^2) _α-2(c)



The derivative of J_α(θ) turns out to be quite simple.

	Suppose α > 0. Then
	
    J̇_α(θ)    = -α^2(2 α-1)^-1 J_α-1(θ) sinθ
    _α(c)    = α^2 (2 α-1)^-1_α-1(c)



	We will prove the first formula.
	The second follows from chain rule.
	By <ref>,
	
    J_α(θ)    =  1 2π_αsin^2α+1θ∫ x (x) e^xcosθ x^α
    J̇_α(θ)    =  1 2π_α[(2α+1)sin^2αθcosθ∫ x (x) e^xcosθ x^α
       = -sin^2α+2θ∫ x (x) e^xcosθ x^α+1]
       = (2α+1) θ J_α(θ) - c_α+1c_αθ J_α+1(θ)
       = (2α+1) θ [cosθ J_α(θ) - J_α+1(θ)].

	As α + 1 > 1, by <ref>, this is
	
    =-(2α+1)θ [(α-1)^2 (2α+1)^-1(2α-1)^-1sin^2 θ J_α-1(θ)]
       = -(α-1)^2(2α-1)^-1sinθ J_α-1(θ).



Thus _α(1) = α^2(2α-1)^-1_α-1(1) = α^2(2α-1)^-1 for any α > 0 by <ref>.
For 1/2 < α≤ 1, _α(1) ≥ 1 with equality iff α = 1, and for α = 1/2, _α(1) =∞ > 1 by continuity of _α(c) in α.
Because for α > -1/2, _α is increasing and convex on [0, 1] and _α(0) > 0 by <ref>, _α intersects identity at a unique point away from 1 when α∈ [1/2, 1).
We record this as a theorem.


	For α∈ [1/2, 1), _α(c) = c has two solutions: an unstable solution at 1 ("unstable" meaning _α(1) > 1) and a stable solution in ^* ∈ (0, 1) ("stable" meaning _α(^*) < 1).




This result confirms that pictures presented in <ref>b,c are qualitatively correct, that there are indeed stable fixed points of _α away from 1.

*

If < 1, then

    = σ_w^2 γ + σ_b^2σ_w^2  + σ_b^2   ≥
    _1()    ≥_1()
     = σ_v^2 _α^α_1() + σ_b^2σ_v^2 _α^α + σ_b^2   ≥_1()

but ≥_1() > as noted above.
Thus by monotone convergence  converges, and ^* = 1 is the only possible fixed point.


By <ref>, = (1 + Θ(^-1)) = 1 -  + Θ() = 1 - u where u := 1 - Θ(^-1).
Using the asymptotic expansion _1(1 - ) = 1 - + U ^3/2 + Θ(^5/2), we have

    (1 - )   = σ_v^2  2 _1(1 - u) + σ_a^2 + (1 - ) 
    
    -   = σ_v^2  2 (_1(1 - u) - 1) - 
       = σ_v^2  2 [- u  + U u^3/2^3/2 + Θ(u^5/2^5/2)] - 
       =  1  [ + σ_v^2  2 (u - U u^3/2^1/2 + Θ(u^5/2^3/2))]
       =   1  [ - σ_a^2 + σ_v^2  2 (Θ() - U u^3/2^1/2 + Θ(u^5/2^3/2))]
       =   [1 + -σ_a^2 + σ_v^2  2 (Θ() - U u^3/2^1/2 + Θ(u^5/2^3/2))]
       =   [1 + -σ_a^2 +  1 2 σ_v^2 (Θ() - U u^3/2^1/2 + Θ(u^5/2^3/2))]

Let the content of the bracket on the RHS be ℵ.
We have = (1+ o(1))B/σ_w^2.
If = O(^-1), then ℵ = 1 - O(), but because  is exponentially decreasing, this means = Θ(1) and does not converge to 0 — this is a contradiction.
Therefore, = ω(^-1), and

    =  [1 -  1 2 B^-1σ_v^2 σ_w^2 U^1/2(1 + o(1))]
     -    = -  1 2 B^-1σ_v^2 σ_w^2 U^3/2(1 + o(1))

Using <ref> to upper and lower bound our dynamics, we get that l ∼ [ 1 4 σ_v^2 σ_w^2  B U l]^-2.




Let ϕ be any nonlinearity.
Suppose ϕ( r, rd) = ϕ( r)(d) for some twice differentiable function (d) independent of , where (1) = 1 naturally.
Suppose further that 

    
  * (d) = d has a solution d = ^* > 0 where (^*) = δ < 1;
    
  * (d) > d for all d < ^* and (d) < d for all 1 > d > ^*; and
    
  *  is nondecreasing. 

Let l :=  l - ^* and suppose 0 < 1.
If γ l →∞ and ϕ(  l) →∞, then l → 0 and satisfies

    = (1 - σ_a^2 + (1 - δ + O())σ_v^2 ϕ( )) + ϕ( )Θ(γ).




	First we note that because ^* is the only stable fixed point of the dynamics x ↦(x), with the basin of attraction [0, 1), we can show l →^* as in the proof of <ref> (using <ref>).
	
Write V l := ϕ(  l).
We first show that l →^*.
When l is large,


    = σ_w^2 γ + σ_b^2σ_w^2  + σ_b^2    =  ( 1 + O(γ))
     = σ_v^2V() + σ_a^2σ_v^2 V + σ_a^2   = ()(1 + O( V ()^-1)).

If γ l is bounded for all l, then → 0 because l →∞.
Since () > 0 for ∈ [0, 1] and V l→∞, we have that in the limit l→∞, lim_l →∞ = 0 = (lim_l →∞) = (0) (by the continuity of ), which is impossible by our assumptions.
Thus γ l →∞, and we have lim_l →∞ = (lim_l →∞).
By our assumptions, ^* is the only stable fixed point of  with basin of attraction [0, 1), so this shows that →^* as desired.

Now we derive the equation in question.
Note that = (1 + Θ(γ)) because ^* < 1.
We use the Taylor expansion (^* + ) = ^* + δ + O(^2).



    (^* + )   =  V (^* + )(1 + Θ(γ)) +  + (^* + ) 
       =  V(^* + δ( + Θ(γ)) + O(^2)) +  + (^* + ) 
       =  V(δ( + Θ(γ)) + O(^2)) + 
       =  (1 -  + (1 - δ + O()) V) + Θ(V γ)




*

We apply <ref>.
We first check the conditions of the lemma, with = _α.
The following conditions were already verified.

    
  * _α has a fixed point ^* less than but very close to 1, where its slope is υ := _α(^*) < 1. (<ref>)
    
    
  * _α(d) > d for all d < ^* and _α(d) < d for all d > ^*. (By the convexity shown in <ref>)
    
  * _α is nondecreasing (<ref>).
    Furthermore, from its integral formula (<ref>), we see easily that _α is smooth at ^* < 1.

We also proved the following

    
  * l ∼ [σ_v^2 σ_w^2α_α (1 - α)]^ 1 1-αl^ 1 1-α (<ref>) and γ l is asymptotically a constant fraction of l (<ref>), so both go to ∞. 
    
  * ψ_α( ) = _α^α = _α (σ_w^2  + σ_b^2)^α = Θ(l^α/(1-α)), so goes to ∞. (<ref>)

Thus, for υ = (^*),

    σ_a^2 + (1 - υ + O())σ_v^2 ϕ( )   ∼(1 - υ)σ_v^2σ_w^2α_α^1-α
       =  l (1-υ)/(1 - α).


Now, ϕ( )γ = Θ(l^-11-α - 1).
By using the dynamics of <ref> to upper and lower bound our dynamics, we have l = Ω(l^-μ - ), O(l^-μ + ) for any > 0, where μ = min((1-υ)/(1-α), 1/(1-α)) = (1-υ)/(1-α).





  §.§.§ Backward Dynamics




	Suppose random variable X ∼(0, σ^2), and Y = ψ_-β(X) for some β > 0, where ψ_α is α-ReLU.
	Then for ξ > 0, Y has density
	
    [Y ∈ [ξ, ξ+ξ]] =  1 β√(2πσ^2)ξ^- 1 β - 1 e^-ξ^-2/β/2σ^2.

	At ξ = 0, Y has density given by a Dirac delta of mass 1 2.
	
	Furthermore, Y has finite second moment iff β <  1 2.


	We have
	
    [Y ∈ [ξ, ∞)]    = [X ∈ [0, ξ^-1/β]]
       =  1 √(2πσ^2)∫_0^ξ^-1/β e^-x^2/2σ^2 x.

	Differentiating the RHS against ξ using Leibniz's rule, we get
	
    d[Y ∈ [ξ, ∞)]/dξ   =  1 √(2πσ^2)  e^-ξ^-2/β/2σ^2 d dξξ^-1/β
       = -1β√(2πσ^2)ξ^- 1 β - 1 e^-ξ^-2/β/2σ^2.

	Negating both sides gives the density f_Y of Y for ξ > 0.
	For ξ = 0, observe that lim_ξ→ 0 f_Y(ξ) = 0 because, while ξ^- 1 β -1 blows up polynomially, e^-ξ^-2/β/2σ^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) =  1 2 δ_0, where δ_0 is the Dirac delta.
	
	For the second assertion, observe that
	
    f_Y(ξ) ∼ 1 β√(2πσ^2)ξ^- 1 β - 1   as ξ→∞.

	Thus, ξ^2 f_Y(ξ) is integrable iff 2 - 1 β - 1 < -1 β <  1 2.

*

	Note that ψ̇_α∝ψ_α - 1, so it suffices to show that (ψ_α-1(ζ)^2) = (ψ_2α - 2(ζ)) is infinite for ζ∼(0, σ^2).
	By <ref> with β = 2 - 2α, ψ_2α - 2(ζ) has finite variance iff β <  1 2 α >  3 4.



*

If α = 1, then

    ℸ = ℸ (1 +  1 2σ_v^2 σ_w^2).

So ℸl-m / ℸ l = Θ(1) B^m for B = 1 +  1 2σ_v^2 σ_w^2.

If 1 2 < α < 1, then ℸ/ℸ - 1 is

    =σ_v^2 σ_w^2 ϕ̇( )
       = σ_v^2 σ_w^2 α^2 _α - 1^α - 1
       = σ_v^2 σ_w^2 α^2 _α - 1 (σ_w^2 )^α - 1 + Θ(^α - 2)
       = σ_v^2 σ_w^2αα^2 _α - 1(K_1 l^ - K_2 l^log l + o(l^log l))^α - 1 + Θ(l^α - 21 - α)    by <ref>
       = σ_v^2 σ_w^2αα^2 _α - 1 [K_1^α - 1 l + Θ(l^-2log l)] + O(l^-3)
       = σ_v^2 σ_w^2αα^2 _α - 1 K_1^α - 1 l + Θ(l^-2log l)
       = R  l + Θ(l^-2log l)

where R = σ_v^2 σ_w^2αα^2 _α - 1 K_1^α - 1 = α^2(1-α)(2 α - 1) and K_1 = [σ_v^2 σ_w^2α_α (1 - α)]^ 1 1-α.
So

    ℸ   = ℸexp(R  l + Θ(l^-2log l))
    ℸl -m   = Θ(1)ℸl l l-m^R

as desired.


*

	The proof is similar to that of <ref>.

