Title: A MEAN FIELD THEORY OF BATCH NORMALIZATION

Abstract: We develop a mean field theory for batch normalization in fully-connected feedforward neural networks. In so doing, we provide a precise characterization of signal propagation and gradient backpropagation in wide batch-normalized networks at initialization. Our theory shows that gradient signals grow exponentially in depth and that these exploding gradients cannot be eliminated by tuning the initial weight variances or by adjusting the nonlinear activation function. Indeed, batch normalization itself is the cause of gradient explosion. As a result, vanilla batch-normalized networks without skip connections are not trainable at large depths for common initialization schemes, a prediction that we verify with a variety of empirical simulations. While gradient explosion cannot be eliminated, it can be reduced by tuning the network close to the linear regime, which improves the trainability of deep batch-normalized networks without residual connections. Finally, we investigate the learning dynamics of batch-normalized networks and observe that after a single step of optimization the networks achieve a relatively stable equilibrium in which gradients have dramatically smaller dynamic range. Our theory leverages Laplace, Fourier, and Gegenbauer transforms and we derive new identities that may be of independent interest.

1. INTRODUCTION:

Deep neural networks have been enormously successful across a broad range of disciplines. These successes are often driven by architectural innovations. For example, the combination of convolutions (LeCun et al., 1990) , residual connections (He et al., 2015) , and batch normalization (Ioffe & Szegedy, 2015) has allowed for the training of very deep networks and these components have become essential parts of models in vision (Zoph et al.) , language (Chen & Wu, 2017) , and reinforcement learning (Silver et al., 2017) . However, a fundamental problem that has accompanied this rapid progress is a lack of theoretical clarity. An important consequence of this gap between theory and experiment is that two important issues become conflated. In particular, it is generally unclear whether novel neural network components improve generalization or whether they merely increase the fraction of hyperparameter configurations where good generalization can be achieved. Resolving this confusion has the promise of allowing researchers to more effectively and deliberately design neural networks.

Recently, progress has been made (Poole et al., 2016; Schoenholz et al., 2016; Daniely et al., 2016; Pennington et al., 2017; Hanin & Rolnick, 2018; Yang, 2019) in this direction by considering neural networks at initialization, before any training has occurred. In this case, the parameters of the network are random variables which induces a distribution of the activations of the network as well as the gradients. Studying these distributions is equivalent to understanding the prior over functions that these random neural networks compute. Picking hyperparameters that correspond to well-conditioned priors ensures that the neural network will be trainable and this fact has been extensively verified experimentally. However, to fulfill its promise of making neural network design less of a black box, these techniques must be applied to neural network architectures that are used in practice. Over the past year, this gap has closed significantly and theory for networks with skip connections (Yang & Schoenholz, 2017; 2018) , convolutional networks (Xiao et al., 2018) , and gated recurrent networks (Chen et al., 2018; Gilboa et al., 2019) have been developed. More recently, Yang (2019) devised a formalism that extends this approach to include an even wider range of architectures.

Before state-of-the-art models can be analyzed in this framework, a slowly-decreasing number of architectural innovations must be studied. One particularly important component that has thus-far remained elusive is batch normalization.

Our Contributions. In this paper, we develop a theory of fully-connected networks with batch normalization whose weights and biases are randomly distributed. A significant complication in the case of batch normalization (compared to e.g. layer normalization or weight normalization) is that the statistics of the network depend non-locally on the entire batch. Thus, our first main result is to recast the theory for random fully-connected networks so that it can be applied to batches of data. We then extend the theory to include batch normalization explicitly and validate this theory against Monte-Carlo simulations. We show that as in previous cases we can leverage our theory to predict valid hyperparameter configurations.

In the process of our investigation, we identify a number of previously unknown properties of batch normalization that make training unstable. In particular, for most nonlinearities used in practice, batchnorm in a deep, randomly initialized network induces high degree of symmetry in the embeddings of the samples in the batch (Thm 3.4) . Whenever this symmetry takes hold, we show that for any choice of nonlinearity, gradients of fully-connected networks with batch normalization explode exponentially in the depth of the network (Thm 3.9). This imposes strong limits on the maximum trainable depth of batch normalized networks. This limit can be lifted partially but not completely by pushing activation functions to be more linear at initialization. It might seem that such gradient explosion ought to lead to learning dynamics that are unfavorable. However, we show that networks with batch normalization causes the scale of the gradients to naturally equilibrate after a single step of gradient descent (provided the initial gradients are not so large as to cause numerical instabilities). For shallower networks, this equilibrating effect is sufficient to allow adequate training.

Finally, we note that there is a related vein of research that has emerged that leverages the prior over functions induced by random networks to perform exact Bayesian inference (Lee et al., 2017; de G. Matthews et al., 2018; Novak et al., 2019; Garriga-Alonso et al., 2019; Yang, 2019) . One of the natural consequences of this work is that the prior for networks with batch normalization can be computed exactly in the wide network limit. As such, it is now possible to perform exact Bayesian inference in the case of wide neural networks with batch normalization.

2. RELATED WORK:

Batch normalization has rapidly become an essential part of the deep learning toolkit. Since then, a number of similar modifications have been proposed including layer normalization (Ba et al., 2016) and weight normalization (Salimans & Kingma, 2016) . Comparisons of performance between these different schemes have been challenging and inconclusive (Gitman & Ginsburg, 2017) . The original introduction of batchnorm in Ioffe & Szegedy (2015) proposed that batchnorm prevents "internal covariate shift" as an explanation for its effectiveness. Since then, several papers have approached batchnorm from a theoretical angle, especially following Ali Rahimi's catalyzing call to action at NIPS 2017. Balduzzi et al. (2017) found that batchnorm in resnets allow deep gradient signal propagation in contrast to the case without batchnorm. Santurkar et al. (2018) found that batchnorm does not help covariate shift but helps by smoothing loss landscape. Bjorck et al. (2018) reached the opposite conclusion as our paper for residual networks with batchnorm, that batchnorm works in this setting because it induces beneficial gradient dynamics and thus allows a much bigger learning rate. Luo et al. (2018) explores similar ideas that batchnorm allows large learning rates and likewise uses random matrix theory to support their claims. Kohler et al. (2018) identified situations in which batchnorm can provably induce acceleration in training. Of the above that mathematically analyze batchnorm, all but Santurkar et al. (2018) make simplifying assumptions on the form of batchnorm and typically do not have gradients flowing through the batch variance. Even Santurkar et al. (2018) only analyzes a vanilla network which gets added a single batchnorm at a single moment in training. Our analysis here on the other hand works for networks with arbitrarily many batchnorm layers with very general activation functions, and indeed, this deep stacking of batchnorm is precisely what leads to gradient explosion. It is an initialization time analysis for networks of infinite width, but we show experimentally that the resulting insights predict both training and test time behavior.

We remark that Philipp & Carbonell (2018) has also empirically noted gradient explosion happens in deep batchnorm networks with various nonlinearities. In this work, in contrast, we develop a precise theoretical characterization of gradient statistics and, as a result, we are able to make significantly stronger conclusions.

3. THEORY:

We begin with a brief recapitulation of mean field theory in the fully-connected setting. In addition to recounting earlier results, we rephrase the formalism developed previously to compute statistics of neural networks over a batch of data. Later, we will extend the theory to include batch normalization. We consider a fully-connected network of depth L whose layers have width N l , activation function1 φ, weights W l ∈ R N l-1 ×N l , and biases b l ∈ R N l . Given a batch of B inputs2 {x i : x i ∈ R N0 } i=1,••• ,B , the pre-activations of the network are defined by the recurrence relation,

EQUATION

At initialization, we choose the weights and biases to be i.i.d. as W l αβ ∼ N (0, σ 2 w /N l-1 ) and b l α ∼ N (0, σ 2 b ). In the following, we use α, β, . . . for the neuron index, and i, j, . . . for the batch index. We will be concerned with understanding the statistics of the pre-activations and the gradients induced by the randomness in the weights and biases. For ease of exposition we will typically take the network to have constant width N l = N .

In the mean field approximation, we iteratively replace the pre-activations in Eq. ( 2) by Gaussian random variables with matching first and second moments. In the infinite width limit this approximation becomes exact (Lee et al., 2017; de G. Matthews et al., 2018; Yang, 2019) . Since the weights are i.i.d. with zero mean it follows that the mean of each pre-activation is zero and the covariance between distinct neurons are zero. The pre-activation statistics are therefore given by (h l α11 , • • • , h l α B B )

N l-1 →∞ ------→ N (0, Σ l δ α1•••α B ) where Σ l are B × B covariance matrices and δ α1,•••α B is the Kronecker-δ that is one if α 1 = α 2 = • • • = α B and zero otherwise. Definition 3.1. Let V be the operator on functions φ such that V φ (Σ) = E[φ(h)φ(h) T : h ∼ N (0, Σ)] computes the matrix of uncentered second moments of φ(h) for h ∼ N (0, Σ).

Using the above notation, we can express the covariance matrices by the recurrence relation,

EQUATION

At first Eq. ( 2) may seem challenging since the expectation involves a Gaussian integral in R B . However, each term in the expectation of V φ involves at most a pair of pre-activations and so the expectation may be reduced to the evaluation of O(B 2 ) two-dimensional integrals. These integrals can either be performed analytically (Cho & Saul, 2009; Williams, 1997) or efficiently approximated numerically (Lee et al., 2017) , and so Eq. ( 2) defines a computationally efficient method for computing the statistics of neural networks after random initialization. This theme of dimensionality reduction will play a prominent role in the forthcoming discussion on batch normalization.

Eq. ( 2) defines a dynamical system over the space of covariance matrices. Studying the statistics of random feed-forward networks therefore amounts to investigating this dynamical system and is an enormous simplification compared with studying the pre-activations of the network directly. As is common in the dynamical systems literature, a significant amount of insight can be gained by investigating the behavior of Eq. (2) in the vicinity of its fixed points. For most common activation functions, Eq. (2) has a fixed point at some Σ * . Moreover, when the inputs are non-degenerate, this fixed point generally has a simple structure with Σ * = q * [(1c * )I + c * 11 T ] owing to permutation symmetry among elements of the batch. We refer to fixed points with such symmetry as BSB1 (Batch Symmetry Breaking 1 or 1 Block Symmetry Breaking) fixed points. As we will discuss later, in the context of batch normalization other fixed points with fewer symmetries ("BSBk fixed points") may become preferred. In the fully-connected setting fixed points may efficiently be computed by solving the fixed point equation induced by Eq. ( 2) in the special case B = 2. The structure of this fixed point implies that in asymptotically deep feed-forward neural networks all inputs yield pre-activations of identical norm with identical angle between them. Neural networks that are deep enough so that their pre-activation statistics lie in this regime have been shown to be untrainable (Schoenholz et al., 2016) .

Notation As we often talk about matrices and also linear operators over matrices, we write T {Σ} for an operator T applied to a matrix Σ, and matrix multiplication is still written as juxtaposition.

Composition of matrix operators are denoted with T 1 • T 2 .

Local Convergence to BSB1 Fixed Point. To understand the behavior of Eq. ( 2) near its BSB1 fixed point we can consider the Taylor series in the deviation from the fixed point, ∆Σ l = Σ l -Σ * .

To lowest order we generically find,

EQUATION

where J = dV φ dΣ Σ=Σ * is the B 2 × B 2 Jacobian of V φ . In most prior work where φ was a pointwise non-linearity, Eq. ( 3) reduces to the case B = 2 which naturally gave rise to linearized dynamics in q l = E[(h l i ) 2 ] and c l = E[h l i z l j ]/q l . However, in the case of batch normalization we will see that one must consider the evolution of Eq. ( 3) as a whole. This is qualitatively reminiscent of the case of convolutional networks studied in Xiao et al. (2018) where the evolution of the entire pixel × pixel covariance matrix had to be evaluated.

The dynamics induced by Eq. (3) will be controlled by the eigenvalues of J: Suppose J has eigenvalues λ i -ordered such that λ 1 ≥ λ 2 ≥ • • • ≥ λ B 2 -with associated eigen"vectors" e i (note that the e i will themselves be B × B matrices). It follows that if ∆Σ 0 = i c i e i for some choice of constants c i then ∆Σ l = i c i λ l i e i . Thus, if λ i < 1 for all i, ∆Σ l will approach zero exponentially and the fixed-point will be stable. The number of layers over which Σ will approach Σ * will be given by -1/ log(λ 1 ). By contrast if λ i > 1 for any i then the fixed point will be unstable. In this case, there is typically a different, stable, fixed point that must be identified. It follows that if the eigenvalues of J can be computed then the dynamics will follow immediately.

While J may appear to be a complicated object at first, a moment's thought shows that each diagonal element J{∆Σ} ii is only affected by the corresponding diagonal element ∆Σ ii , and each off-diagonal element J{∆Σ} ij is only affected by ∆Σ ii , ∆Σ jj , and ∆Σ ij . Such a Diagonal-Offdiagonal Semidirect (DOS) operator has a simple eigendecomposition with two eigenspaces corresponding to changes in the off-diagonal and changes in the diagonal (Thm E.72). The associated eigenvalues are precisely those calculated by Schoenholz et al. (2016) in a simplified analysis. DOS operators are a particularly simple form of a more general operator possessing an abundance of symmetries called ultrasymmetric operators, which will play a prominent role in the analysis of batchnorm below.

Gradient Dynamics. Similar arguments allow us to develop a theory for the statistics of gradients. The backpropogation algorithm gives an efficient method of propagating gradients from the end of the network to the earlier layers as,

EQUATION

)

Here L is the loss function and δ l i = ∇ h l i L are N l -dimensional vectors that describe the error signal from neurons in the l'th layer due to the i'th element of the batch. The preceding discussion gave a precise characterization of the statistics of the h l i that we can leverage to understand the statistics of δ l i . Assuming gradient independence, that is, that an iid set of weights are used during backpropagation (see Appendix B for more discussions), it is easy to see that E[δ l αi ] = 0 and E[δ l αi δ l βj ] = Π l ij if α = β and 0 otherwise, where Π l is a covariance matrix and we may once again drop the neuron index. We can construct a recurrence relation to compute Π l ,

EQUATION

Typically, we will be interested in understanding the dynamics of Π l when Σ l has converged exponentially towards its fixed point. Thus, we study the approximation,

EQUATION

Since these dynamics are linear (and in fact componentwise), explosion and vanishing of gradients will be controlled by V φ (Σ * ).

3.1. BATCH NORMALIZATION:

We now extend the mean field formalism to include batch normalization. Here, the definition for the neural network is modified to be the coupled equations,

h l i = W l φ( hl-1 i ) + b l hl αi = γ α h l αi -µ α σ α + β α (7)

where γ α and β α are parameters, and

µ α = 1 N l i h αi and σ 2 α = 1 N l i (h αi -µ α ) 2

+ are the per-neuron batch statistics. In practice ≈ 10 -5 or so to prevent division by zero, but in this paper, unless stated otherwise (in the last few sections), is assumed to be 0. Unlike in the case of vanilla fully-connected networks, here the pre-activations are invariant to σ 2 w and σ 2 b . Without a loss of generality, we therefore set σ 2 w = 1 and σ 2 b = 0 for the remainder of the text. In principal, batch normalization additionally yields a pair of hyperparameters γ and β which are set to be constants. However, these may be incorporated into the nonlinearity and so without a loss of generality we set γ = 1 and β = 0. In order to avoid degenerate results, we assume B ≥ 4 unless stated otherwise; we shall discuss the small B regime in Appendix J.

If one treats batchnorm as a "batchwise nonlinearity", then the arguments from the previous section can proceed identically and we conclude that as the width of the network grows, the pre-activations will be jointly Gaussian with identically distributed neurons. Thus, we arrive at an analogous expression to Eq. (2),

Σ l = V B φ (Σ l-1 )

where

EQUATION

)

Here we have introduced the projection operator G = I -1 B 11 T which is defined such that Gx = x -µ1 with µ = i x i /B. Unlike φ, B φ does not act component-wise on h. It is therefore not obvious whether V B φ can be evaluated without performing a B-dimensional Gaussian integral.

Theoretical tools. In this paper, we present several ways to analyze high dimensional integrals like the above: 1. the Laplace method 2. the Fourier method 3. spherical integration 4. and the Gegenbauer method. The former two use the Laplace and Fourier transforms to simplify expressions like the above, where the Laplace method requires that φ be positive homogeneous. Often the Laplace method will give clean, closed form answers for such φ. Because batchnorm can be thought of as a linear projection (G) followed by projection to the sphere of radius √ B, spherical integration techniques are often very useful and in fact is typically the most straightforward way of numerically evaluating quantities. Lastly, the Gegenbauer method expresses objects in terms of the Gegenbauer coefficients of φ. Briefly, Gegenbauer polynomials {C

(α) l (x)} ∞ l=0 are orthogonal polynomials with respect to the measure (1 -x 2 ) α-1 2 on [-1, 1].

They are intimately related to spherical harmonics (a natural basis for functions on a sphere), which explains their appearance in this context. The Gegenbauer method is the most illuminating amongst them all, and is what allows us to conclude that gradient explosion happens regardless of nonlinearity under general conditions. See Appendix D for a more in-depth discussion of these techniques. In what follows, we will mostly present results by the Laplace method for φ = relu and by the Gegenbauer method for general φ, but give pointers to the appendix for others.

Back to the topic of Eq. ( 8), we present a pair of results that expresses Eq. ( 8) in a more manageable form. From previous work (Poole et al., 2016) , V φ can be expressed in terms of a two-dimensional Gaussian integrals independent of B. When φ is degree-α positive homogeneous (e.g. rectified linear activations) we can relate V φ and V B φ by the Laplace transform. (see Thm E.5).

Theorem 3.2. Suppose φ : R → R is degree-α positive homogeneous. For any positive semi-definite matrix Σ define the projection Σ G = GΣG. Then

EQUATION

)

whenever the integral exists.

Using this parameterization, when V φ has a closed form solution (like ReLU), V B φ involves only a single integral. A similar expression can be derived for general φ by applying Fourier transform; see Appendix E.2. Next we express Eq. ( 8) as a spherical integral (see Proposition E.34)

Theorem 3.3. Let e ∈ R B×B-1 have as columns a set of orthonormal basis vectors of {x ∈ R B :

EQUATION

)

Together these theorems provide analytic recurrence relations for random neural networks with batch normalization over a wide range of activation functions. By analogy to the fully-connected case we would like to study the dynamical system over covariance matrices induced by these equations.

We begin by investigating the fixed point structure of Eq. ( 8). As in the case of feed-forward networks, permutation symmetry implies that there exist BSB1 fixed points

Σ * = q * [(1-c * )I+c * 11 T ].

We will see that this fixed point is in fact unique, and a clean expression of q * and c * can be obtained in terms of Gegenbauer basis (see Thm F.13).

Theorem 3.4 (Gegenbauer expansions of BSB1 fixed point).

If φ( √ B -1x) has Gegenbauer ex- pansion ∞ l=0 a l 1 c B-1,l C ( B-3 2 ) l (x) where c B-1,l = B-3 B-3+2l , then q * = ∞ l=0 a 2 l 1 c B-1,l C ( B-3 2 ) l (1), q * c * = ∞ l=0 a 2 l 1 c B-1,l C ( B-3 2 ) l -1 B -1 .

Thus the entries of the BSB1 fixed point are diagonal quadratic forms of the Gegenbauer coefficients of φ( √ B -1x). Even more concise closed forms are available when the activation functions are degree α positive homogeneous (see Thm F.8). In particular, for ReLU we arrive at the following Theorem 3.5 (BSB1 fixed point for ReLU). When φ = relu, then

q * = 1 2 , c * = J 1 -1 B -1 = 1 π - 1 2(B -1) + O 1 (B -1) 2 (11)

where

J 1 (c) = 1 π ( √ 1 -c 2 + (π -arccos(c))c

) is the arccosine kernel (Cho & Saul, 2009) .

In the appendix, Thm F.5 also describes a trigonometric integral formula for general nonlinearities.

In the presence of batch normalization, when the activation function grows quickly, a winner-takeall phenomenon can occur where a subset of samples in the batch have much bigger activations than others. This causes the covariance matrix to form blocks of differing magnitude, breaking the BSB1 symmetry. One notices this, for example, as the degree α of α-ReLU (i.e. αth power of ReLU) increases past a point α transition (B) depending on the batch size B (see Fig. 1 ). However, we observe, through simulations, that by far most of the nonlinearities used in practice, like ReLU, leaky ReLU, tanh, sigmoid, etc, all lead to BSB1 fixed points, and we can prove this rigorously for φ = id (see Corollary F.3). We discuss symmetry-broken fixed points (BSB2 fixed points) in Appendices J and K, but in the main text, from here on, Assumption 1. Unless stated otherwise, we assume that any nonlinearity φ mentioned induces Σ l to converge to a BSB1 fixed point under the dynamics of Eq. (8).

3.1.1. LINEARIZED DYNAMICS:

With the fixed point structure for batch normalized networks having been described, we now investigate the linearized dynamics of Eq. ( 8) in the vicinity of these fixed points.

To determine the eigenvalues of

dV B φ

dΣ Σ=Σ * it is helpful to consider the action of batch normalization in more detail. In particular, we notice that B φ can be decomposed into the composition of three separate operations, B φ = φ • n • G. As discussed above, Gh subtracts the mean from h and we introduce the new function n(h) = √ Bh/||h|| which normalizes a centered h by its standard deviation. Applying the chain rule, we can rewrite the Jacobian as,

EQUATION

where • denotes composition and G ⊗2 is the natural extension of G to act on matrices as

G ⊗2 {Σ} = GΣG = Σ G . It ends up being advantageous to study G ⊗2 • dV [φ•n] dΣ Σ=Σ * • G ⊗2 =: G ⊗2 • Ĵ • G ⊗2

and to note that the nonzero eigenvalues of this object are identical to the nonzero eigenvalues of the Jacobian (see Lemma F.17).

At face value, this is a complicated object since it simultaneously has large dimension and possesses an intricate block structure. However, the permutation symmetry of the BSB1 Σ * induces strong symmetries in Ĵ that significantly simplify the analysis (see Appendix F.3). In particular while Ĵijkl is a four-index object, we have Ĵijkl = Ĵπ(i)π(j)π(k)π(l) for all permutations π on B and Ĵijkl = Ĵjilk . We call linear operators possessing such symmetries ultrasymmetric (Defn E.53) and show that all ultrasymmetric operators conjugated by G ⊗2 admit an eigendecomposition that contains three distinct eigenspaces with associated eigenvalues (see Thm E.62).

Theorem 3.6. Let T be an ultrasymmetric matrix operator. Then on the space of symmetric matrices, G ⊗2 • T • G ⊗2 has the following orthogonal (under trace inner product) eigendecomposition, 1. an eigenspace {Σ : Σ G = 0} with eigenvalue 0. M . The specific forms of the eigenvalues can be obtained as linear functions of the entries of T ; see Thm E.62 for details. Note that, as the eigenspaces are orthogonal, this implies that G ⊗2 • T • G ⊗2 is self-adjoint (even when T is not).

In our context with T = Ĵ, the eigenspaces can be roughly interpreted as follows: The deviation ∆Σ = Σ -Σ * from the fixed point decomposes as a linear combination of components in each of the eigenspaces. The RG-component captures the average norm of elements of the batch (the trace of ∆Σ), the L-component captures the fluctuation of such norms, and the M-component captures the covariances between elements of the batch.

Because of the explicit normalization of batchnorm, one sees immediately that the RG-component goes to 0 after 1 step. For positive homogeneous φ, we can use the Laplace method to obtain closed form expressions for the other eigenvalues (see Thm F.33). The below theorem shows that, as the batch size becomes larger, a deep ReLU-batchnorm network takes more layers to converge to a BSB1 fixed point.

EQUATION

where

µ * = q * (1 -c * ) = 1 2 1 -J 1 -1 B-1

, and denotes increasing limit as B → ∞.

Published as a conference paper at ICLR 2019 More generally, we can evaluate them for general nonlinearity using spherical integration (Appendix F.3.1) and, more enlightening, using the Gegenbauer method, is the following

EQUATION

where the coefficients w B-1,l , u B-1,l , wB-1,l , ũB-1,l , v B-1,l are given in Thms F.22 and F.24.

A BSB1 fixed point is not locally attracting if λ ↑ L > 1 or λ ↑ M > 1. Thus Thm 3.8 yields insight on the stability of the BSB1 fixed point, which we can interpret heuristically as follows. The specific forms of the coefficients w

B-1,l , u B-1,l , wB-1,l , ũB-1,l , v B-1,l show that λ ↑

M is typically much smaller than λ ↑ L (but there are exceptions like φ = sin), and

w B-1,1 < v B-1,1 but w B-1,l ≥ v B-1

,l for all l ≥ 2. Thus one expects that, the larger a 1 is, i.e. the "more linear" and less explosive φ is, the smaller λ ↑ L is and the more likely that Eq. ( 8) converges to a BSB1 fixed point. This is consistent with the "winner-take-all" intuition for the emergence of BSB2 fixed point explained above. See Appendix F.3.2 for more discussion.

3.1.2. GRADIENT BACKPROPAGATION:

With a mean field theory of the pre-activations of feed-forward networks with batch normalization having been developed, we turn our attention to the backpropagation of gradients. In contrast to the case of networks without batch normalization, we will see that exploding gradients at initialization are a severe problem here. To this end, one of the main results from this section will be to show that fully-connected networks with batch normalization feature exploding gradients for any choice of nonlinearity such that Σ l converges to a BSB1 fixed point. Below, by rate of gradient explosion we mean the β such that the gradient norm squared grows as β L+o(L) with depth L. As before, all computations below assumes gradient independence (see Appendix B for a discussion).

As a starting point we seek an analog of Eq. ( 6) in the case of batch normalization. However, because the activation functions no longer act point-wise on the pre-activations, the backpropagation equation becomes,

EQUATION

)

where h l α = (h l α1 , . . . , h l αB ) and we observe the additional sum over the batch. Computing the resulting covariance matrix Π l , we arrive at the recurrence relation,

Π l = E dB φ (h) dh T Π l+1 dB φ (h) dh : h ∼ N (0, Σ l ) =: V B φ (Σ l ) † {Π l+1 } (17)

where we have defined the linear operator

V F (Σ) † {•} such that V F (Σ) † {Π} = E[F T h ΠF h : h ∼ N (0, Σ)]

for any vector-indexed linear operator F h . As in the case of vanilla feed-forward networks, here we will be concerned with the behavior of gradients when Σ l is close to its fixed point. We therefore study the asymptotic approximation to Eq. ( 17) given by Π l = V B φ (Σ * ) † {Π l+1 }. In this case the dynamics of Π are linear and are therefore naturally determined by the eigenvalues of

V B φ (Σ * ) † .

As in the forward case, batch normalization is the composition of three operations

B φ = φ • n • G.

Applying the chain rule, Eq. ( 17) can be rewritten as,

EQUATION

with F (Σ) appropriately defined. Note that since G ⊗2 is an idempotent operator,

(V B φ (Σ) † ) n = (G ⊗2 • F (Σ)) n = (G ⊗2 • F (Σ) • G ⊗2 ) n-1 • G ⊗2 • F (Σ), so that it suffices to study the eigende- composition of G ⊗2 • F (Σ * ) • G ⊗2 . Due to the symmetry of Σ * , F (Σ * ) is ultrasymmetric, so that G ⊗2 • F (Σ * ) • G ⊗2

has eigenspaces RG, L, M and we can compute its eigenvalues via Thm 3.6. More illuminating, however, is the Gegenbauer expansion (see Thm G.5). It requires a new identity Thm E.47 involving Gegenbauer polynomials integrated over a sphere, which may be of independent interest. Theorem 3.9 (Batchnorm causes gradient explosion). Suppose φ(

√ B -1x), φ ( √ B -1x) ∈ L 2 ((1 -x 2 ) B-3 2 ). If φ( √ B -1x) has Gegenbauer expansion ∞ l=0 a l 1 c B-1,l C ( B-3 2 ) l (x), then gra- dients explode at the rate of λ ↓ G = ∞ l=0 l+B-3 B-3 • l a 2 l r l ∞ l=0 a 2 l r l where r l = c -1 B-1,l C ( B-3 2 ) l (1) -C ( B-3 2 ) l -1 B-1 > 0.

Consequently, for any non-constant φ (i.e. there is a j > 0 such that a j = 0), λ ↓ G > 1; φ minimizes λ ↓ G iff it is linear (i.e. a i = 0, ∀i ≥ 2), in which case gradients explode at the rate of B-2 B-3 . This contrasts starkly with the case of non-normalized fully-connected networks, which can use the weight and bias variances to control its mean field network dynamics (Poole et al., 2016; Schoenholz et al., 2016) . As a corollary, we disprove the conjecture of the original batchnorm paper (Ioffe & Szegedy, 2015) that "Batch Normalization may lead the layer Jacobians to have singular values close to 1" in the initialization setting, and in fact prove the exact opposite, that batchnorm forces the layer Jacobian singular values away from 1.

Appendix G.1 discusses the numerical evaluation of all eigenvalues, and as usual, the Laplace method yields closed forms for positive homogeneous φ (Thm G.12). We highlight the result for ReLU. Theorem 3.10. In a ReLU-batchnorm network, the gradient norm explodes exponentially at the rate of

EQUATION

)

which decreases to π π-1 ≈ 1.467 as B → ∞. In contrast, for a linear batchnorm network, the gradient norm explodes exponentially at the rate of B-2 B-3 , which goes to 1 as B → ∞. Fig. 1 shows theory and simulation for ReLU gradient dynamics.

Weight Gradient While all of the above only study the gradient with respect to the hidden preactivations, Appendix L shows that the weight gradient norms at layer l is just Π l , µ * G = µ * tr Π l , and thus by Thm 3.9, the weight gradients explode as well at the same rate λ ↓ G .

Effect of as a hyperparameter In practice, is usually treated as small constant and is not regarded as a hyperparameter to be tuned. Nevertheless, we can investigate its effect on gradient explosion. A straightforward generalization of the analysis presented above to the case of > 0 suggests somewhat larger values than typically used can ameliorate (but not eliminate) gradient explosion problems. See Fig. 5(c, d ).

3.2. CROSS-BATCH DYNAMICS:

Forward In addition to analyzing the correlation between preactivations of samples in a batch, we also study the correlation between those of different batches. The dynamics Eq. ( 8) can be generalized to simultaneous propagation of k batches (see Eq. ( 54) and Appendix H).

EQUATION

)

Here the domain of the dynamics is the space of block matrices, with diagonal blocks and offdiagonal blocks resp. representing within-batch and cross-batch covariance. We observe empirically A linear network with batch norm is shown acting on two minibatches of size 64 after random orthogonal initialization. The datapoints in the minibatch are chosen to form a 2d circle in input space, except for one datapoint that is perturbed separately in each minibatch (leftmost datapoint at input layer 0). Because the network is linear, for a given minibatch it performs an affine transformation on its inputs -a circle in input space remains an ellipse throughout the network. However, due to batch norm the coefficients of that affine transformation change nonlinearly as the datapoints in the minibatch are changed. (a) Each pane shows a scatterplot of activations at a given layer for all datapoints in the minibatch, projected onto the top two PCA directions. PCA directions are computed using the concatenation of the two minibatches. Due to the batch norm nonlinearity, minibatches that are nearly identical in input space grow increasingly dissimilar with depth. Intuitively, this chaotic input-output map can be understood as the source of exploding gradients when batch norm is applied to very deep networks, since very small changes in an input correspond to very large movements in network outputs. Batch norm leads to a chaotic input-output map with increasing depth. A linear batch norm is shown acting on two minibatches of size 64 after random orthogonal . The datapoints in the minibatch are chosen to form a 2d circle in input space, except oint that is perturbed separately in each minibatch (leftmost datapoint at input layer 0). etwork is linear, for a given minibatch it performs an affine transformation on its inputs input space remains an ellipse throughout the network. However, due to batch norm ts of that affine transformation change nonlinearly as the datapoints in the minibatch (a) Each pane shows a scatterplot of activations at a given layer for all datapoints tch, projected onto the top two PCA directions. PCA directions are computed using ation of the two minibatches. Due to the batch norm nonlinearity, minibatches that entical in input space grow increasingly dissimilar with depth. Intuitively, this chaotic map can be understood as the source of exploding gradients when batch norm is applied networks, since very small changes in an input correspond to very large movements in uts. (b) The correlation between the two minibatches, as a function of layer, for the k. Despite having a correlation near one at the input layer, the two minibatches rapidly ith depth. See ?? for a theoretical treatment. The datapoints in the minibatch are chosen to form a 2d circle in input space, except for one datapoint that is perturbed separately in each minibatch (leftmost datapoint at input layer 0). Because the network is linear, for a given minibatch it performs an affine transformation on its inputs -a circle in input space remains an ellipse throughout the network. However, due to batch norm the coefficients of that affine transformation change nonlinearly as the datapoints in the minibatch are changed. (a) Each pane shows a scatterplot of activations at a given layer for all datapoints in the minibatch, projected onto the top two PCA directions. PCA directions are computed using the concatenation of the two minibatches. Due to the batch norm nonlinearity, minibatches that are nearly identical in input space grow increasingly dissimilar with depth. Intuitively, this chaotic input-output map can be understood as the source of exploding gradients when batch norm is applied to very deep networks, since very small changes in an input correspond to very large movements in network outputs. (b) The correlation between the two minibatches, as a function of layer, for the same network. Despite having a correlation near one at the input layer, the two minibatches rapidly decorrelate with depth. See ?? for a theoretical treatment.

A.1 VGG19 WITH BATCHNORM ON CIFAR100

Even though at initialization time batchnorm causes gradient explosion, after the first few epochs, the relative gradient norms kr ✓ Lk/k✓k for weight parameters ✓ = W or BN scale parameter ✓ = , equilibrate to about the same magnitude. See Figure 2 : Batch norm leads to a chaotic input-output map with increasing depth. A linear network with batch norm is shown acting on two minibatches of size 64 after random orthogonal initialization. The datapoints in the minibatch are chosen to form a 2d circle in input space, except for one datapoint that is perturbed separately in each minibatch (leftmost datapoint at input layer 0). Because the network is linear, for a given minibatch it performs an affine transformation on its inputs -a circle in input space remains an ellipse throughout the network. However, due to batch norm the coefficients of that affine transformation change nonlinearly as the datapoints in the minibatch are changed. (a) Each pane shows a scatterplot of activations at a given layer for all datapoints in the minibatch, projected onto the top two PCA directions. PCA directions are computed using the concatenation of the two minibatches. Due to the batch norm nonlinearity, minibatches that are nearly identical in input space grow increasingly dissimilar with depth. Intuitively, this chaotic input-output map can be understood as the source of exploding gradients when batch norm is applied to very deep networks, since very small changes in an input correspond to very large movements in network outputs. (b) The correlation between the two minibatches, as a function of layer, for the same network. Despite having a correlation near one at the input layer, the two minibatches rapidly decorrelate with depth. See Appendix H for a theoretical treatment.

that for most nonlinearities used in practice like ReLU, or even for the identity function, ie no pointwise nonlinearity, the global fixed point of this dynamics is cross-batch BSB1, with diagonal BSB1 blocks, and off-diagonal entries all equal to the same constant c * CB :

EQUATION

To interpret this phenomenon, note that, after mean centering of the preactivations, this covariance matrix becomes a multiple of identity. Thus, deep embedding of batches loses the mutual information between them in the input space. Qualitatively, this implies that two batches that are similar at the input to the network will become increasingly dissimilar -i.e. chaotic -as the signal propagates deep into the network. This loss in fact happens exponentially fast, as illustrated in Fig. 2 , and as shown theoretically by Thm H.13, at a rate of λ ↑ M < 1. Formally, the linearized dynamics of Eq. ( 20) has an eigenspace M given by all block matrices with zero diagonal blocks. The associated eigenvalue is λ ↑ M . We shall return to λ ↓ M shortly and give its Gegenbauer expansion (Thm 3.11) .

Backward The gradients of two batches of input are correlated throughout the course of backpropagation. We hence also study the generalization of Eq. ( 17) to the k-batch setting (see Appendix I):

EQUATION

where for simplicity of exposition and without loss of generality, we take k = 2. Like in Eq. ( 17), we have assumed that Σ l has converged to its limit Σ * . It is not hard to see that each block of Πl evolves independently, with the diagonal blocks in particular evolving according to Eq. ( 17). Analyzing the off-diagonal blocks might seem unwieldy at first, but it turns out that its dynamics is given simply by scalar multiplication by a constant λ ↓ M (Thm I.4). Even more surprisingly, Theorem 3.11 (Batchnorm causes information loss). The convergence rate of cross-batch preactivation covariance λ ↑ M is equal to the decay rate of the cross-batch gradient covariance λ ↓

M:

. If

EQUATION

)

This quantity is always < 1, and for any fixed B, is maximized by φ = id. Furthermore, for φ = id,

EQUATION

-2 and increase to 1 from below as B → ∞.

Thus, a deep batchnorm network loses correlation information between two input batches, exponentially fast in depth, no matter what nonlinearity (that induces fixed point of the form given by Eq. ( 21)). This again contrasts with the case for vanilla networks which can control this rate of information loss by tweaking the initialization variances for weights and biases Poole et al. (2016); Schoenholz et al. (2016) . The absence of coordinatewise nonlinearity, i.e. φ = id, maximally suppresses both this loss of information as well as the gradient explosion, and in the ideal, infinite batch scenario, can cure both problems.

4. EXPERIMENTS:

Having developed a theory for neural networks with batch normalization at initialization, we now explore the relationship between the properties of these random networks and their learning dynamics. We will see that the trainability of networks with batch normalization is controlled by gradient explosion. We quantify the depth scale over which gradients explode by ξ = 1/ log λ ↓ G where, as above, λ ↓ G is the largest eigenvalue of the jacobian. Across many different experiments we will see strong agreement between ξ and the maximum trainable depth.

We first investigate the relationship between trainability and initialization for rectified linear networks as a function of batch size. The results of these experiments are shown in Fig. 3 where in each case we plot the test accuracy after training as a function of the depth and the batch size and overlay 16ξ in white dashed lines. In Fig. 3 (a) we consider networks trained using SGD on MNIST where we observe that networks deeper than about 50 layers are untrainable regardless of batch size. In (b) we compare standard batch normalization with a modified version in which the batch size is held fixed but batch statistics are computed over subsets of size B. This removes subtle gradient fluctuation effects noted in Smith & Le (2018) . In (c) we do the same experiment with RMSProp and in (d) we train the networks on CIFAR10. In all cases we observe a nearly identical trainable region.

t < l a t e x i t s h a 1 _ b a s e 6 4 = " U S o N v 0 f C I v i 1 l L E l 3 u P / p 3 R Z N / c = " > A A A B + H i c b Z A 9 S w N B E I b n 4 l e M X 1 F L m 8 U g W I U 7 E b Q M 2 F g m Y D 4 g O c L e Z i 9 Z s r t 3 7 M 4 J M e Q X 2 G p v J 7 b + G 1 t / i Z v k C k 1 8 Y e D h n R l m e K N U C o u + / + U V N j a 3 t n e K u 6 W 9 / Y P D o / L x S c s m m W G 8 y R K Z m E 5 E L Z d C 8 y Y K l L y T G k 5 V J H k 7 G t / N + + 1 H b q x I 9 A N O U h 4 q O t Q i F o y i s x r Y L 1 f 8 q r 8 Q W Y c g h w r k q v f L 3 7 1 B w j L F N T J J r e 0 G f o r h l B o U T P J Z q Z d Z n l I 2 p k P e d a i p 4 j a c L h 6 d k Q v n D E i c G F c a y c L 9 v T G l y t q J i t y k o j i y q 7 2 5 + W 8 v U i u X M b 4 N p 0 K n G X L N l o f j T B J M y D w F M h C G M 5 Q T B 5 Q Z 4 X 4 n b E Q N Z e i y K r l Q g t U I 1 q F 1 V Q 0 c N 6 4 r N T + P p w h n c A 6 X E M A N 1 O A e 6 t A E B h y e 4 Q V e v S f v z X v 3 P p a j B S / f O Y U / 8 j 5 / A B K 4 k 4 w = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " U S o N v 0 f C I v i 1 l L E l 3 u P / p 3 R Z N / c = " > A A A B + H i c b Z A 9 S w N B E I b n 4 l e M X 1 F L m 8 U g W I U 7 E b Q M 2 F g m Y D 4 g O c L e Z i 9 Z s r t 3 7 M 4 J M e Q X 2 G p v J 7 b + G 1 t / i Z v k C k 1 8 Y e D h n R l m e K N U C o u + / + U V N j a 3 t n e K u 6 W 9 / Y P D o / L x S c s m m W G 8 y R K Z m E 5 E L Z d C 8 y Y K l L y T G k 5 V J H k 7 G t / N + + 1 H b q x I 9 A N O U h 4 q O t Q i F o y i s x r Y L 1 f 8 q r 8 Q W Y c g h w r k q v f L 3 7 1 B w j L F N T J J r e 0 G f o r h l B o U T P J Z q Z d Z n l I 2 p k P e d a i p 4 j a c L h 6 d k Q v n D E i c G F c a y c L 9 v T G l y t q J i t y k o j i y q 7 2 5 + W 8 v U i u X M b 4 N p 0 K n G X L N l o f j T B J M y D w F M h C G M 5 Q T B 5 Q Z 4 X 4 n b E Q N Z e i y K r l Q g t U I 1 q F 1 V Q 0 c N 6 4 r N T + P p w h n c A 6 X E M A N 1 O A e 6 t A E B h y e 4 Q V e v S f v z X v 3 P p a j B S / f O Y U / 8 j 5 / A B K 4 k 4 w = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " U S o N v 0 f C I v i 1 l L E l 3 u P / p 3 R Z N / c = " > A A A B + H i c b Z A 9 S w N B E I b n 4 l e M X 1 F L m 8 U g W I U 7 E b Q M 2 F g m Y D 4 g O c L e Z i 9 Z s r t 3 7 M 4 J M e Q X 2 G p v J 7 b + G 1 t / i Z v k C k 1 8 Y e D h n R l m e K N U C o u + / + U V N j a 3 t n e K u 6 W 9 / Y P D o / L x S c s m m W G 8 y R K Z m E 5 E L Z d C 8 y Y K l L y T G k 5 V J H k 7 G t / N + + 1 H b q x I 9 A N O U h 4 q O t Q i F o y i s x r Y L 1 f 8 q r 8 Q W Y c g h w r k q v f L 3 7 1 B w j L F N T J J r e 0 G f o r h l B o U T P J Z q Z d Z n l I 2 p k P e d a i p 4 j a c L h 6 d k Q v n D E i c G F c a y c L 9 v T G l y t q J i t y k o j i y q 7 2 5 + W 8 v U i u X M b 4 N p 0 K n G X L N l o f j T B J M y D w F M h C G M 5 Q T B 5 Q Z 4 X 4 n b E Q N Z e i y K r l Q g t U I 1 q F 1 V Q 0 c N 6 4 r N T + P p w h n c A 6 X E M A N 1 O A e 6 t A E B h y e 4 Q V e v S f v z X v 3 P p a j B S / f O Y U / 8 j 5 / A B K 4 k 4 w = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " U S o N v 0 f C I v i 1 l L E l 3 u P / p 3 R Z N / c = " > A A A B + H i c b Z A 9 S w N B E I b n 4 l e M X 1 F L m 8 U g W I U 7 E b Q M 2 F g m Y D 4 g O c L e Z i 9 Z s r t 3 7 M 4 J M e Q X 2 G p v J 7 b + G 1 t / i Z v k C k 1 8 Y e D h n R l m e K N U C o u + / + U V N j a 3 t n e K u 6 W 9 / Y P D o / L x S c s m m W G 8 y R K Z m E 5 E L Z d C 8 y Y K l L y T G k 5 V J H k 7 G t / N + + 1 H b q x I 9 A N O U h 4 q O t Q i F o y i s x r Y L 1 f 8 q r 8 Q W Y c g h w r k q v f L 3 7 1 B w j L F N T J J r e 0 G f o r h l B o U T P J Z q Z d Z n l I 2 p k P e d a i p 4 j a c L h 6 d k Q v n D E i c G F c a y c L 9 v T G l y t q J i t y k o j i y q 7 2 5 + W 8 v U i u X M b 4 N p 0 K n G X L N l o f j T B J M y D w F M h C G M 5 Q T B 5 Q Z 4 X 4 n b E Q N Z e i y K r l Q g t U I 1 q F 1 V Q 0 c N 6 4 r N T + P p w h n c A 6 X E M A N 1 O A e 6 t A E B h y e 4 Q V e v S f v z X v 3 P p a j B S / f O Y U / 8 j 5 / A B K 4 k 4 w = < / l a t e x i

It is counter intuitive that training can occur at intermediate depths, from 10 to 50 layers, where there is significant gradient explosion. To gain insight into the behavior of the network during learning we record the magnitudes of the weights, the gradients with respect to the pre-activations, and the gradients with respect to the weights for the first 10 steps of training for networks of different depths. The result of this experiment is shown in Fig. 4 . Here we see that before learning, as expected, the norm of the weights is constant and independent of layer while the gradients feature exponential explosion. However, we observe that two related phenomena occur after a single step of learning: the weights grow exponentially in the depth and the magnitude of the gradients are stable up to some threshold after which they vanish exponentially in the depth. This is as the result of the scaling property of batchnorm, where

B φ (αh) = α -1 B φ (h):

The first-step gradients dominate the weights due to gradient explosion, hence the exponential growth in weight norms, and thereafter, the gradients are scaled down commensurately. Thus, it seems that although the gradients of batch normalized networks at initialization are ill-conditioned, the gradients appear to quickly reach a stable dynamical equilibrium. While this appears to be beneficial for shallower networks, in deeper ones, the relative gradient vanishing can in fact be so severe as to cause lower layers to mostly stay constant during training. Aside from numerical issues, this seems to be the primary mechanism through which gradient explosion causes training problems for networks deeper than 50 layers.

As discussed in the theoretical exposition above, batch normalization necessarily features exploding gradients for any nonlinearity that converges to a BSB1 fixed point. We performed a number of experiments exploring different ways of ameliorating this gradient explosion. These experiments c, d ) tanh and rectified linear networks respectively as a function of , here we observe a well defined phase transition near ∼ 1. Note that in the case of rectified linear activations we use β = 2 so that the function is locally linear about 0.

We also find initializing β and/or setting > 0 having positive effect on VGG19 with batchnorm. See Figs. 8 and 9 .

are shown in Fig. 5 with theoretical predictions for the maximum trainable depth overlaid; in all cases we see exceptional agreement. In Fig. 5 (a, b ) we explore two different ways of tuning the degree to which activation functions in a network are nonlinear. In Fig. 5 (a) we tune γ ∈ [0, 2] for networks with tanh-activations and note that in the γ → 0 limit the function is linear. In Fig. 5 (b) we tune β ∈ [0, 2] for networks with rectified linear activations and we note, similarly, that in the β → ∞ limit the function is linear. As expected, we see the maximum trainable depth increase significantly with decreasing γ and increasing β. In Fig. 5 (c, d ) we vary for tanh and rectified linear networks respectively. In both cases, we observe a critical point at large where gradients do not explode and very deep networks are trainable.

5. CONCLUSION:

In this work we have presented a theory for neural networks with batch normalization at initialization. In the process of doing so, we have uncovered a number of counterintuitive aspects of batch normalization and -in particular -the fact that at initialization it unavoidably causes gradients to explode with depth. We have introduced several methods to reduce the degree of gradient explosion, enabling the training of significantly deeper networks in the presence of batch normalization. Finally, this work paves the way for future work on more advanced, state-of-the-art, network architectures and topologies.

(a) 

Correlation Correlation between minibatches:

Figure 6 : Batch norm leads to a chaotic input-output map with increasing depth. A linear network with batch norm is shown acting on two minibatches of size 64 after random orthogonal initialization. The datapoints in the minibatch are chosen to form a 2d circle in input space, except for one datapoint that is perturbed separately in each minibatch (leftmost datapoint at input layer 0). Because the network is linear, for a given minibatch it performs an affine transformation on its inputs -a circle in input space remains an ellipse throughout the network. However, due to batch norm the coefficients of that affine transformation change nonlinearly as the datapoints in the minibatch are changed. (a) Each pane shows a scatterplot of activations at a given layer for all datapoints in the minibatch, projected onto the top two PCA directions. PCA directions are computed using the concatenation of the two minibatches. Due to the batch norm nonlinearity, minibatches that are nearly identical in input space grow increasingly dissimilar with depth. Intuitively, this chaotic input-output map can be understood as the source of exploding gradients when batch norm is applied to very deep networks, since very small changes in an input correspond to very large movements in network outputs. (b) The correlation between the two minibatches, as a function of layer, for the same network. Despite having a correlation near one at the input layer, the two minibatches rapidly decorrelate with depth. See Appendix H for a theoretical treatment.

A VGG19 WITH BATCHNORM ON CIFAR100

Even though at initialization time batchnorm causes gradient explosion, after the first few epochs, the relative gradient norms ∇ θ L / θ for weight parameters θ = W or BN scale parameter θ = γ, equilibrate to about the same magnitude. See Fig. 7 . with batchnorm on CIFAR100 with data augmentation. We use 8 random seeds for each combination, and assign to each combination the median training/validation accuracy over all runs. We then aggregate these scores here. In the first row we look at training accuracy with different learning rate vs β initialization at different epochs of training, presenting the max over . In the second row we do the same for validation accuracy. In the third row, we look at the matrix of training accuracy for learning rate vs , taking max over β. In the fourth row, we do the same for validation accuracy.

B GRADIENT INDEPENDENCE ASSUMPTION:

Following prior literature Schoenholz et al. ( 2016); Yang & Schoenholz (2017) ; Xiao et al. (2018) , in this paper, in regards to computations involving backprop, we assume Assumption 2. During backpropagation, whenever we multiply by W T for some weight matrix W , we multiply by an iid copy instead. 8 , except we don't take the max over the unseen hyperparameter but rather set it to 0 (the default value).

As in these previous works, we find excellent agreement between computations made under this assumption and the simulations (see Fig. 1 ). Yang (2019) in fact recently rigorously justified this assumption as used in the computation of moments, for a wide variety of architectures, like multilayer perceptron (Schoenholz et al., 2016) , residual networks (Yang & Schoenholz, 2017) , and convolutional networks (Xiao et al., 2018) studied previously, but without batchnorm. The reason that their argument does not extend to batchnorm is because of the singularity in its Jacobian at 0. However, as observed in our experiments, we expect that a proof can be found to extend Yang & Schoenholz (2017) 's work to batchnorm.

C NOTATIONS:

Definition C.1. Let S B be the space of PSD matrices of size B × B. Given a measurable function Φ : R B → R A , define the integral transform V Φ : S B → S A by V Φ (Σ) = E[Φ(h) ⊗2 : h ∼ N (0, Σ)].

3 When φ : R → R and B is clear from context, we also write V φ for V applied to the function acting coordinatewise by φ.

Definition C.2. For any φ : R → R, let B φ : R B → R B be batchnorm (applied to a batch of neuronal activations) followed by coordinatewise applications of φ, When B is clear from context, we will suppress the subscript/superscript B. In short, for h ∈ R B , Gh zeros the sample mean of h. G is a projection matrix to the subspace of vectors h ∈ R B of zero coordinate sum. With the above definitions, we then have B φ (h) = φ( √ BGh/ Gh ).

B φ (x) j = φ (x j -Avg x) / 1 B B i=1 (x i -Avg x) 2 (here Avg x = 1 B B i=1 x i ). When φ = id

We use Γ to denote the Gamma function, P (a, b) := Γ(a + b)/Γ(a) to denote the Pochhammer symbol, and Beta(a, b) = Γ(a)Γ(b) Γ(a+b) to denote the Beta function. We use •, • to denote the dot product for vectors and trace inner product for matrices. Sometimes when it makes for simpler notation, we will also use • for dot product.

We adopt the matrix convention that, for a multivariate function f : R n → R m , the Jacobian is 

df dx =       ∂f1 ∂x1 ∂f1 ∂x2 • • • ∂f1 ∂xn ∂f2 ∂x1 ∂f2 ∂x2

D A GUIDE TO THE REST OF THE APPENDIX:

As discussed in the main text, we are interested in several closely related dynamics induced by batchnorm in a fully-connected network with random weights. One is the forward propagation equation Eq. ( 43)

Σ l = V B φ (Σ l-1 )

studied in Appendix F. Another is the backward propagation equation Eq. ( 51)

Π l = V B φ (Σ * ) † {Π l+1 }

studied in Appendix G. We will also study their generalizations to the simultaneous propagation of multiple batches in Appendices H and I.

In general, we discover several ways to go about such analyses, each with their own benefits and drawbacks:

• The Laplace Method (Appendix E.1): If the nonlinearity φ is positive homogeneous, we can use Schwinger's parametrization to turn V B φ and V B φ into nicer forms involving only 1 or 2 integrals, as mentioned in the main text. All relevant quantities, such as eigenvalues of the backward dynamics, can be further simplified to closed forms. Lemma E.2 gives the master equation for the Laplace method. Because the Laplace method requires φ to be positive homogeneous to apply, we review relevant results of such functions in Appendices E.1.1 and E.1.2.

• The Fourier Method (Appendix E.2): For polynomially bounded and continuous nonlinearity φ (for the most general set of conditions, see Thm E.25), we can obtain similar simplifications using the Fourier expansion of the delta function, with the penalty of an additional complex integral.

• Spherical Integration (Appendix E.3): Batchnorm can naturally be interpreted as a linear projection followed by projection to a sphere of radius √ B. Assuming that the forward dynamics converges to a BSB1 fixed point, one can express V B φ and V B φ as spherical integrals, for all measurable φ. One can reduce this (B -1)-dimensional integral, in spherical angular coordinates, into 1 or 2 dimensions by noting that the integrand depending on φ only depends on 1 or 2 angles. Thus this method is very suitable for numerical evaluation of quantities of interest for general φ.

• The Gegenbauer Method (Appendix E.4): Gegenbauer polynomials are orthogonal polynomials under the weight (1x 2 ) α for some α. They correspond to a special type of spherical harmonics called zonal harmonics (Defn E.44) that has the very useful reproducing property (Fact E.46). By expressing the nonlinearity φ in this basis, we can see easily that the relevant eigenvalues of the forward and backward dynamics are ratios of two quadratic forms of φ's Gegenbauer coefficients. In particular, for the largest eigenvalue of the backward dynamics, the two quadratic forms are both diagonal, and in a way that makes apparent the necessity of gradient explosion (under the BSB1 fixed point assumption). In regards to numerical computations, the Gegenbauer method has the benefit of only requiring 1-dimensional integrals to obtain the coefficients, but if φ has slowlydecaying coefficients, then a large number of such integrals may be required to get an accurate answer.

In each of the following sections, we will attempt to conduct an analysis with each method when it is feasible.

When studying the convergence rate of the forward dynamics and the gradient explosion of the backward dynamics, we encounter linear operators on matrices that satisfy an abundance of symmetries, called ultrasymmetric operators. We study these operators in Appendix E.5 and contribute several structural results Thms E.61, E.62, E.72 and E.72 on their eigendecomposition that are crucial to deriving the asymptotics of the dynamics.

All of the above will be explained in full in the technique section Appendix E. We then proceed to study the forward and backward dynamics in detail (Appendices F to I), now armed with tools we need. Up to this point, we only consider the dynamics assuming Eq. ( 43) converges to a BSB1 fixed point, and B ≥ 4. In Appendix J we discuss what happens outside this regime, and in Appendix K we show our current understanding of BSB2 fixed point dynamics.

Note that the backward dynamics studied in these sections is that of the gradient with respect to the hidden preactivations. Nevertheless, we can compute the moments of the weight gradients from this, which is of course what is eventually used in gradient descent. We do so in the final section Appendix L.

Main Technical Results Corollary F.3 establishes the global convergence of Eq. ( 43) to BSB1 fixed point for φ = id (as long as the initial covariance is nondegenerate), but we are not able to prove such a result for more general φ. We thus resort to studying the local convergence properties of BSB1 fixed points.

First, Thms F.5, F.8 and F.13 respectively compute the BSB1 fixed point from the perspectives of spherical integration, the Laplace method, and the Gegenbauer method. Thms F.22 and F.24 give the Gegenbauer expansion of the BSB1 local convergence rate, and Thm F.33 gives a more succinct form of the same for positive homogeneous φ. Appendix F.3.1 discusses how to do this computation via spherical integration.

Next, we study the gradient dynamics, and Thms G.5 and G.12 yield the gradient explosion rates with the Gegenbauer method and the Laplace method (for positive homogeneous φ), respectively. Appendix G.1.1 discusses how to compute this for general φ using spherical integration.

We then turn to the dynamics of cross batch covariances. Thm H.3 gives the form of the cross batch fixed point via spherical integration and the Gegenbauer method, and Thm H.6 yields a more specific form for positive homogeneous φ via the Laplace method. The local convergence rate to this fixed point is given by Thms H.13 and H.14, respectively using the Gegenbauer method and the Laplace method. Finally, the correlation between gradients of the two batches is shown to decrease exponentially in Corollary I.6 using the Gegenbauer method, and the decay rate is given succinctly for positive homogeneous φ in Thm I.3.

E TECHNIQUES:

We introduce the key techniques and tools in our analysis in this section.

E.1 LAPLACE METHOD:

As discussed above, the Laplace method is useful for deriving closed form expressions for positive homogeneous φ. The key insight here is to apply Schwinger parametrization to deal with normalization.

Lemma E.1 (Schwinger parametrization). For z > 0 and c > 0,

c -z = Γ(z) -1 ∞ 0 x z-1 e -cx dx

The following is the key lemma in the Laplace method.

Lemma E.2 (The Laplace Method Master Equation).

For A, A ∈ N, let f : R A → R A and let k ≥ 0. Suppose f (y) ≤ h( y ) for some nondecreasing function h : R ≥0 → R ≥0 such that E[h(r y ) : y ∈ N (0, I A )] exists for every r ≥ 0. Define ℘(Σ) := E[ y -2k f (y) : y ∼ N (0, Σ)].

Then on {Σ ∈ S A : rank Σ > 2k}, ℘(Σ) is well-defined and continuous, and furthermore satisfies

EQUATION

Proof. ℘(Σ) is well-defined for full rank Σ because the y -2k singularity at y = 0 is Lebesgueintegrable in a neighborhood of 0 in dimension A > 2k.

We prove Eq. ( 25) in the case when Σ is full rank and then apply a continuity argument.

Proof of Eq. (25) for full rank Σ. First, we will show that we can exchange the order of integration

∞ 0 ds R A dy s k-1 f (y)e -1 2 y T (Σ -1 +2sI)y = R A dy ∞ 0 ds s k-1 f (y)e -1 2 y T (Σ -1 +2sI)y

by Fubini-Tonelli's theorem. Observe that

(2π) -A/2 det Σ -1/2 ∞ 0 ds s k-1 R A dy f (y) e -1 2 y T (Σ -1 +2sI)y = ∞ 0 ds s k-1 det(Σ(Σ -1 + 2sI)) -1/2 E[ f (y) : y ∼ N (0, (Σ -1 + 2sI) -1 )] = ∞ 0 ds s k-1 det(I + 2sΣ) -1/2 E[ f (y) : y ∼ N (0, Σ(I + 2sΣ) -1 )] For λ = Σ 2 , f ( Σ(I + 2sΣ) -1 y) ≤ h( Σ(I + 2sΣ) -1 y ≤ h λ 1 + 2sλ y ≤ h( √ λ y ) Because E[h( √ λ y ) : y ∈ N (0, I)] exists by assumption, E[ f (y) : y ∼ N (0, Σ(I + 2sΣ) -1 )] = E[ f ( Σ(I + 2sΣ) -1 y) : y ∼ N (0, I)] → E[ f (0) : y ∼ N (0, I)] = f (0)

as s → ∞, by dominated convergence with dominating function h( √ λ y )e -1 2 y 2 (2π) -A/2 . By the same reasoning, the function s → E[ f (y) : y ∼ N (0, Σ(I + 2sΣ) -1 )] is continuous. In particular this implies that sup 0≤s≤∞ E[ f (y) : y ∼ N (0, Σ(I + 2sΣ) -1 )] < ∞. Combined with the fact that det(I + 2sΣ

) -1/2 = Θ(s -A/2 ) as s → ∞, ∞ 0 ds s k-1 det(I + 2sΣ) -1/2 E[ f (y) : y ∼ N (0, Σ(I + 2sΣ) -1 )] ≤ ∞ 0 ds Θ(s k-1-A/2 )

which is bounded by our assumption that A/2 > k. This shows that we can apply Fubini-Tonelli's theorem to allow exchanging order of integration.

Thus,

E[ y -2k f (y) : y ∼ N (0, Σ)] = E[f (y) ∞ 0 ds Γ(k) -1 s k-1 e -y 2 s : y ∼ N (0, Σ)] = (2π) -A/2 det Σ -1/2 dy e -1 2 y T Σ -1 y f (y) ∞ 0 ds Γ(k) -1 s k-1 e -y 2 s = (2π) -A/2 Γ(k) -1 det Σ -1/2 ∞ 0 ds s k-1 dy f (y)e -1 2 y T (Σ -1 +2sI)y (by Fubini-Tonelli) = Γ(k) -1 ∞ 0 ds s k-1 det(Σ(Σ -1 + 2sI)) -1/2 E[f (y) : y ∼ N (0, (Σ -1 + 2sI) -1 )] = Γ(k) -1 ∞ 0 ds s k-1 det(I + 2sΣ) -1/2 E[f (y) : y ∼ N (0, Σ(I + 2sΣ) -1 )]

Domain and continuity of ℘(Σ). The LHS of Eq. ( 25), ℘(Σ), is defined and continuous on rank Σ/2 > k.

Indeed, if Σ = M I C M T , where M is a full rank A × C matrix with rank Σ = C ≤ A, then E[ y -2k f (y) : y ∼ N (0, Σ)] = E[ M z -2k f (M z) : z ∼ N (0, I C )].

This is integrable in a neighborhood of 0 iff C > 2k, while it's always integrable outside a ball around 0 because f by itself already is. So ℘(Σ) is defined whenever rank Σ > 2k. Its continuity can be established by dominated convergence.

Proof of Eq. (25) for rank Σ > 2k. Observe that det(I + 2sΣ) -1/2 is continuous in Σ and, by an application of dominated convergence as in the above, E[f (y) : y ∼ N (0, Σ(I + 2sΣ) -1 )] is continuous in Σ. So the RHS of Eq. ( 25) is continuous in Σ whenever the integral exists. By the reasoning above, E[f (y) : y ∼ N (0, Σ(I + 2sΣ) -1 )] is bounded in s and det(I + 2sΣ) -1/2 = Θ(s -rank Σ/2 ), so that the integral exists iff rank Σ/2 > k.

To summarize, we have proved that both sides of Eq. ( 25) are defined and continous for rank Σ > 2k. Because the full rank matrices are dense in this set, by continuity Eq. ( 25) holds for all rank Σ > 2k.

If φ is degree-α positive homogeneous, i.e. φ(ru) = r α φ(u) for any u ∈ R, r ∈ R + , we can apply Lemma E.2 with k = α, (Cho & Saul, 2009) .

V B φ (Σ) = E[φ( √ BGh/ Gh ) ⊗2 : h ∼ N (0, Σ)] = B α E[φ(Gh) ⊗2 / Gh 2α : h ∼ N (0, Σ)] = B α E[φ(z) ⊗2 / z 2α : z ∼ N (0, GΣG)] = B α Γ(α) -1 ∞ 0 ds s α-1 det(I + 2sΣG) -1/2 V φ (G(I + 2sΣG) -1 ΣG). If φ = relu, then V φ (Σ) ij = 1 2 Σ ii if i = j 1 2 J 1 (Σ ij / Σ ii Σ jj ) Σ ii Σ jj otherwise and, more succinctly, V φ (Σ) = D 1/2 V φ (D -1/2 ΣD -1/2 )D 1/2 where D = Diag(Σ). Here J 1 (c) := 1 π ( √ 1 -c 2 + (π -arccos(c))c)

Matrix simplification. We can simplify the expression G(I + 2sΣG) -1 ΣG, leveraging the fact that G is a projection matrix.

Definition E.3. Let e be an B × (B -1) matrix whose columns form an orthonormal basis of

im G := {Gv : v ∈ R B } = {w ∈ R B : i w i = 0}.

Then the B × B matrix ẽ = (e|B -1/2 1) is an orthogonal matrix. For much of this paper e can be any such basis, but at certain sections we will consider specific realizations of e for explicit computation.

From easy computations it can be seen that G = ẽ

I B-1 0 0 0 ẽT . Suppose Σ = ẽ Σ v v T a ẽT

where Σ is (B -1) × (B -1), v is a column vector and a is a scalar. Then Σ = e T Σe and

ΣG = ẽ Σ 0 v T 0

ẽT is block lower triangular, and

(I + 2sΣG) -1 = ẽ (I + 2sΣ ) -1 0 * 1 ẽT (I + 2sΣG) -1 ΣG = ẽ (I + 2sΣ ) -1 Σ 0 * 0 ẽT G(I + 2sΣG) -1 ΣG = ẽ (I + 2sΣ ) -1 Σ 0 0 0 ẽT = e(I + 2sΣ ) -1 Σ e T = GΣG(I + 2sGΣG) -1 =: Σ G (I + 2sΣ G ) -1

where Definition E.4. For any matrix Σ, write Σ G := GΣG.

Similarly, det(I B + 2sΣG) = det(I B-1 + 2sΣ ) = det(I B + 2sΣ G ). So, altogether, we have Theorem E.5. Suppose φ : R → R is degree-α positive homogeneous. Then for any B × (B -1) matrix e whose columns form an orthonormal basis of im

EQUATION

Now, in order to use Laplace's method, we require φ to be positive homogeneous. What do these functions look like? The most familiar example to a machine learning audience is most likely ReLU. It turns out that 1-dimensional positive homogeneous functions can always be described as a linear combination of powers of ReLU and its reflection across the y-axis. Below, we review known facts about these functions and their integral transforms, starting with the powers of ReLU in Appendix E.1.1 and then onto the general case in Appendix E.1.2.

E.1.1 α-RELU:

Recall that α-ReLUs (Yang & Schoenholz, 2017) are, roughly speaking, the αth power of ReLU. Definition E.6. The α-ReLU function ρ α : R → R sends x → x α when x > 0 and x → 0 otherwise. This is a continuous function for α > 0 but discontinuous at 0 for all other α.

We briefly review what is currently known about the V and W transforms of ρ α (Cho & Saul, 2009; Yang & Schoenholz, 2017) .

Definition E.7. For any α > -1 2 , define c α = 1 √ π 2 α-1 Γ(α + 1 2 ).

Published as a conference paper at ICLR 2019 When considering only 1-dimensional Gaussians, V ρα is very simple. Proposition E.8. If α > -1 2 , then for any q ∈ S 1 = R ≥0 , V ρα (q) = c α q α To express results of V ρα on S B for higher B, we first need the following Definition E.9. Define

J α (θ) := 1 2πc α (sin θ) 2α+1 Γ(α + 1) π/2 0 dη cos α η (1 -cos θ cos η) 1+α and J α (c) = J α (arccos c) for α > -1/2.

Then Proposition E.10. For any Σ ∈ S B , let D be the diagonal matrix with the same diagonal as Σ. Then

V ρα (Σ) = c α D α/2 J α (D -1/2 ΣD -1/2 )D α/2

where J α is applied entrywise.

For example, J α and J α for the first few integral α are

J 0 (θ) = π -θ π J 1 (θ) = sin θ + (π -θ) cos θ π J 2 (θ) = 3 sin θ cos θ + (π -θ)(1 + 2 cos 2 θ) 3π J 0 (c) = π -arccos c π J 1 (c) = √ 1 -c 2 + (π -arccos c)c π J 2 (c) = 3c √ 1 -c 2 + (π -arccos c)(1 + 2c 2 ) 3π

One can observe very easily that (Daniely et al., 2016; Yang & Schoenholz, 2017 ) Proposition E.11. For each α > -1/2, J α (c) is an increasing and convex function on c ∈ [0, 1],

and is continous on c ∈ [0, 1] and smooth on c ∈ (0, 1). J α

(1) = 1, J α (0) = 1 2 √ π Γ( α 2 + 1 2 ) 2 Γ(α+ 1

2 ) , and J α (-1) = 0.

Yang & Schoenholz (2017) also showed the following fixed point structure Theorem E.12. For α ∈ [1/2, 1), J α (c) = c has two solutions: an unstable solution at 1 ("unstable" meaning J α (1) > 1) and a stable solution in c * ∈ (0, 1) ("stable" meaning J α (c * ) < 1).

The α-ReLUs satisfy very interesting relations amongst themselves. For example, Lemma E.13. Suppose α > 1. Then

J α (θ) = cos θJ α-1 (θ) + (α -1) 2 (2α -1) -1 (2α -3) -1 sin 2 θJ α-2 (θ) J α (c) = cJ α-1 (c) + (α -1) 2 (2α -1) -1 (2α -3) -1 (1 -c 2 )J α-2 (c)

In additional, surprisingly, one can use differentiation to go from α to α + 1 and from α to α -1! Proposition E.14 (Yang & Schoenholz (2017) ). Suppose α > 1/2. Then

J α (θ) = -α 2 (2α -1) -1 J α-1 (θ) sin θ J α (c) = α 2 (2α -1) -1 J α-1 (c) so that J α-n (c) =   α β=α-n+1 β -2 (2β -1)   (∂/∂c) n J α (c)

We have the following from Cho & Saul (2009) Proposition E.15 (Cho & Saul (2009) ). For all α ≥ 0 and integer n ≥ 1

J n+α (c) = c α c n+α (1 -c 2 ) n+α+ 1 2 (∂/∂c) n (J α (c)/(1 -c 2 ) α+ 1 2 ) =   α+n-1 β=α (2β + 1)   -1 (1 -c 2 ) n+α+ 1 2 (∂/∂c) n (J α (c)/(1 -c 2 ) α+ 1 2 )

This implies in particular that we can obtain J α from J α and J α+1 . Proposition E.16. For all α ≥ 0,

J α (c) = (2α + 1)(1 -c 2 ) -1 (J 1+α (c) -cJ α (c)))

Proof.

J 1+α (c) = c α c 1+α (1 -c 2 ) 1+α+ 1 2 (∂/∂c)(J α (c)/(1 -c 2 ) α+ 1 2 ) = (2α + 1) -1 (1 -c 2 ) α+3/2 (J α (c)/(1 -c 2 ) α+1/2 + 2c(α + 1/2)J α (c)/(1 -c 2 ) α+3/2 ) = (2α + 1) -1 J α (c)(1 -c 2 ) + cJ α (c) J α (c) = (2α + 1)(1 -c 2 ) -1 (J 1+α (c) -cJ α (c))

Note that we can also obtain this via Lemma E.13 and Proposition E.14.

E.1.2 POSITIVE-HOMOGENEOUS FUNCTIONS IN 1 DIMENSION:

Suppose for some α ∈ R, φ : R → R is degree α positive-homogeneous, i.e. φ(rx) = r α φ(x) for any x ∈ R, r > 0. The following simple lemma says that we can always express φ as linear combination of powers of α-ReLUs. Proposition E.17. Any degree α positive-homogeneous function φ : R → R with φ(0) = 0 can be written as x → aρ α (x)bρ α (-x).

Proof. Take a = φ(1) and b = φ(-1). Then positive-homogeneity determines the value of φ on R \ {0} and it coincides with x → aρ α (x)bρ α (-x).

As a result we can express the V and W transforms of any positive-homogeneous function in terms of those of α-ReLUs. Proposition E.18. Suppose φ : R → R is degree α positive-homogeneous. By Proposition E.17, φ restricted to R \ {0} can be written as x → aρ α (x)bρ α (-x) for some a and b. Then for any PSD

2 × 2 matrix M , V φ (M ) 11 = (a 2 + b 2 )V ρα (M ) 11 = c α (a 2 + b 2 )M α 11 V φ (M ) 22 = (a 2 + b 2 )V ρα (M ) 22 = c α (a 2 + b 2 )M α 22 V φ (M ) 12 = V φ (M ) 21 = (a 2 + b 2 )V ρα (M ) 12 -2abV ρα (M ) 12 = c α (M 11 M 22 ) α/2 ((a 2 + b 2 )J α (M 12 / M 11 M 22 ) -2abJ α (-M 12 / M 11 M 22 ))

where

M := 1 0 0 -1 M 1 0 0 -1 .

Proof. We directly compute, using the expansion of φ into ρ α s:

V φ (M ) 11 = E[φ(x) 2 : x ∼ N (0, M 11 )] = E[a 2 ρ α (x) 2 + b 2 ρ α (-x) 2 + 2abρ α (x)ρ α (-x)] = a 2 E[ρ α (x) 2 ] + b 2 E[ρ α (-x) 2 ] = (a 2 + b 2 ) E[ρ α (x) 2 : x ∼ N (0, M 11 )] (28) = c α (a 2 + b 2 )M α 11

where in Eq. ( 28) we used negation symmetry of centered Gaussians. The case of V φ (M ) 22 is similar.

V

φ (M ) 12 = E[φ(x)φ(y) : (x, y) ∼ N (0, M )] = E[a 2 ρ α (x)ρ α (y) + b 2 ρ α (-x)ρ α (-y) -abρ α (x)ρ α (-y) -abρ α (-x)ρ α (y)] = (a 2 + b 2 )V ρα (M ) 12 -2abV ρα (M ) 12 = c α (M 11 M 22 ) α/2 ((a 2 + b 2 )J α (M 12 / M 11 M 22 ) -2abJ α (-M 12 / M 11 M 22 ))

where in the last equation we have applied Proposition E.10.

This then easily generalizes to PSD matrices of arbitrary dimension: Corollary E.19. Suppose φ : R → R is degree α positive-homogeneous. By Proposition E.17, φ restricted to R \ {0} can be written as x → aρ α (x)bρ α (-x) for some a and b. Let Σ ∈ S B . Then

V φ (Σ) = c α D α/2 J φ (D -1/2 ΣD -1/2 )D α/2

where J φ is defined below and is applied entrywise. Explicitly, this means that for all i,

V φ (Σ) ii = c α Σ α ii J φ (1) V φ (Σ) ij = c α J φ (Σ ij / Σ ii Σ jj )Σ α/2 ii Σ α/2 jj

Definition E.20. Suppose φ : R → R is degree α positive-homogeneous. By Proposition E.17, φ restricted to R \ {0} can be written as x → aρ α (x)bρ α (-x) for some a and b. Define

J φ (c) := (a 2 + b 2 )J α (c) -2abJ α (-c).

Let us immediately make the following easy but important observations.

Proposition E.21. J φ (1) = a 2 + b 2 and J φ (0) = (a 2 + b 2 -2ab)J α (0) = (a -b) 2 1 2 √ π Γ( α 2 + 1 2 ) 2 Γ(α+ 1 2 ) ,

Proof. Use the fact that J α (-1) = 0 and J α (0

) = 1 2 √ π Γ( α 2 + 1 2 ) 2 Γ(α+ 1

2 ) by Proposition E.11.

As a sanity check, we can easily compute that J id (c) = 2J 1 (c) -2J 1 (-c) = 2c because id(x) = relu(x)-relu(-x). By Corollary E.19 and c 1 = 1 2 , this recovers the obvious fact that V (id) (Σ) = Σ. We record the partial derivatives of V φ . Proposition E.22. Let φ be positive homogeneous of degree α. Then for all i with Σ ii = 0,

∂V φ (Σ) ii ∂Σ ii = c α αΣ α-1 ii J φ (1)

For all i = j with Σ ii , Σ jj = 0,

∂V φ (Σ) ii ∂Σ ij = 0 ∂V φ (Σ) ij ∂Σ ij = c α Σ (α-1)/2 ii Σ (α-1)/2 jj J φ (c ij ) ∂V φ (Σ) ij ∂Σ ii = 1 2 c α Σ 1 2 (α-2) ii Σ 1 2 α jj (αJ φ (c ij ) -c ij J φ (c ij ))

where c ij = Σ ij / Σ ii Σ jj and J φ denotes its derivative.

Proposition E.23. If φ(c) = aρ α (c) -bρ α (-c) with α > 1/2 on R \ {0}, then J φ (c) = (2α -1) -1 J φ (c) = (2α + 1)(1 -c 2 ) -1 ((a 2 + b 2 )J α+1 (c) + 2abJ α+1 (-c)) -cJ φ (c) J φ (0) = (a + b) 2 1 √ π Γ( α 2 + 1) 2 Γ(α + 1 2 )

Proof. We have φ (c) = aαρ α-1 (c) + bαρ α-1 (-c). On the other hand,

J φ = (a 2 + b 2 )J α (c) + 2abJ α (-c) which by Proposition E.14 is α 2 (2α -1) -1 ((a 2 + b 2 )J α-1 (c) + 2abJ α-1 (-c)) = (2α -1) -1 J φ (c).

This proves the first equation.

With Proposition E.16,

J φ (c) = (a 2 + b 2 )J α (c) + 2abJ α (-c) = (a 2 + b 2 )(2α + 1)(1 -c 2 ) -1 (J α+1 (c) -cJ α (c)) + 2ab(2α + 1)(1 -c 2 ) -1 (J α+1 (-c) + cJ α (-c)) = (2α + 1)(1 -c 2 ) -1 ((a 2 + b 2 )J α+1 (c) + 2abJ α+1 (-c)) -cJ φ (c)

This gives the second equation.

Expanding J α+1 (0) With Proposition E.11, we get

J φ (0) = (2α + 1)((a 2 + b 2 )J α+1 (0) + 2abJ α+1 (0)) = (2α + 1)(a + b) 2 1 2 √ π Γ( α 2 + 1) 2 Γ(α + 3

2 ) Unpacking the definition of Γ(α + 3

2 ) then yields the third equation.

In general, we can factor diagonal matrices out of V φ .

Proposition E.24. For any Σ ∈ S B , D any diagonal matrix, and φ positive-homogeneous with degree α,

V φ (DΣD) = D α V φ (Σ)D α

Proof. For any i,

V φ (DΣD) ii = E[φ(x) 2 : x ∼ N (0, D 2 ii Σ ii )] = E[φ(D ii x) 2 : x ∼ N (0, Σ ii )] = E[D 2α ii φ(x) 2 : x ∼ N (0, Σ ii )] = D 2α ii V φ (Σ) ii

For any i = j,

V φ (DΣD) ij = E[φ(x)φ(y) : (x, y) ∼ N (0, D 2 ii Σ ii D ii D jj Σ ij D ii D jj Σ ji D 2 jj Σ jj ) = E[φ(D ii x)φ(D jj y) : (x, y) ∼ N (0, Σ ii Σ ij Σ ji Σ jj )] = D α ii D α jj E[φ(x)φ(y) : (x, y) ∼ N (0, Σ ii Σ ij Σ ji Σ jj )] = D α ii D α jj V φ (Σ) ij E.2 FOURIER METHOD

The Laplace method crucially used the fact that we can pull out the norm factor Gh out of φ, so that we can apply Schwinger parametrization. For general φ this is not possible, but we can apply some wishful thinking and proceed as follows

V B φ (Σ) = E[φ( √ BGh/ Gh ) ⊗2 : h ∼ N (0, Σ)] = E[φ(Gh/r) ⊗2 ∞ 0 d(r 2 )δ(r 2 -Gh 2 /B) : h ∼ N (0, Σ)] = E[φ(Gh/r) ⊗2 ∞ 0 d(r 2 ) ∞ -∞ dλ 2π e iλ(r 2 -Gh 2 /B) : h ∼ N (0, Σ)] " = " R B dh det(2πΣ) e -1 2 h T Σ -1 h φ(Gh/r) ⊗2 ∞ 0 d(r 2 ) ∞ -∞ dλ 2π e iλ(r 2 -Gh 2 /B) " = " ∞ 0 d(r 2 ) ∞ -∞ dλ 2π R B dh det(2πΣ) e -1 2 h T Σ -1 h φ(Gh/r) ⊗2 e iλ(r 2 -Gh 2 /B) (29) = 1 det(2πΣ) ∞ 0 d(r 2 ) ∞ -∞ dλ 2π e iλr 2 R B dhφ(Gh/r) ⊗2 e -1 2 (h T Σ -1 h+2iλ Gh 2 /B) (30) = 1 det(2πΣ) ∞ 0 d(r 2 ) ∞ -∞ dλ 2π e iλr 2 R B dhφ(Gh/r) ⊗2 e -1 2 (h T (Σ -1 +2iλG/B)h) " = " 1 det(2πΣ) ∞ 0 d(r 2 ) ∞ -∞ dλ 2π e iλr 2 det(2π(Σ -1 + 2iλG/B) -1 ) E[φ(Gh/r) ⊗2 : h ∼ N (0, Σ -1 + 2iλG/B)] = ∞ 0 d(r 2 ) ∞ -∞ dλ 2π e iλr 2 det(I + 2iλGΣ/B) -1/2 E h∼N (0,G(Σ -1 +2iλG/B)Gr -2 )) φ(h) ⊗2 (31) = ∞ 0 d(r 2 ) ∞ -∞ dλ 2π e iλr 2 det(I + 2iλGΣ/B) -1/2 V φ (G(Σ -1 + 2iλG/B)Gr -2 )

Here all steps are legal except for possibly Eq. ( 29), Eq. ( 30), and Eq. ( 31). In Eq. ( 29), we "Fourier expanded" the delta function -this could possibly be justified as a principal value integral, but this cannot work in combination with Eq. ( 30) where we have switched the order of integration, integrating over h first. Fubini's theorem would not apply here because e iλ(r 2 -Gh 2 /B) has norm 1 and is not integrable. Finally, in Eq. ( 31), we need to extend the definition of Gaussian to complex covariance matrices, via complex Gaussian integration. This last point is no problem, and we will just define

E h∼N (0,Σ) f (h) := det(2πΣ) -1/2 R B dh f (h)e -1 2 h T Σ -1 h

for general complex Σ, whenever Σ is nonsingular and this integral is exists, and

E h∼N (0,GΣG) f (h) := det(2π(e T Σe)) -1/2 R B-1 dz f (ez)e -1 2 z T (e T Σe) -1 z

if e T Σe is nonsingular (as in the case above). The definition of V φ stays the same given the above extension of the definition of Gaussian expectation.

Nevertheless, the above derivation is not correct mathematically due to the other points. However, its end result can be justified rigorously by carefully expressing the delta function as the limit of mollifiers. 2. there exists an > 0 such that any γ ∈ [-, ] and for each a, b ∈ [D],

z 2 ∈[a,b] dz e -1 2 z T Σ -1 z |φ(z a / z 2 + γ)φ(z b / z 2 + γ)|

exists and is uniformly bounded by some number possibly depending on a and b.

3. for each a, b ∈ [D], R D dz b a d(r 2 )e -1 2 z T Σ -1 z |φ(z a /r)φ(z b /r)| exists and is finite. 4. ∞ -∞ dλ det(I + 2iλΣ) -1/2 V φ ((Σ -1 + 2iλI) -1 r -2

) exists and is finite for each r > 0.

If e -1 2 z T Σ -1 z φ(z/ z ) ⊗2 is integrable over z ∈ R D , then (2π) -D/2 det Σ -1/2 dz e -1 2 z T Σ -1 z φ(z/ z ) ⊗2 = 1 2π lim a 0,b ∞ b a ds ∞ -∞ dλ e iλs det(I + 2iλΣ) -1/2 V φ ((Σ -1 + 2iλI) -1 s -1 ).

Similarly,

(2π) -D/2 det Σ -1/2 dz e -1 2 z T Σ -1 z φ( √ Dz/ z ) ⊗2 = 1 2π lim a 0,b ∞ b a ds ∞ -∞ dλ e iλs det(I + 2iλΣ) -1/2 V φ ((Σ -1 + 2iλI/D) -1 s -1 ).

If G is the mean-centering projection matrix and Σ G = GΣG, then by the same reasoning as in Thm E.5,

EQUATION

Note that assumption (1) is satisfied if φ is continuous; assumption (4) is satisfied if D ≥ 3 and for all Π, V φ (Π) ≤ Π α for some α ≥ 0; this latter condition, as well as assumptions ( 2) and

(3), will be satisfied if φ is polynomially bounded. Thus the common coordinatewise nonlinearities ReLU, identity, tanh, etc all satisfy these assumptions.

Warning: in general, we cannot swap the order of integration as

∞ -∞ dλ ∞ -∞ ds. For example, if φ = id, then (2π) -D/2 det Σ -1/2 dz e -1 2 z T Σ -1 z φ(z/ z ) ⊗2 = 1 2π lim a 0,b ∞ b a ds ∞ -∞ dλ e iλs det(I + 2iλΣ) -1/2 (Σ -1 + 2iλI) -1 s -1 = 1 2π ∞ -∞ dλ ∞ 0 ds e iλs det(I + 2iλΣ) -1/2 (Σ -1 + 2iλI) -1 s -1

because the s-integral in the latter diverges (in a neighborhood of 0).

Proof. We will only prove the first equation; the others follow similarly.

By dominated convergence,

dze -1 2 z T Σ -1 z φ(z/ z ) ⊗2 = lim a 0,b ∞ dze -1 2 z T Σ -1 z φ(z/ z ) ⊗2 I( z 2 ∈ [a, b]).

Let η : R → R be a nonnegative bump function (i.e. compactly supported and smooth) with support [-1, 1] and integral 1 such that η (0) = η (0) = 0. Then its Fourier transform η(t) decays like O(t -2 ). Furthermore, η (x) := -1 η(x/ ) is a mollifier, i.e. for all f ∈ L 1 (R), f * η → f in L 1 and pointwise almost everywhere. Now, we will show that

dze -1 2 z T Σ -1 z φ(z a /r)φ(z b /r)I( z 2 ∈ [a, b]) = lim 0 dze -1 2 z T Σ -1 z b a d(r 2 )φ(z a /r)φ(z b /r)η ( z 2 -r 2 ) by dominated convergence. Pointwise convergence is immediate, because b a d(r 2 )φ(z a /r)φ(z b /r)η (r 2 -z 2 ) = ∞ -∞ d(r 2 )φ(z a /r)φ(z b /r)I(r 2 ∈ [a, b])η ( z 2 -r 2 ) = [(r 2 → φ(z a /r)φ(z b /r)I(r 2 ∈ [a, b])) * η ]( z 2 ) → φ(z a / z )φ(z b / z )I( z 2 ∈ [a, b])

as → 0 (where we used assumption (1) that φ(z a /r)φ(

z b /r)I(r 2 ∈ [a, b]) is L 1 ).

Finally, we construct a dominating integrable function. Observe

dze -1 2 z T Σ -1 z b a d(r 2 )|φ(z a /r)φ(z b /r)|η ( z 2 -r 2 ) = dze -1 2 z T Σ -1 z b a d(r 2 )|φ(z a /r)φ(z b /r)|η ( z 2 -r 2 ) = dze -1 2 z T Σ -1 z b a d( z 2 + γ)|φ(z a / z 2 + γ)φ(z b / z 2 + γ)|η (-γ) = dze -1 2 z T Σ -1 z b-z 2 a-z 2 dγ |φ(z a / z 2 + γ)φ(z b / z 2 + γ)|η (-γ) = z 2 ∈[a-,b+ ] dze -1 2 z T Σ -1 z b-z 2 a-z 2 dγ |φ(z a / z 2 + γ)φ(z b / z 2 + γ)|η (-γ) since supp η = [-, ] ≤ z 2 ∈[a-,b+ ] dze -1 2 z T Σ -1 z - dγ |φ(z a / z 2 + γ)φ(z b / z 2 + γ)|η (-γ) = - dγ η (-γ) z 2 ∈[a-,b+ ] dze -1 2 z T Σ -1 z |φ(z a / z 2 + γ)φ(z b / z 2 + γ)|

For small enough then, this is integrable by assumption (2), and yields a dominating integrable function for our application of dominated convergence.

In summary, we have just proven that

dze -1 2 z T Σ -1 z φ(z/ z ) ⊗2 = lim a 0,b ∞ lim 0 dze -1 2 z T Σ -1 z b a d(r 2 )φ(z a /r)φ(z b /r)η ( z 2 -r 2 ) Now, dze -1 2 z T Σ -1 z b a d(r 2 )φ(z a /r)φ(z b /r)η ( z 2 -r 2 ) = dze -1 2 z T Σ -1 z b a d(r 2 )φ(z a /r)φ(z b /r) 1 2π ∞ -∞ dλ η (λ)e -iλ( z 2 -r 2 ) = 1 2π dze -1 2 z T Σ -1 z b a d(r 2 )φ(z a /r)φ(z b /r) ∞ -∞ dλ η( λ)e -iλ( z 2 -r 2 )

Note that the absolute value of the integral is bounded above by

1 2π dze -1 2 z T Σ -1 z b a d(r 2 )|φ(z a /r)φ(z b /r)| ∞ -∞ dλ |η( λ)e -iλ( z 2 -r 2 ) | = 1 2π dze -1 2 z T Σ -1 z b a d(r 2 )|φ(z a /r)φ(z b /r)| ∞ -∞ dλ |η( λ)| ≤ C dze -1 2 z T Σ -1 z b a d(r 2 )|φ(z a /r)φ(z b /r)|

for some C, by our construction of η that η(t) = O(t -2 ) for large |t|. By assumption (3), this integral exists. Therefore we can apply the Fubini-Tonelli theorem and swap order of integration.

1 2π dze -1 2 z T Σ -1 z b a d(r 2 )φ(z a /r)φ(z b /r) ∞ -∞ dλ η( λ)e -iλ( z 2 -r 2 ) = 1 2π b a d(r 2 ) ∞ -∞ dλ η( λ)e -iλ( z 2 -r 2 ) dzφ(z a /r)φ(z b /r)e -1 2 z T Σ -1 z = 1 2π b a d(r 2 ) ∞ -∞ dλ η( λ)e iλr 2 dzφ(z a /r)φ(z b /r)e -1 2 z T (Σ -1 +2iλI)z = 1 2π b a d(r 2 )r D ∞ -∞ dλ η( λ)e iλr 2 dzφ(z a )φ(z b )e -1 2 r 2 z T (Σ -1 +2iλI)z = (2π) D/2-1 det Σ 1/2 b a d(r 2 ) ∞ -∞ dλ η( λ)e iλr 2 det(I + 2iλΣ) -1/2 V φ ((Σ -1 + 2iλI) -1 r -2 ) By assumption (4), b a d(r 2 ) ∞ -∞ dλ η( λ)e iλr 2 det(I + 2iλΣ) -1/2 V φ ((Σ -1 + 2iλI) -1 r -2 ) ≤ C b a d(r 2 ) ∞ -∞ dλ det(I + 2iλΣ) -1/2 V φ ((Σ -1 + 2iλI) -1 r -2 )

is integrable (where we used the fact that η is bounded -which is true for all η ∈ C ∞ 0 (R)), so by dominated convergence (applied to 0),

dze -1 2 z T Σ -1 z φ(z/ z ) ⊗2 = lim a 0,b ∞ lim 0 dze -1 2 z T Σ -1 z b a d(r 2 )φ(z a /r)φ(z b /r)η ( z 2 -r 2 ) = (2π) D/2-1 det Σ 1/2 lim a 0,b ∞ b a d(r 2 ) ∞ -∞ dλ e iλr 2 det(I + 2iλΣ) -1/2 V φ ((Σ -1 + 2iλI) -1 r -2 )

where we used the fact that η(0) = R η(x) dx = 1 by construction. This gives us the desired result after putting back in some constants and changing r 2 → s.

E.3 SPHERICAL INTEGRATION:

As mentioned above, the scale invariance of batchnorm combined with the BSB1 fixed point assumption will naturally lead us to consider spherical integration. Here we review several basic facts. Definition E.26. We define the spherical angular coordinates in R B , (r, θ 1 , . .

. , θ B-2 ) ∈ [0, ∞) × A =: [0, π] B-3 × [0, 2π], by x 1 = r cos θ 1 x 2 = r sin θ 1 cos θ 2 x 3 = r sin θ 1 sin θ 2 cos θ 3 . . . x B-2 = r sin θ 1 • • • sin θ B-3 cos θ B-2 x B-1 = r sin θ 1 • • • sin θ B-3 sin θ B-2 . Then we set v :=         cos θ 1 sin θ 1 cos θ 2 sin θ 1 sin θ 2 cos θ 3 . . . sin θ 1 • • • sin θ B-3 cos θ B-2 sin θ 1 • • • sin θ B-3 sin θ B-2        

as the unit norm version of x. The integration element in this coordinate satisfies

dx 1 . . . dx B-1 = r B-2 sin B-3 θ 1 sin B-4 θ 2 • • • sin θ B-3 dr dθ 1 • • • dθ B-2 .

Fact E.27. If F : S B-2 → R, then its expectation on the sphere can be written as

E v∼S B-2 F (v) = 1 2 π 1-B 2 Γ B -1 2 A F (v) sin B-3 θ 1 • • • sin θ B-3 dθ 1 • • • dθ B-2

where A = [0, π] B-3 × [0, 2π] and v are as in Defn E.26.

The following integrals will be helpful for simplifying many expressions involving spherical integration using angular coordinates. Lemma E.28. For j, k ≥ 0 and 0 ≤ s ≤ t ≤ π/2, t s sin j θ cos k θ dθ = 1 2 Beta cos 2 t, cos 2 s; k + 1 2 , j + 1 2 .

By antisymmetry of cos with respect to:

θ → π -θ, if π/2 ≤ s ≤ t ≤ π, t s sin j θ cos k θ dθ = 1 2 (-1) k Beta cos 2 s, cos 2 t; k + 1 2 , j + 1 2 . Proof. Set x := cos 2 θ =⇒ dx = -2 cos θ sin θ dθ. So the integral in question is -1 2 cos 2 t cos 2 s (1 -x) j-1 2 x k-1 2 dx = 1 2 (Beta(cos 2 s; k + 1 2 , j + 1 2 ) -Beta(cos 2 t; k + 1 2 , j + 1 2 ))

As consequences, Lemma E.29. For j, k ≥ 0,

π 0 sin j θ cos k θ dθ = 1 + (-1) k 2 Beta( j + 1 2 , k + 1 2 ) = 1 + (-1) k 2 Γ( 1+j 2 )Γ( 1+k 2 ) Γ( 2+j+k 2 )

.

Lemma E.30.

π 0 sin θ 1 dθ 1 π 0 sin 2 θ 2 dθ 2 • • • π 0 sin k θ k dθ k = π k/2 /Γ( k + 2 2 ).

Proof. By Lemma E.29,

EQUATION

2 ) . Thus this product of integrals is equal to

k r=1 √ π Γ( r+1 2 ) Γ( r+2 2 ) = π k/2 Γ(1)/Γ( k + 2 2 ) = π k/2 /Γ( k + 2 2 ). Lemma E.31. For B > 1, µ > 0, ∞ 0 r B e -r 2 /2µ dr = 2 B-1 2 µ B+1 2 Γ( B+1 2 )

Proof. Apply change of coordinates z = r 2 /2µ.

Definition E.32. Let K(v; Σ, β) := (det e T Σe) -1/2 (v T (e T Σe) -1 v) -B-1+β 2 and K(v; Σ) := K(v; Σ, 0).

The following proposition allows us to express integrals involving batchnorm as expectations over the uniform spherical measure.

Proposition E.33. If f : R B-1 → R B , B > 1 -β, and Σ ∈ R (B-1)×(B-1) is nonsingular, then E h∼N (0,Σ) h β f (h/ h ) = 2 β/2 P B -1 2 , β 2 E v∼S B-2 f (v)(det Σ) -1/2 (v T Σ -1 v) -B-1+β 2 .

If f : R B → R B and Σ ∈ R B×B is such that e T Σe is nonsingular (where e is as in Defn E.3) and

GΣG = Σ, E h∼N (0,Σ) Gh β f (Gh/ Gh ) = 2 β/2 P B -1 2 , β 2 E v∼S B-2 f (ev)K(v; Σ, β).

In particular,

E[f (B(h)) : h ∼ N (0, Σ)] = E v∼S B-2 f (ev)K(v; Σ)

Proof. We will prove the second statement. The others follow trivially.

E[ Gh β f (B(h)) : h ∼ N (0, Σ)] = E[ x β f (ex/ x ) : x ∼ N (0, e T Σe)] = (2π) 1-B 2 (det e T Σe) -1/2 R B-1 x β f (ex/ x )e -1 2 x T (e T Σe) -1 x dx = (2π) 1-B 2 (det e T Σe) -1/2 A f (ev) sin B-3 θ 1 • • • sin θ B-3 dr dθ 1 • • • dθ B-2 × ∞ 0 r B-2+β e -1 2 r 2 •v T (e T Σe) -1 v (spherical coordinates Defn E.26) = (2π) 1-B 2 (det e T Σe) -1/2 × 2 B-3 2 Γ B -1 2 × 2 B-3+β 2 Γ( B-1+β 2 ) 2 B-3 2 Γ B-1 2 × A f (ev) sin B-3 θ 1 • • • sin θ B-3 dθ 1 • • • dθ B-2 (v T (e T Σe) -1 v) B-1+β 2 = 1 2 π 1-B 2 Γ B -1 2 × 2 β/2 Γ B -1 + β 2 Γ B -1 2 -1 × A f (ev) sin B-3 θ 1 • • • sin θ B-3 dθ 1 • • • dθ B-2 (v T (e T Σe) -1 v) B-1+β 2 = 2 β/2 Γ B -1 + β 2 Γ B -1 2 -1 E v∼S B-2 f (ev) (det e T Σe) 1/2 (v T (e T Σe) -1 v) B-1+β 2 = 2 β/2 Γ B -1 + β 2 Γ B -1 2 -1 E v∼S B-2 f (ev)K(v; Σ, β)

As a straightforward consequence, Proposition E.34. If e T Σe is nonsingular,

V B φ (Σ) = E v∼S B-2 φ( √ Bev) ⊗2 K(v; Σ)

For numerical calculations, it is useful to realize e as the matrix whose columns are e B-m := (m(m + 1)) -1/2 (0, . . . , 0, -m, m 1, 1, . . . , 1) T , for each m = 1, . . . , B -1.

e T =         1 √ B(B-1) 1 √ (B-1)(B-2)

. . .

EQUATION

We have

(ev) 1 = - B -1 B cos θ 1 (ev) 2 = cos θ 1 - B -2 B -1 sin θ 1 cos θ 2 (ev) 3 = cos θ 1 + sin θ 1 cos θ 2 - B -3 B -2 sin θ 1 sin θ 2 cos θ 3 (ev) 4 = cos θ 1 + sin θ 1 cos θ 2 + sin θ 1 sin θ 2 cos θ 3 - B -4 B -3 sin θ 1 sin θ 2 sin θ 3 cos θ 4 (34)

so in particular they only depend on θ 1 , . . . , θ 4 .

Angular coordinates In many situations, the sparsity of this realization of e combined with sufficient symmetry allows us to simplify the high-dimensional spherical integral into only a few dimensions, using one of the following lemmas Lemma E.35. Let f : R 4 → R A for some A ∈ N. Suppose B ≥ 6. Then for any k < B -1

E[r -k f (v 1 , . . . , v 4 ) : x ∼ N (0, I B-1 )] = Γ((B -1 -k)/2) Γ((B -5)/2) 2 -k/2 π -2 π 0 dθ 1 • • • π 0 dθ 4 f (v θ 1 , . . . , v θ 4 ) sin B-3 θ 1 • • • sin B-6 θ 4

where v i = x i / x and r = x , and

v θ 1 = cos θ 1 v θ 2 = sin θ 1 cos θ 2 v θ 3 = sin θ 1 sin θ 2 cos θ 3 v θ 4 = sin θ 1 sin θ 2 sin θ 3 cos θ 4 .

Lemma E.36. Let f : R 2 → R A for some A ∈ N. Suppose B ≥ 4. Then for any k < B -1

E[r -k f (v 1 , v 2 ) : x ∼ N (0, I B-1 )] = Γ((B -1 -k)/2) Γ((B -3)/2) 2 -k/2 π -1 π 0 dθ 1 π 0 dθ 2 f (v θ 1 , v θ 2 ) sin B-3 θ 1 sin B-4 θ 2

where v i = x i / x and r = x , and

v θ 1 = cos θ 1 v θ 2 = sin θ 1 cos θ 2 Lemma E.37. Let f : R → R A for some A ∈ N. Suppose B ≥ 3. Then for any k < B -1 E[r -k f (v) : x ∼ N (0, I B-1 )] = Γ((B -1 -k)/2) Γ((B -2)/2) 2 -k/2 π -1/2 π 0 dθ f (cos θ) sin B-3 θ

where v i = x i / x and r = x .

We will prove Lemma E.35; those of Lemma E.36 and Lemma E.37 are similar.

Proof of Lemma E.35.

E[r -k f (v 1 , . . . , v 4 ) : x ∼ N (0, I B-1 )] = (2π) -(B-1)/2 dx r -k f (v 1 , . . . , v 4 )e -r 2 /2 = (2π) -(B-1)/2 r B-2 sin B-3 θ 1 • • • sin θ B-3 dr dθ 1 • • • dθ B-2 r -k f (v θ 1 , . . . , v θ 4 )e -r 2 /2 = (2π) -(B-1)/2 f (v θ 1 , . . . , v θ 4 ) sin B-3 θ 1 • • • sin θ B-3 dθ 1 • • • dθ B-2 ∞ 0 r B-2-k e -r 2 /2 dr = (2π) -(B-1)/2 f (v θ 1 , . . . , v θ 4 ) sin B-3 θ 1 • • • sin θ B-3 dθ 1 • • • dθ B-2 2 (B-3-k)/2 Γ((B -1 -k)/2) (change of coordinate r ← r 2 and definition of Γ function) = Γ((B -1 -k)/2)2 -(2+k)/2 π -(B-1)/2 f (v θ 1 , . . . , v θ 4 ) sin B-3 θ 1 • • • sin θ B-3 dθ 1 • • • dθ B-2

Because v θ 1 , . . . , v θ 4 only depends on θ 1 , . . . , θ 4 , we can integrate out θ 5 , . . . , θ B-2 . By applying Lemma E.30 to θ 5 , . . . θ B-3 and noting that

2π 0 dθ B-2 = 2π, we get E[r -k f (v θ 1 , . . . , v θ 4 ) : x ∼ N (0, I B-1 )] = Γ((B -1 -k)/2)2 -(2+k)/2 π -(B-1)/2 f (v θ 1 , . . . , v θ 4 ) sin B-3 θ 1 • • • sin B-6 θ 4 π (B-7)/2 Γ((B -5)/2) -1 (2π) = Γ((B -1 -k)/2) Γ((B -5)/2) 2 -k/2 π -2 f (v θ 1 , . . . , v θ 4 ) sin B-3 θ 1 • • • sin B-6 θ 4

Cartesian coordinates. We can often also simplify the high dimensional spherical integrals without trigonometry, as the next two lemmas show. Both are proved easily using change of coordinates. Lemma E.38. For any u ∈ S n ,

E v∼S n f (u • v) = ω n-1 ω n 1 -1 dh (1 -h 2 ) n-2 2 f (h)

where ω n-1 is hypersurface of the (n -1)-dimensional unit hypersphere. Lemma E.39. For any u 1 , u 2 ∈ S n ,

E v∼S n f (u 1 • v, u 2 • v) = ω n-2 ω n v ≤1 v∈span(u1,u2) dv (1 -v 2 ) n-3 2 f (u 1 • v, u 2 • v)

where ω n-1 is hypersurface of the (n -1)-dimensional unit hypersphere.

E.4 GEGENBAUER EXPANSION:

Definition E.40. The Gegenbauer polynomials {C

(α) l (x)} ∞ l=0 , with deg C (α)

l (x) = l, are the set of orthogonal polynomials with respect to the weight function

(1 -x 2 ) α-1/2 (i.e. 1 -1 (1 - x 2 ) α-1/2 C (α) l (x)C (α) m (x) = 0 if l = m). By convention, they are normalized so that 1 -1 C (α) n (x) 2 (1 -x 2 ) α-1 2 dx = π2 1-2α Γ(n + 2α) n!(n + α)Γ(α) 2 .

Here are several examples of low degree Gegenbauer polynomials

C (α) 0 (x) = 1 C (α) 1 (x) = 2αx C (α) 2 (x) = -α + 2α(1 + α)x 2 C (α) 3 (x) = -2α(1 + α)x + 4 3 α(1 + α)(2 + α)x 3 .

They satisfy a few identities summarized below

Fact E.41 (Suetin). C (α) n (±1) = (±1) n n+2α-1

n Fact E.42 (Weisstein) .

C (λ) n (x) = 2λC (λ+1) n-1 (x) (35) xC (λ) n (x) = nC (λ) n (x) + C (λ) n-1 (x) (36) C (λ) n+1 (x) = (n + 2λ)C (λ) n (x) + xC (λ) n (x) (37) nC (λ) n (x) = 2(n -1 + λ)xC (λ) n-1 (x) -(n -2 + 2λ)C (λ) n-2 (x), ∀n ≥ 2 (38)

By repeated applications of Eq. ( 38), Proposition E.43.

xC (λ) n (x) = 1 2(n + λ) (n + 1)C (λ) n+1 (x) + (n -1 + 2λ)C (λ) n-1 (x) , ∀n ≥ 1 x 2 C (λ) n (x) = κ 2 (n, λ)C (λ) n+2 (x) + κ 0 (n, λ)C (λ) n (x) + κ -2 (n, λ) C (λ) n-2 (x), ∀n ≥ 2 xC (λ) 0 (x) = x = 1 2λ C (λ) 1 (x) x 2 C (λ) 0 (x) = 1 4λ(1 + λ) 2C (λ) 2 (x) + 2λC (λ) 0 (x) = κ 2 (0, λ)C (λ) 2 (x) + κ 0 (0, λ)C (λ) 0 (x) x 2 C (λ) 1 (x) = κ 2 (1, λ)C (λ) 3 (x) + κ 0 (1, λ)C (λ) 1 (x) where κ 2 (n, λ) := (n+1)(n+2) 4(n+λ)(n+1+λ) , κ 0 (n, λ) := (n+1)(n+2λ) 4(n+λ)(n+1+λ) + (n-1+2λ)n 4(n+λ)(n-1+λ) , κ -2 (n, λ) := (n-1+2λ)(n-2+2λ) 4(n+λ)(n-1+λ

) . This proposition is useful for the Gegenbauer expansion of the local convergence rate of Eq. ( 43).

Zonal harmonics. The Gegenbauer polynomials are closely related to a special type of spherical harmonics called zonal harmonics, which are intuitively those harmonics that only depend on the "height" of the point along some fixed axis.

Theorem E.47 (Dirichlet-like form of φ diagonalizes over Gegenbauer basis:

). Let u 1 , u 2 ∈ S B-2 . Suppose φ( √ B -1x), φ ( √ B -1x) ∈ L 2 ((1 -x 2 ) B-3 2 ) and φ( √ B -1x) has Gegenbauer ex- pansion ∞ l=0 a l 1 c B-1,l C ( B-3 2 ) l (x). Then E v∈S B-2 (u 1 • u 2 -(u 1 • v)(u 2 • v))φ ( √ B -1u 1 • v)φ ( √ B -1u 2 • v) = ∞ l=0 a 2 l (l + B -3)lc -1 B-1,l C ( B-3 2 ) l (u 1 • u 2 )

This is proved using the following lemmas.

Lemma E.48.

xC (α) l (x) = lC (α) l (x) + j∈[l-2]2 2(j + α)C (α) j (x) C (α) l (x) = j∈[l-1]2 2(j + α)C (α) j (x)

where [n] 2 is the sequence {(n%2), (n%2) + 1, . . . , n -2, n} and n%2 is the remainder of n/2.

Proof. Apply Eq. ( 36) and Eq. ( 37) alternatingly.

Lemma E.49.

lC (α) l (t) + j∈[l-2]2 2(j + α)C (α) j (t) -t j∈[l-1]2 2(j + α)C (α) j (t) = 0

Proof. Apply Eq. ( 38) to lC 

φ ( √ B -1x) = 1 √ B -1 ∞ l=0 a l 1 c B-1,l C ( B-3 2 ) l (x) = 1 √ B -1 ∞ l=0 a l 1 c B-1,l j∈[l-1]2 (2j + B -3)C ( B-3 2 ) j (x) = 1 √ B -1 ∞ l=0 j∈[l-1]2 a l c B-1,j c B-1,l (2j + B -3)c -1 B-1,j C ( B-3 2 ) j (x) = 1 √ B -1 ∞ j=0 c -1 B-1,j C ( B-3 2 ) j (x) l∈[j+1,∞]2 a l c B-1,j c B-1,l (2j + B -3) where [n, m] 2 is the sequence n, n + 2, . . . , m. So for any u 1 , u 2 ∈ S B-2 , E v∈S B-2 φ ( √ B -1u 1 • v)φ ( √ B -1u 2 • v) = 1 B -1 ∞ j=0 c -1 B-1,j C ( B-3 2 ) j (u 1 • u 2 )   l∈[j+1,∞]2 a l c B-1,j c B-1,l (2j + B -3)   2 = 1 B -1 ∞ j=0 c B-1,j (2j + B -3) 2 C ( B-3 2 ) j (u 1 • u 2 )   l,m∈[j+1,∞]2 a l c B-1,l a m c B-1,m   .

For l < m and ml even, the coefficient of a l a m is

EQUATION

)

For each l, the coefficient of

EQUATION

EQUATION

2 ), we have

xφ ( √ B -1x) = 1 √ B -1 ∞ l=0 a l 1 c B-1,l xC ( B-3 2 ) l (x) = 1 √ B -1 ∞ l=0 a l 1 c B-1,l   lC ( B-3 2 ) l (x) + j∈[l-2]2 (2j + B -3)C ( B-3 2 ) j (x)   = 1 √ B -1 ∞ j=0 c -1 B-1,j C ( B-3 2 ) j (x)   a j j + l∈[j+2,∞]2 a l c B-1,j c B-1,l (2j + B -3)   Then E v∈S B-2 (u 1 • v)(u 2 • v)φ ( √ B -1u 1 • v)φ ( √ B -1u 2 • v) = 1 B -1 ∞ j=0 c -1 B-1,j C ( B-3 2 ) j (u 1 • u 2 )   a j j + l∈[j+2,∞]2 a l c B-1,j c B-1,l (2j + B -3)   2

For l < m and ml even, the coefficient of a l a m is

1 B -1 2l(2l + B -3) c B-1,l c B-1,m c -1 B-1,l C ( B-3 2 ) l (u 1 • u 2 ) +2 j∈[l-2]2 c B-1,j c B-1,l c B-1,j c B-1,m (2j + B -3) 2 c -1 B-1,j C ( B-3 2 ) j (u 1 • u 2 )   = 2c -1 B-1,l c -1 B-1,m B -3 B -1   lC ( B-3 2 ) l (u 1 • u 2 ) + j∈[l-2]2 (2j + B -3)C ( B-3 2 ) j (u 1 • u 2 )   (41) For each l, the coefficient of a 2 l is 1 B -1   l 2 c -1 B-1,l C ( B-3 2 ) l (u 1 • u 2 ) + j∈[l-2]2 c B-1,j c B-1,l 2 (2j + B -3) 2 c -1 B-1,j C ( B-3 2 ) j (u 1 • u 2 )   = c -2 B-1,l B -3 B -1   l 2 2l + B -3 C ( B-3 2 ) l (u 1 • u 2 ) + j∈[l-2]2 (2j + B -3)C ( B-3 2 ) j (u 1 • u 2 )   (42)

By Lemma E.49 and Eqs. ( 39) and ( 41), the coefficient of a l a m with m > l in

E v∈S B-2 (u 1 •u 2 -(u 1 • v)(u 2 •v))φ ( √ B -1u 1 •v)φ ( √ B -1u 2 •v) is 0.

Similarly, by Lemma E.49 and Eqs. ( 40) and ( 42), the coefficient of

a 2 l is c -2 B-1,l B-3 B-1 l+B-3 2l+B-3 lC ( B-3 2 ) l (u 1 • u 2 ) = c -1 B-1,l l+B-3 B-1 lC ( B-3 2 ) l (u 1 • u 2 ).

E.5 SYMMETRY AND ULTRASYMMETRY:

As remarked before, all commonly used nonlinearities ReLU, tanh, and so on induce the dynamics Eq. ( 43) to converge to BSB1 fixed points, which we formally define as follows. Definition E.50. We say a matrix Σ ∈ S B is BSB1 (short for "1-Block Symmetry Breaking") if Σ has one common entry on the diagonal and one common entry on the off-diagonal, i.e.

    a b b • • • b a b • • • b b a • • • . . . . . . . . . . . .    

We will denote such a matrix as BSB1(a, b). Note that BSB1(a, b) can be written as (a-b)I +b11 T . Its spectrum is given below. Lemma E.51. A B × B matrix of the form µI + ν11 T has two eigenvalues µ and Bν + µ, each with eigenspaces {x : i x i = 0} and {x :

x 1 = • • • = x B }.

Equivalently, if it has a on the diagonal and b on the off-diagonal, then the eigenvalues are ab and (B -1)b + a.

The following simple lemma will also be very useful Lemma E.52.

Let Σ := BSB1 B (a, b). Then GΣG = (a -b)G.

Proof. G and BSB1 B (a, b) can be simultaneously diagonalized by Lemma E.51. Note that G zeros out the eigenspace R1, and is identity on its orthogonal complement. The result then follows from easy computations.

The symmetry of BSB1 covariance matrices, especially the fact that they are isotropic in a codimension 1 subspace, is crucial to much of our analysis.

Ultrasymmetry. Indeed, as mentioned before, when investigating the asymptotics of the forward or backward dynamics, we encounter 4-tensor objects, such as

dV B φ

dΣ or V B φ , that would seem daunting to analyze at first glance. However, under the BSB1 fixed point assumption, these 4tensors contain considerable symmetry that somehow makes this possible. We formalize the notion of symmetry arising from such scenarios: Definition E.53. Let T : R B×B → R B×B be a linear operator. Let δ i ∈ R B be the vector with 0 everywhere except 1 in coordinate i; then δ i δ T j is the matrix with 0 everywhere except 1 in position (i, j). Write [kl|ij] := T (δ i δ T j ) kl . Suppose T has the property that for all i, j, k, l ∈ [B] = {1, . . . , B}

• [kl|ij] = [π(k)π(l)|π(i)π(j)

] for all permutation π on [B], and

• [ij|kl] = [ji|lk].

Then we say T is ultrasymmetric. Remark E.54. In what follows, we will often "normalize" the representation "[ij|kl]" to the unique "[i j |k l ]" that is in the same equivalence class according to Defn E.53 and such that i , j , k , l ∈ [4] and i ≤ j , k ≤ l unless i = l , j = k , in which case the normalization is [12|21]. Explicitly, we have the following equivalence classes and their normalized representations We will study the eigendecomposition of an ultrasymmetric T as well as the projection Theorem E.61 (Eigendecomposition of an ultrasymmetric operator). Let T : R B×B → R B×B be an ultrasymmetric linear operator. Then T has the following eigendecomposition.

G ⊗2 • T • G ⊗2 : Σ → G(T {GΣG})G,

1. Two 1-dimensional eigenspace R • BSB1(λ T BSB1,iα 22 , α 21 ) with eigenvalue λ T BSB1,i for i = 1, 2, where

α 11 = [11|11] + [11|22](B -1) α 12 = 2(B -1)[11|12] + (B -2)(B -1)[11|23] α 21 = 2[12|11] + (B -2)[12|33] α 22 = [12|12] + 4(B -2)[12|13] + [12|21] + (B -2)(B -3)[12|34]

and λ T BSB1,1 and λ T BSB1,2 are the roots to the quadratic

x 2 -(α 11 + α 22 )x + α 11 α 22 -α 12 α 21 . 2. Two (B -1)-dimensional eigenspaces S B • L(λ T L,i -β 22 , β 21 ) with eigenvalue λ T L,i for i = 1, 2.

Here S B • W denotes the linear span of the orbit of matrix W under simultaneous permutation of its column and rows (by the same permutation), and 

β 11 = [11|11] -[11|22] β 12 = 2(B -2)([11|23] -[11|12]) β 21 = -[12|11] + [12|33] β 22 = [12|21] + [12|12] + 2(B -4)[12|13] -2(B -3)[12|34].

3.. Eigenspace M (dimension:

B(B -3)/2) with eigenvalue λ T M := [12|12] + [12|21] - 4[12|13] + 2[12|34].

The proof is by careful, but ultimately straightforward, computation.

Proof. We will use the bracket notation of Defn E.53 to denote entries of T , and implicitly simplify it according to Remark E.54. Item 1. Let U ∈ R B×B be the BSB1 matrix. By ultrasymmetry of T and BSB1 symmetry of A, T {U } is also BSB1. So we proceed to calculate the diagonal and off-diagonal entries of T {G}.

We have:

T {BSB1(a, b)} 11 = [11|11]a + 2(B -1)[11|12]b + [11|22](B -1)a + [11|23](B -2)(B -1)b = ([11|11] + [11|22](B -1))a + (2(B -1)[11|12] + (B -2)(B -1)[11|23])b T {BSB1(a, b)} 12 = [12|12]b + 2[12|11]a + (B -2)[12|33]a + 2(B -2)[12|13]b + [12|21]b + 2(B -2)[12|23]b + (B -2)(B -3)[12|34]b = (2[12|11] + (B -2)[12|33])a + ([12|12] + 2(B -2)[12|13] + [12|21] + 2(B -2)[12|13] + (B -2)(B -3)[12|34])b = (2[12|11] + (B -2)[12|33])a + ([12|12] + 4(B -2)[12|13] + [12|21] + (B -2)(B -3)[12|34])b

Thus BSB1(ω 1 , γ 1 ) and BSB1(ω 2 , γ 2 ) are the eigenmatrices of T , where (ω 1 , γ 1 ) and (ω 2 , γ 2 ) are the eigenvectors of the matrix α 11 α 12 α 21 α 22 with

α 11 = [11|11] + [11|22](B -1) α 12 = 2(B -1)[11|12] + (B -2)(B -1)[11|23] α 21 = 2[12|11] + (B -2)[12|33] α 22 = [12|12] + 4(B -2)[12|13] + [12|21] + (B -2)(B -3)[12|34]

The eigenvalues are the two roots λ T BSB1,1 , λ T BSB1,2 to the quadratic x 2 -(α 11 + α 22 )x + α 11 α 22α 12 α 21 and the corresponding eigenvectors are

(ω 1 , γ 1 ) = (λ 1 -α 22 , α 21 ) (ω 2 , γ 2 ) = (λ 2 -α 22 , α 21 ).

Item 2. We will study the image of L B (a, b) (Defn E.57) under T . We have

T {L(a, b)} 11 = -T {L(a, b)} 22 = [11|11]a + 2(B -2)[11|12](-b) + [11|22](-a) + 2(B -2)[11|23]b = ([11|11] -[11|22])a + 2(B -2)([11|23] -[11|12])b T {L(a, b)} 12 = T {L(a, b)} 21 = [12|12]0 + [12|21]0 + [12|11](a -a) + [12|13](b -b) + [12|31](b -b) + [12|33]0 + [12|34]0 = 0 T {L(a, b)} 33 = T {L(a, b)} ii , ∀i ≥ 3 = [11|11]0 + [11|12](2b -2b) + [11|22](a -a) + [11|23](2(B -3)b -2(B -3)b) = 0 T {L(a, b)} 34 = T {L(a, b)} ij , ∀i = j & i, j ≥ 3 = [12|12]0 + [12|11]0 + [12|13](2b -2b) + [12|21]0 + [12|31](2b -2b) + [12|33](a -a) + [12|34](2(B -4)b -2(B -4)b) = 0 T {L(a, b)} 13 = T {L(a, b)} 1j , ∀j ≥ 3 = T {L(a, b)} j1 , ∀j ≥ 3 = [12|12](-b) + [12|13](B -3 -1)(-b) + [12|11]a + [12|21](-b) + [12|31](B -3 -1)(-b) + [12|33](-a) + [12|34]2(B -3)b = [12|11] -[12|33] a + 2(B -3)[12|34] -2(B -4)[12|13] -[12|12] -[12|21] b

Thus L(a, b) transforms under T by the matrix

L(a, b) → L( β 11 β 12 β 21 β 22 a b )

with

β 11 = [11|11] -[11|22] β 12 = 2(B -2)([11|23] -[11|12]) β 21 = -[12|11] + [12|33] β 22 = [12|21] + [12|12] + 2(B -4)[12|13] -2(B -3)[12|34].

So if λ T L,1 and λ T L,2 are the roots of the equation

x 2 -(β 11 + β 22 )x + β 11 β 22 -β 12 β 21 then T {L(λ T L,1 -β 22 , β 21 )} = λ T L,1 L(λ T L,1 -β 22 , β 21 ) T {L(λ T L,2 -β 22 , β 21 )} = λ T L,2 L(λ T L,2 -β 22 , β 21 )

Proof. Item 1 As in the proof of Thm E.61, we find

T {BSB1(a, b)} = BSB1 α 11 α 12 α 21 α 22 a b

where

α 11 = [11|11] + [11|22](B -1) α 12 = 2(B -1)[11|12] + (B -2)(B -1)[11|23] α 21 = 2[12|11] + (B -2)[12|33] α 22 = [12|12] + 4(B -2)[12|13] + [12|21] + (B -2)(B -3)[12|34]. For a = B -1, b = -1 so that BSB1(B -1, -1) = BG, we get T {BSB1(B -1, -1)} = BSB1 ((B -1)α 11 -α 12 , (B -1)α 21 -α 22 ) G ⊗2 • T {BSB1(B -1, -1)} = G BSB1 ((B -1)α 11 -α 12 , (B -1)α 21 -α 22 ) G = ((B -1)(α 11 -α 21 ) -(α 12 -α 22 )) G = B -1 ((B -1)(α 11 -α 21 ) -(α 12 -α 22 )) BSB1(B -1, -1) = λ G,T G BSB1(B -1, -1) by Lemma E.52. Item 2. It suffices to show that L(B -2, 1) is an eigenmatrix with the eigenvalue λ G,T L .

As in the proof of Thm E.61, we find

T {L(a, b)} = L β 11 β 12 β 21 β 22 a b

where

β 11 = [11|11] -[11|22] β 12 = 2(B -2)([11|23] -[11|12]) β 21 = -[12|11] + [12|33] β 22 = [12|21] + [12|12] + 2(B -4)[12|13] -2(B -3)[12|34].

So with a = B -2, b = 1, we have

T {L(B -2, 1)} = L((B -2)β 11 + β 12 , (B -2)β 21 + β 22 ) G ⊗2 • T {L(B -2, 1)} = G L((B -2)β 11 + β 12 , (B -2)β 21 + β 22 ) G = B -1 ((B -2)β 11 + β 12 + 2((B -2)β 21 + β 22 ))L(B -2, 1) = λ G,T L L(B -2, 1)

by Lemma E.60.

Item 3. The proof is exactly the same as in that of Thm E.61.

Noting that the eigenspaces of Thm E.62 are orthogonal, we have Proposition E.63. G ⊗2 • T H G for any ultrasymmetric operator T is self-adjoint.

DOS Operators. In some cases, such as the original study of vanilla tanh networks by Schoenholz et al. ( 2016), the ultrasymmetric operator involved has a much simpler form: Definition E.64. Let T : H B → H B be such that for any Σ ∈ H B and any i = j ∈ [B],

T {Σ} ii = uΣ ii T {Σ} ij = vΣ ii + vΣ jj + wΣ ij

Then we say that T is diagonal-off-diagonal semidirect, or DOS for short. We write more specifically T = DOS B (u, v, w).

Proof. Let λ 1 l ≥ λ 2 l ≥ • • • ≥ λ A l

be the eigenvalues of Σ l . It's easy to see that i λ i l = 1 for all l ≥ 1. So WLOG we assume l ≥ 1 and this equality holds.

We will show the "exponential convergence" statement; that lim l→∞ Σ l = 1 A I then follows from the trace condition above.

Proof of Item 2. For brevity, we write suppress the superscript index l, and use • to denote • l . We will now compute the eigenvalues

λ 1 ≥ • • • ≥ λ A of Σ.

First, notice that Σ and Σ can be simultaneously diagonalized, and by induction all of {Σ l } l≥0 can be simultaneous diagonalized. Thus we will WLOG assume that Σ is diagonal, Σ = Diag(λ 1 , . . . , λ A ), so that Σ = Diag(γ 1 , . . . , γ A ) for some {γ i } i . These {γ i } i form the eigenvalues of Σ but a priori we don't know whether they fall in decreasing order; in fact we will soon see that they do, and

γ i = λ i .

We have

EQUATION

Therefore,

γ i -γ k = ∞ 0 A j=1 (1 + 2sλ j ) -1/2 (λ i (1 + 2sλ i ) -1 -λ k (1 + 2sλ k ) -1 ) ds = (λ i -λ k ) ∞ 0 A j=1 (1 + 2sλ j ) -1/2 (1 + 2sλ i ) -1 (1 + 2sλ k ) -1 ds Since the RHS integral is always positive, λ i ≥ λ k =⇒ γ i ≥ γ k and thus γ i = λ i for each i. Define T (λ) := ∞ 0 A j=1 (1 + 2sλ 1 ) -1/2 (1 + 2sλ 1 ) -1 (1 + 2sλ A ) -1 ds, so that λ 1 -λ A = (λ 1 -λ A )T (λ).

Note first that:

A j=1 (1 + 2sλ j ) -1/2 (1 + 2sλ i ) -1 (1 + 2sλ k ) -1

is (strictly) log-convex and hence (strictly) convex in λ. Furthermore, T (λ) is (strictly) convex because it is an integral of (strictly) convex functions. Thus T is maximized over any convex region by its extremal points, and only by its extremal points because of strict convexity.

The convex region we are interested in is given by

A := {(λ i ) i : λ 1 ≥ λ 2 ≥ • • • ≥ λ A ≥ 0 & i λ i = 1}.

The unit sum condition follows from the normalization h/ h of the iteration map. The extremal points of A are {ω k := (

EQUATION

where the last line follows because ω A is the only point with last coordinate nonzero so that a A = λ A .

We now show that the gap λ 1 lλ A l → 0 as l → ∞. There are two cases: if λ A l is bounded away from 0 infinitely often, then T (λ) < 1 -infinitely often for a fixed > 0 so that the gap indeed vanishes with l. Now suppose otherwise, that λ A l converges to 0; we will show this leads to a contradiction. Notice that

λ A = ∞ 0 A j=1 (1 + 2sλ j ) -1/2 λ A (1 + 2sλ A ) -1 ds λ A /λ A ≥ ∞ 0 A j=1 (1 + 2s(1 -λ A )/(A -1)) -(A-1)/2 (1 + 2sλ A ) -3/2 ds

where the first lines is Eq. ( 44) and the 2nd line follows from the convexity of

A-1 j=1 (1 + 2sλ j ) -1/2

as a function of (λ 1 , . . . , λ A-1 ). By a simple application of dominated convergence, as λ A → 0, this integral converges to a particular simple form,

λ A /λ A ≥ ∞ 0 A j=1 (1 + 2s/(A -1)) -(A-1)/2 ds = 1 + 2 A -3

Thus for large enough l, λ A l+1 /λ A l is at least 1 + for some > 0, but this contradicts the convergence of λ A to 0.

Altogether, this proves that λ 1 lλ A l → 0 and therefore λ l → ω A as l → ∞. Consequently, λ A l is bounded from below, say by K (where K > 0 because by Eq. ( 44), λ A l is never 0), for all l, and by Eq. ( 45), we prove Item 2 by taking K = 1 -K 2 A+2 . In addition, asymptotically, the gap decreases exponentially as

T (ω A ) l = A A+2 l , proving Item 3. Theorem F.2. Consider the dynamics Σ l = E[(h/ h ) ⊗2 : h ∼ N (0, Σ l-1 )] on Σ l ∈ S A . Suppose Σ 0 = M T DM where M is orthogonal and D is a diagonal matrix. If D has rank C and D ii = 0, ∀1 ≤ i ≤ C, then 1. lim l→∞ Σ l = 1 C M T D M

where D is the diagonal matrix with D ii = I(D ii = 0). 2. This convergence is exponential in the sense that, for any Σ 0 of rank C, there is a constant

K < 1 such that λ 1 (Σ l ) -λ C (Σ l ) < K(λ 1 (Σ l-1 ) -λ C (Σ l-1

)) for all l ≥ 2. Here λ i denotes the ith largest eigenvalue.

3. Asymptotically, λ 1 (Σ l ) -λ C (Σ l ) = O((1 -2 C+2 ) l ).

Proof. Note that Σ l can always be simultaneously diagonalized with Σ 0 , so that

Σ l = M T D l M for some diagonal D l which has D l ii = 0, ∀0 ≤ i ≤ C. Then we have (D l ) = E[(h/ h ) ⊗2 : h ∼ N (0, (D l-1 ) )],

where (D l ) means the diagonal matrix obtained from D l by deleting the dimensions with zero eigenvalues. The proof then finishes by Lemma F.1.

From this it easily follows the following characterization of the convergence behavior. Corollary F.3. Consider the dynamics of Eq. ( 43)

for φ = id: Σ l = B E[(Gh/ Gh ) ⊗2 : h ∼ N (0, Σ l-1 )] on Σ l ∈ S B .

Suppose GΣG has rank C < B and factors as êDê T where D ∈ R C×C is a diagonal matrix with no zero diagonal entries and ê is an B × C matrix whose columns form an orthonormal basis of a subspace of im G. Then

1. lim l→∞ Σ l = B C êI C êT .

2. This convergence is exponential in the sense that, for any GΣ 0 G of rank C, there is a constant

K < 1 such that λ 1 (Σ l ) -λ C (Σ l ) < K(λ 1 (Σ l-1 ) -λ C (Σ l-1

)) for all l ≥ 2.

Here λ i denotes the ith largest eigenvalue.

3.. Asymptotically:

, λ 1 (Σ l ) -λ C (Σ l ) = O((1 -2 C+2 ) l ).

General nonlinearity. We don't (currently) have a proof of any characterization of the basin of attraction for Eq. ( 43) for general φ. Thus we are forced to resort to finding its fixed points manually and characterize their local convergence properties.

F.2 LIMIT POINTS:

Batchnorm B φ is permutation-equivariant, in the sense that B φ (πh) = πB φ (h) for any permtuation matrix π. Along with the case of φ = id studied above, this suggests that we look into BSB1 fixed points Σ * . What are the BSB1 fixed points of Eq. (43)?

F.2.1 SPHERICAL INTEGRATION:

The main result (Thm F.5) of this section is an expression of the BSB1 fixed point diagonal and offdiagonal entries in terms of 1-and 2-dimensional integrals. This allows one to numerically compute such fixed points.

By a simple symmetry argument, we have the following Lemma F.4. Suppose X is a random vector in R B , symmetric in the sense that for any permutation matrix π and any subset U of R B measurable with respect to the distribution of X, P (X ∈ U ) = P (X ∈ π(U )). Let Φ : R B → R B be a symmetric function in the sense that Φ(πx) = πΦ(x) for any permutation matrix π. Then E Φ(X)Φ(X) T = µI + ν11 T for some µ and ν.

This lemma implies that V B φ (Σ) is BSB1 whenever Σ is BSB1. In the following, We in fact show that for any BSB1 Σ, V B φ (Σ) is the same, and this is the unique BSB1 fixed point of Eq. ( 43).

Theorem F.5. The BSB1 fixed point Σ * = µ * I + ν * 11 T to Eq. ( 43) is unique and satisfies

µ * + ν * = Γ( B-1 2 ) Γ( B-2 2 ) √ π π 0 dθ 1 sin B-3 θ 1 φ( √ Bζ 1 (θ 1 )) 2 ν * = B-3 2π π 0 dθ 1 π 0 dθ 2 sin B-3 θ 1 sin B-4 θ 2 φ( √ Bζ 1 (θ 1 ))φ( √ Bζ 2 (θ 1 , θ 2 )), if B ≥ 4 1 2π 2π 0 dθφ(- √ 2 sin θ)φ(sin θ/ √ 2 -cos θ √ 3/ √ 2), if B = 3

where

ζ 1 (θ) = - B -1 B cos θ ζ 2 (θ 1 , θ 2 ) = 1 B(B -1) cos θ 1 - B -2 B -1 sin θ 1 cos θ 2

Proof. For any BSB1 Σ, e T Σe = µI B-1 for some µ. Thus we can apply Proposition E.33 to GΣG.

We observe that K(v; GΣG) = 1, so that

V B φ (Σ) = V B φ (GΣG) = E v∼S B-2 φ( √ Bev) ⊗2 .

Note that this is now independent of the specific values of Σ. Now, using the specific form of e given by Eq. ( 33), we see that V B φ (Σ) 11 only depends on θ 1 and V B φ (Σ) 12 only depends on θ 1 and θ 2 . Then applying Proposition E.33 followed by Lemma E.36 and Lemma E.37 yield the desired result.

F.2.2 LAPLACE METHOD:

In the case that φ is positive-homogeneous, we can apply Laplace's method Appendix E.1 and obtain the BSB1 fixed point in a much nicer form.

Definition F.6. For any α ≥ 0, define K α,B :=

c α P B-1 2 , α -1 B-1 2 α where P (a, b) := Γ(a + b)/Γ(a)

is the Pochhammer symbol.

Note that Proposition F.7. lim B→∞ K α,B = c α .

We give a closed form description of the BSB1 fixed point for a positive homogeneous φ in terms of its J function. Theorem F.8. Suppose φ : R → R is degree α positive-homogeneous. For any BSB1 Σ ∈ S B , V B φ (Σ) is BSB1. The diagonal entries are K α,B J φ (1) and the off-diagonal entries are

K α,B J φ -1 B-1 .

Here c α is as defined in Defn E.7 and J φ is as defined in Defn E.20. Thus a BSB1 fixed point of Eq. (43) exists and is unique.

Proof. Let e be an B × (B -1) matrix whose columns form an orthonormal basis of im G := {Gv : a, b ) denote a matrix with a on the diagonal and b on the off-diagnals. Note that V φ is positive homogeneous of degree α, so for any µ,

v ∈ R B } = {w ∈ R B : i w i = 0}. Let BSB1(

V φ (eµI B-1 (I B-1 + 2sµI B-1 ) -1 e T ) = V φ ( µ 1 + 2sµ ee T ) = µ 1 + 2sµ α V φ (ee T )

by Proposition E.24

= µ 1 + 2sµ α V φ (G) = µ 1 + 2sµ B -1 B α V φ BSB1 1, -1 B -1 = µ 1 + 2sµ B -1 B α c α J φ BSB1 1, -1 B -1 by Corollary E.19

So by Eq. ( 26),

V B φ (eµIe T ) = B α Γ(α) -1 ∞ 0 ds s α-1 (1 + 2sµ) -(B-1)/2 µ 1 + 2sµ B -1 B α c α J φ BSB1 1, -1 B -1 = Γ(α) -1 (B -1) α c α J φ BSB1 1, -1 B -1 ∞ 0 ds s α-1 (1 + 2sµ) -(B-1)/2-α µ α = Γ(α) -1 (B -1) α c α J φ BSB1 1, -1 B -1 Beta(α, (B -1)/2)2 -α = c α Γ((B -1)/2) Γ(α + (B -1)/2) B -1 2 α J φ BSB1 1, -1 B -1

Corollary F.9. Suppose φ : R → R is degree α positive-homogeneous. If Σ * is the BSB1 fixed point of Eq. (43) as given by Thm F.8, then

EQUATION

)

By setting α = 1 and φ = relu, we get

Corollary F.10. For any BSB1 Σ ∈ S B , V B relu (Σ) is BSB1 with diagonal entries 1 2 and off-diagonal entries 1 2 J 1 -1 B-1 , so that G ⊗2 {V B relu (Σ)} = G(V B relu (Σ))G = 1 2 -1 2 J 1 -1 B-1 G.

By setting α = 1 and φ = id = x → relu(x)relu(-x), we get

Corollary F.11. For any BSB1 Σ ∈ S B , V B id (Σ) is BSB1 with diagonal entries 1 and off-diagonal entries -1 B-1 , so that G ⊗2 {V B id (Σ)} = B/(B -1)G.

Remark F.12. One might hope to tweak the Laplace method for computing the fixed point to work for the Fourier method, but because there is no nice relation between V φ (cΣ) and V φ (Σ) in general, we cannot simplify Eq. (32) as we can Eq. ( 26) and Eq. ( 27).

F.2.3 GEGENBAUER EXPANSION:

It turns out that the BSB1 fixed point can be described very cleanly using the Gegenbauer coefficients of φ. -2 , so that by Proposition E.34,

When Σ is BSB1, K(v, Σ) = 1 for all v ∈ S B

Σ * = E v∼S B-2 φ( √ Bev) ⊗2 Σ * ab = E v∼S B-2 φ( √ Be a,: v)φ( √ Be b,: v) for any a, b ∈ [B]

, independent of the actual values of Σ. Here e a,: is the ath row of e.

Theorem F.13. If φ( √ B -1x) ∈ L 2 ((1 -x 2 ) B-3 2 ) has Gegenbauer expansion ∞ l=0 a l 1 c B-1,l C ( B-3 2 ) l (x), then Σ * = µ * I + ν * 11 T with µ * = ∞ l=0 a 2 l 1 c B-1,l C ( B-3 2 ) l (1) -C ( B-3 2 ) l -1 B -1 = ∞ l=0 a 2 l 1 c B-1,l B -4 + l l -C ( B-3 2 ) l -1 B -1 ν * = ∞ l=0 a 2 l 1 c B-1,l C ( B-3 2 ) l -1 B -1

Proof. By the reproducing property Fact E.46, we have

Σ * ab = φ( √ Be a,: •), φ( √ Be b,: •) S B-2 = φ( √ B -1ê a,: •), φ( √ B -1ê b,: •) S B-2 = ∞ l=0 a l 1 c B-1,l C ( B-3 2 ) l (ê a,: •), ∞ l=0 a l 1 c B-1,l C ( B-3 2 ) l (ê b,: •) S B-2 = ∞ l=0 a l Z B-2,(l) êa,: , ∞ l=0 a l Z B-2,(l) êb,: S B-2 = ∞ l=0 a 2 l Z B-2,(l) êa,: (ê b,: ) = ∞ l=0 a 2 l 1 c B-1,l C ( B-3 2 ) l ( êa,: , êb,: ) = ∞ l=0 a 2 l 1 c B-1,l C ( B-3 2 ) l 1 if a = b -1/(B -1) else

This gives the third equation. Since µ * = Σ * 11 -Σ * 12 , straightforward arithmetic yields the first equation. The second follows from Fact E.41.

We will see later that the rate of gradient explosion also decomposes nicely in terms of Gegenbauer coefficients (Thm G.5). However, this is not quite true for the eigenvalues of the forward convergence (Thms F.22 and F.24).

F.3 LOCAL CONVERGENCE:

In this section we consider linearization of the dynamics given in Eq. ( 43). Thus we must consider linear operators on the space of PSD linear operators S B . To avoid confusion, we use the following notation: If T : S B → S B (for example the Jacobian of V B φ ) and Σ ∈ S B , then write T {Σ} for the image of Σ under T .

A priori, the Jacobian of V B φ at its BSB1 fixed point may seem like a very daunting object. But a moment of thought shows that Proposition F.14. The Jacobian dV B φ dΣ Σ : H B → H B is ultrasymmetric for any BSB1 Σ. Now we prepare several helper reults in order to make progress understanding batchnorm.

Definition F.15. Define n(x) = √ Bx/ x , i.e. division by sample standard deviation.

Batchnorm B φ can be decomposed as the composition of three steps, φ • n • G, where G is meancentering, n is division by standard deviation, and φ is coordinate-wise application of nonlinearity.

We have, as operators H → H,

dV B φ (Σ) dΣ = d E[B φ (z) ⊗2 : z ∈ N (0, Σ)] dΣ = d E[(φ • n)(x) ⊗2 : x ∈ N (0, GΣG)] dΣ = dV (φ•n) (Σ G ) dΣ G • dΣ G dΣ = dV (φ•n) (Σ G ) dΣ G • G ⊗2

Definition F.16. With Σ * being the unique BSB1 fixed point, write

U := G ⊗2 • dV (φ•n) (Σ G ) dΣ G Σ G =GΣ * G : H G → H G .

It turns out to be advantageous to study U first and relate its eigendecomposition back to that of dV B φ (Σ) dΣ Σ=Σ * , where Σ * is the BSB1 fixed point, by applying Lemma F.17. Lemma F.17. Let X and Y be two vector spaces. Consider linear operators A : X → Y, B : Y → X. Then 1. rank AB = rank BA 2. If v ∈ Y is an eigenvector of AB with nonzero eigenvalue, then X Bv = 0 and Bv is an eigenvector of BA of the same eigenvalue 3. Suppose AB has k = rank AB linearly independent eigenvectors {v i } k i=1 of nonzero eigenvalues {λ i } k i=1 . Then BA has k linearly independent eigenvectors {Bv i } k i=1 with the same eigenvalues {λ i } k i=1 which are all eigenvectors of BA with nonzero eigenvalues, up to linear combinations within eigenvectors with the same eigenvalue.

With A = G ⊗2 and B = dV (φ•n) (Σ G ) dΣ G

, Thm E.61 implies that AB and BA can both be diagonalized, and this lemma implies that all nonzero eigenvalues of dV B φ (Σ) dΣ can be recovered from those of U.

Proof. (Item 1) Observe rank AB = rank ABAB ≤ rank BA. By symmetry the two sides are in fact equal.

(Item 2) Bv cannot be zero or otherwise ABv = A0 = 0, contradicting the fact that v is an eigenvector with nonzero eigenvalue. Suppose λ is the eigenvalue associated to v. Then BA(Bv) = B(ABv) = B(λv) = λBv, so Bv is an eigenvector of BA with the same eigenvalue.

(Item 3) Item 2 shows that {Bv i } k i=1 are eigenvectors BA with the same eigenvalues {λ i } k i=1 . The eigenspaces with different eigenvalues are linearly independent, so it suffices to show that if {Bv ij } j are eigenvectors of the same eigenvalue λ s , then they are linearly independent. But j a j Bv ij = 0 =⇒ j a j v ij = 0 because B is injective on eigenvectors by Item 2, so that a j = 0 identically. Hence {Bv j } j is linearly independent.

Since rank BA = k, these are all of the eigenvectors with nonzero eigenvalues of BA up to linearly combinations.

Lemma F.18. Let f : R B → R A be measurable, and Σ ∈ S B be invertible. Then for any Λ ∈ R B×B , with •, • denoting trace inner product,

d dΣ E[f (z) : z ∼ N (0, Σ)]{Λ} = 1 2 E[f (z) Σ -1 zz T Σ -1 -Σ -1 , Λ : z ∼ N (0, Σ)]

If f is in addition twice-differentiable, then

d dΣ E[f (z) : z ∼ N (0, Σ)] = 1 2 E d 2 f (z) dz 2 : z ∼ N (0, Σ)

whenever both sides exist.

Proof. Let Σ t , t ∈ (-, ) be a smooth path in S B , with Σ 0 = Σ. Write d dt Σ t = Σt . Then

d dt E[f (z) : z ∼ N (0, Σ t )] = d dt (2π) -B/2 det Σ -1/2 t dz e -1 2 z T Σ -1 t z f (z) t=0 = (2π) -B/2 -1 2 det Σ -1/2 0 tr Σ -1 0 Σ0 dz e -1 2 z T Σ -1 0 z/2 f (z) + (2π) -B/2 det Σ -1/2 0 dz e -1 2 z T Σ -1 0 z -1 2 z T (-Σ -1 0 Σ0 Σ -1 0 )zf (z) = -1 2 (2π) -B/2 det Σ -1/2 0 dz e -1 2 z T Σ -1 0 z f (z)(tr Σ -1 0 Σ0 -z T Σ -1 0 Σ0 Σ -1 0 z) = 1 2 (2π) -B/2 det Σ -1/2 dz e -1 2 z T Σ -1 z f (z) Σ -1 zz T Σ -1 -Σ -1 , Σ0 = -1 2 E[f (z) Σ -1 -Σ -1 zz T Σ -1 , Σ0 : z ∼ N (0, Σ)] Note that v T d dz e -1 2 z T Σ -1 z = -v T Σ -1 ze -1 2 z T Σ -1 z w T d 2 dz 2 e -1 2 z T Σ -1 z v = w T Σ -1 zz T Σ -1 v -w T Σ -1 v e -1 2 z T Σ -1 z = Σ -1 zz T Σ -1 -Σ -1 , vw T e -1 2 z T Σ -1 z

so that as a cotensor,

d 2 dz 2 e -1 2 z T Σ -1 z {Λ} = Σ -1 zz T Σ -1 -Σ -1 , Λ e -1 2 z T Σ -1 z for any Λ ∈ R B×B . Therefore, d E[f (z) : z ∼ N (0, Σ)] dΣ {Λ} = 1 2 (2π) -B/2 det Σ -1/2 dz f (z) d 2 e -1 2 z T Σ -1 z d 2 z {Λ} = 1 2 (2π) -B/2 det Σ -1/2 dz e -1 2 z T Σ -1 z d 2 f (z) d 2 z {Λ}

(by integration by parts)

= 1 2 E d 2 f (z) d 2 z : z ∼ N (0, Σ) {Λ} Note that for any Λ ∈ H B with Λ op < µ * , E[φ • n(z) ⊗2 : z ∈ N (0, µ * G + Λ)] = E[φ • n(Gz) ⊗2 : z ∈ N (0, µ * I + Λ)] = E[B φ (z) ⊗2 : z ∈ N (0, µ * I + Λ)]

so that we have, for any Λ ∈ H B , U{Λ} = dV B φ (Σ) dΣ

Σ=µ * I {Λ} = 1 2 E[B φ (z) ⊗2 µ * -1 zz T -µ * -1 I, Λ : z ∼ N (0, µ * I)] (by Lemma F.18) = 1 2 µ * -2 E[B φ (z) ⊗2 zz T , Λ : z ∼ N (0, µ * I)] - 1 2 µ * -1 Σ * I, Λ

(by Thm F.8)

= (2µ * ) -1 E[B φ (z) ⊗2 zz T , Λ : z ∼ N (0, I)] -Σ * I, Λ (B φ is scale-invariant) = (2µ * ) -1 E[B φ (z) ⊗2 Gzz T G, Λ : z ∼ N (0, I)] -Σ * I, Λ (Λ = GΛG)

Let's extend to all matrices by this formula: Definition F.19. Define

Ũ : R B×B → S B , Λ → (2µ * ) -1 E[B φ (z) ⊗2 Gzz T G, Λ : z ∼ N (0, I)] -Σ * I, Λ F.3.1 SPHERICAL INTEGRATION So Ũ H B = U H B .

Ultimately we will apply Thm E.61 to

G ⊗2 • Ũ H B = G ⊗2 • U Definition F.20. Write τ ij = 1 2 δ i δ T j + δ j δ T i . Then Ũ{τ ij } kl = (2µ * ) -1 (E[B φ (z) k B φ (z) l (Gz) i (Gz) j : z ∼ N (0, I)] -Σ * kl I(i = j)) = (2µ * ) -1 (E[φ(n(y)) k φ(n(y)) l y i y j : y ∼ N (0, G)] -Σ * kl I(i = j)) = (2µ * ) -1 E[φ(n(y)) k φ(n(y)) l y i y j : y ∼ N (0, ee T )] -Σ * kl I(i = j) (See Defn E.3 for defn of e) = (2µ * ) -1 (E[φ(n(ex)) k φ(n(ex)) l (ex) i (ex) j : x ∼ N (0, I B-1 )] -Σ * kl I(i = j))

Here we will realize e as the matrix in Eq. ( 33).

Definition F.21.

Define W ij|kl := E[φ(n(ex)) k φ(n(ex)) l (ex) i (ex) j : x ∼ N (0, I B-1 )].

Then

Ũ{τ ij } kl = (2µ * ) -1 (W ij|kl -Σ * kl I(i = j)).

If we can evaluate W ij|kl then we can use Thm E.61 to compute the eigenvalues of G ⊗2 • U. It's easy to see that W ij|kl is ultrasymmetric Thus WLOG we can take i, j, k, l from {1, 2, 3, 4}.

By Lemma E.35, and the fact that x → ex is an isometry,

W ij|kl = E[r 2 φ( √ Bev) k φ( √ Bev) l (ev) i (ev) j : x ∼ N (0, I B-1 )] = (B -5)(B -3)(B -1)(2π) -2 × π 0 dθ 1 • • • π 0 dθ 4 φ( √ Bev) k φ( √ Bev) l (ev) i (ev) j sin B-3 θ 1 • • • sin B-6 θ 4

If WLOG we further assume that k, l ∈ {1, 2} (by ultrasymmetry), then there is no dependence on θ 3 and θ 4 inside φ. So we can expand (ev) i and (ev) j in trigonometric expressions as in Eq. ( 34) and integrate out θ 3 and θ 4 via Lemma E.29. We will not write out this integral explicitly but instead focus on other techniques for evaluating the eigenvalues.

Thus, if φ( √ B -1x) ∈ L 2 ((1 -x 2 ) B-3 2 ) has Gegenbauer expansion ∞ l=0 a l 1 c B-1,l C ( B-3 2 ) l (x), then by Proposition E.43, x 2 φ( √ B -1x) = ∞ l=0 a l 1 c B-1,l x 2 C ( B-3 2 ) l (x) = ∞ l=0 a l 1 c B-1,l κ 2 (l, B -3 2 )C ( B-3 2 ) l+2 (x) + κ 0 (l, B -3 2 )C ( B-3 2 ) l (x) + κ -2 (l, B -3 2 )C ( B-3 2 ) l-2 (x)

where C

( B-3 2 ) l (x) is understood to be 0 for l < 0 = ∞ l=0 C ( B-3 2 ) l (x) a l-2 c B-1,l-2 κ 2 (l -2, B -3 2 ) + a l c B-1,l κ 0 (l, B -3 2 ) + a l+2 c B-1,l+2 κ -2 (l + 2, B -3 2 )

where

κ i (l, B -3 2 ) is understood to be 0 for l < 0 = ∞ l=0 c -1 B-1,l C ( B-3 2 ) l (x) a l-2 c B-1,l c B-1,l-2 κ 2 (l -2, B -3 2 ) + a l κ 0 (l, B -3 2 ) + a l+2 c B-1,l c B-1,l+2 κ -2 (l + 2, B -3 2 )

We can evaluate, for any a, b (possibly equal),

EQUATION

By Eq. ( 48), the eigenvalue associated with L(B -2, 1) is

B -1 2 µ * -1 ( B B -1 A(1, 1; 1) - 1 B -1 Σ * aa -2 B B -2 A(1, 3; 1) + 2 B -2 Σ * ab ) = B -1 2 µ * -1 ∞ l=0 c -1 B-1,l (γ l C ( B-3 2 ) l (1) -τ l C ( B-3 2 ) l ( -1 B -1 ))

where

γ l := a l a l-2 c B-1,l c B-1,l-2 κ 2 (l -2, B -3 2 ) + a l κ 0 (l, B -3 2 ) + a l+2 c B-1,l c B-1,l+2 κ -2 (l + 2, B -3 2 ) - 1 B -1 a 2 l τ l := 2(B -1) B -2 γ l .

Thus,

Theorem F.22. Suppose φ( √ B -1x) has Gegenbauer expansion ∞ l=0 a l 1 c B-1,l C ( B-3 2 ) l (x).

Then the eigenvalue of U with respect to the eigenspace L is a ratio of quadratic forms

λ ↑ L = ∞ l=0 a 2 l w B-1,l + a l a l+2 u B-1,l ∞ l=0 a 2 l v B-1,l where v B-1,l := c -1 B-1,l C ( B-3 2 ) l (1) -C ( B-3 2 ) l ( -1 B -1 ) w B-1,l := B -1 2 κ 0 l, B -3 2 - 1 B -1 c -1 B-1,l C ( B-3 2 ) l (1) - 2(B -1) B -2 C ( B-3 2 ) l ( -1 B -1 ) = l(B -3 + l)(B -3) (B -5 + 2l)(B -1 + 2l) c -1 B-1,l C ( B-3 2 ) l (1) - 2(B -1) B -2 C ( B-3 2 ) l ( -1 B -1 ) = l(B -3 + l)(B -3 + 2l) (B -5 + 2l)(B -1 + 2l) C ( B-3 2 ) l (1) - 2(B -1) B -2 C ( B-3 2 ) l ( -1 B -1 ) u B-1,l := B -1 2 c -1 B-1,l+2 κ -2 l + 2, B -3 2 C ( B-3 2 ) l (1) - 2(B -1) B -2 C ( B-3 2 ) l ( -1 B -1 ) + c -1 B-1,l κ 2 l, B -3 2 C ( B-3 2 ) l+2 (1) - 2(B -1) B -2 C ( B-3 2 ) l+2 ( -1 B -1 ) = B -1 2 (B -3 + l)(B -2 + l) (B -3)(B -1 + 2l) C ( B-3 2 ) l (1) - 2(B -1) B -2 C ( B-3 2 ) l ( -1 B -1 ) + (l + 1)(l + 2) (B -3)(B -1 + 2l) C ( B-3 2 ) l+2 (1) - 2(B -1) B -2 C ( B-3 2 ) l+2 ( -1 B -1 ) .

Note that v B-1,0 = w B-1,0 = u B-1,0 = 0, so that there is in fact no dependence on a 0 , as expected since batchnorm is invariant under additive shifts of φ.

We see that

λ ↑ L ≥ 1 iff 0 ≤ ∞ l=1 a 2 l (w B-1,l -v B-1,l ) + a l a l+2 u B-1,l .

This is a quadratic form on the coefficients {a l } l . We now analyze it heuristically and argue that the eigenvalue is ≥ 1 typically when φ( √ B -1x) explodes sufficiently as x → 1 or x → -1; in other words, the more explosive φ is, the less likely it is to induce Eq. ( 43) to converge to a BSB1 fixed point.

Heuristically, C

( B-3 2 ) l+2 ( -1 B-1 ) is negligible compared to C ( B-3

2 ) l+2

(1) for sufficiently large B and l, so that w

B-1,l -v B-1,l ≈ l(B-3+l)(B-3) (B-5+2l)(B-1+2l) -1 c -1 B-1,l C ( B-3 2 ) l+2 (1) = B-7 4 + O(B -2 ) + O(l -2 ) c -1 B-1,l C ( B-3 2 ) l+2

(1) and is positive. For small l, we can calculate

w B-1,1 -v B-1,1 = -2B 1 + B < 0 w B-1,2 -v B-1,2 = B(B + 1)(B 2 -9B + 12) 2(B -1)(B + 3) ≥ 0, ∀B ≥ 10 w B-1,3 -v B-1,3 = (B -3)(B + 3)B 2 (2B 3 -19B 2 + 16B + 13) 6(B -1) 2 (B + 1)(B + 5) ≥ 0, ∀B ≥ 10.

In fact, plotting w B-1,l -v B-1,l for various values of l suggests that for B ≥ 10, w B-1,l -v B-1,l ≥ 0 for all l ≥ 1 (Fig. 10 ).

Thus, the more "linear" φ is, the larger a 2 1 is compared to the rest of {a l } l , and the more likely that the eigenvalue is < 1. In the case that a l a l+2 = 0 ∀l, then indeed the eigenvalue is < 1 precisely when a 2 1 is sufficiently large. Because higher degree Gegenbauer polynomials explodes more violently as x → ±1, this is consistent with our claim.

These relations allow us to simplify

λ ↑ M = B -1 (B -3)µ * -1 B -2 Σ * 12 + 2(B -1) B -2 Ã12 (11) + (B -1) Ã12 (12) .

By Eq. ( 49), we know the Gegenbauer expansion of Ã12 (11). The following shows the Gegenbauer expansion of Ã12 (12).

EQUATION

2 ), and φ(

√ B -1x) has Gegenbauer expansion ∞ l=0 a l 1 c B-1,l C ( B-3 2 ) l (x). Then Ã12 (12) = E v∼S B-2 φ( √ Be 1,: v)φ( √ Be 2,: v)(e 1,: v)(e 2,: v) = ∞ l=0 α l a 2 l + β l a l a l+2

where

α l = B -1 B ×      l+1 B-1+2l 2 c -1 B-1,l+1 C ( B-3 2 ) l+1 -1 B-1 + B+l-4 B+2l-5 2 c -1 B-1,l-1 C ( B-3 2 ) l-1 -1 B-1 if l ≥ 1 1 B-1 2 c -1 B-1,1 C ( B-3 2 ) 1 -1 B-1 otherwise β l = B -1 B × 2(l + 1)(B + l -2) (B + 2l -1) 2 c -1 B-1,l+1 C ( B-3 2 ) l+1 -1 B -1 . Proof. Let ψ(x) = xφ( √ B -1x).

Then by Proposition E.43, we have

ψ(x) = ∞ l=0 a l c -1 B-1,l 1 B -3 + 2l (l + 1)C ( B-3 2 ) l+1 (x) + (B + l -4)C ( B-3 2 ) l-1 (x) = ∞ l=0 a l-1 l B + 2l -3 + a l+1 B + l -3 B + 2l -3 c -1 B-1,l C ( B-3 2 ) l (x). Then Ã12 (12) = B -1 B E v∼S B-2 ψ(ê 1,: v)ψ(ê 2,: v) = B -1 B ∞ l=0 a l-1 l B + 2l -3 + a l+1 B + l -3 B + 2l -3 2 c -1 B-1,l C ( B-3 2 ) l -1 B -1 .

Rearranging the sum in terms of a l gives the result.

Combining all of the above, we obtain the following Theorem F.24. Suppose φ( √ B -1x) has Gegenbauer expansion

∞ l=0 a l 1 c B-1,l C ( B-3 2 ) l (x).

Then the eigenvalue of U with respect to the eigenspace M is a ratio of quadratic forms

λ ↑ M = ∞ l=0 a 2 l wB-1,l + a l a l+2 ũB-1,l ∞ l=0 a 2 l v B-1,l where ũB-1,l := 2(B -1) 3 (B -3) 2 (B -2)B (B -3 + l)(B -2 + l) B -1 + 2l C ( B-3 2 ) l -1 B -1 + (B -2)(l + 1)(B -2 + l) B -1 + 2l C ( B-3 2 ) l+1 -1 B -1 + (l + 3)(l + 4)(B -3 + 2l) (B + 2l + 1)(B + 2l + 3) C ( B-3 2 ) l+2 -1 B -1 wB-1,l := B -1 (B -3) 2 (B -2) 2(B -1) 2 (B + 2l -3) 2Bl + B + 2l 2 -6l -5 B(B + 2l -5)(B + 2l -1) -(B -3 + 2l) C ( B-3 2 ) l -1 B -1 + (B -2)(B -1) 2 (l + 1) 2 B(B + 2l -1) C ( B-3 2 ) l+1 -1 B -1 + (B -2)(B -1) 2 (B -4 + l) 2 B(B + 2l -5) C ( B-3 2 ) l-1 -1 B -1 where C ( B-3

2 ) -1 (x) = 0 by convention.

Note that ũB-1,l and wB-1,l depends only on C

( B-3 2 ) • -1 B-1

which is much smaller than

C ( B-3 2 ) •

(1) for degree larger than 0. Thus the eigenvalue for M is typically much smaller than the eigenvalue for L (though there are counterexamples, like sin).

F.3.3 LAPLACE METHOD:

Differentiating Eq. ( 26) In what follows, let φ be a positive-homogeneous function of degree α.

We begin by studying

dV (φ•n) (Σ G ) dΣ G

. We will differentiate Eq. ( 26) directly at G ⊗2 {Σ * } = GΣ * G where Σ * is the BSB1 fixed point given by Thm F.8. To that end, consider a smooth path Σ t ∈ S G B , t ∈ (-, ) for some > 0, with Σ 0 = GΣ * G. Set Σ t = e T Σ t e ∈ S B-1 , so that Σ t = eΣ t e T and Σ 0 = µ * I B-1 where µ * = K α,B (J φ (1) -J φ ( -1 B-1 )) as in Thm F.8. If we write Σ and Σ for the time derivatives, we have

• M B has eigenvalue w = c α p α-1 J φ (c) • L B has eigenvalue (u -2v)(B -2) + 2w B = 1 B (B -2) c α αp α-1 J φ (1) -c α p α-1 αJ φ (c) -cJ φ (c) + 2c α p α-1 J φ (c) = c α p α-1 1 B (B -2)α [J φ (1) -J φ (c)] + (2 + c(B -2)) J φ (c) • RG has eigenvalue (B -1)(u -2v) + w B = 1 B (B -1) c α αp α-1 J φ (1) -c α p α-1 αJ φ (c) -cJ φ (c) + c α p α-1 J φ (c) = c α p α-1 1 B (B -1)α [J φ (1) -J φ (c)] + (1 + c(B -1)) J φ (c)

We record the following consequence which will be used frequently in the sequel.

Theorem F.31. Let φ : R → R be positive-homogeneous of degree α.

Then G ⊗2 • dV φ dΣ Σ=G : H G B → H G B has the following eigendecomposition: • M B has eigenvalue λ G,φ M (B, α) := c α B-1 B α-1 J φ -1 B-1 • L B has eigenvalue λ G,φ L (B, α) := c α B-1 B α-1 1 B (B -2)α J φ (1) -J φ -1 B-1 + B B-1 J φ -1 B-1 • RG has eigenvalue λ G,φ G (B, α) := c α B-1 B α α J φ (1) -J φ -1 B-1

Proof. Plug in p = B-1 B and c = -1/(B -1) for Thm F.30. Theorem F.32. Let φ : R → R be a degree α positive-homogeneous function. Then for p = 0, c ∈ R, dV φ dΣ Σ=BSB1(p,cp) : H B → H B has the following eigendecomposition:

• M B has eigenvalue c α p α-1 J φ (c) • L B (a, b) has eigenvalue c α αp α-1 J φ (1) where a = 2 J φ (c) -αJ φ (1) and b = αJ φ (c) -cJ φ (c)

Proof. Use Thm E.72 with the computations from the proof of Thm F.30 as well as the following computation

w -u = c α p α-1 J φ (c) -c α αp α-1 J φ (1) = c α p α-1 J φ (c) -αJ φ (1) v = 1 2 c α p α-1 αJ φ (c) -cJ φ (c)

Theorem F.33. Let φ be positive-homogeneous with degree α. Assume B-1 2 > α. The operator

U = G ⊗2 • dV (φ•n) (Σ G ) dΣ G Σ G =µ * G : H G B → H G B has 3 distinct eigenvalues.

They are as follows:

1

. λ ↑ G (B, α) := 0 with dimension 1 eigenspace RG. 2. λ ↑ M (B, α) := B α µ * -1 2 -α P B + 1 2 , α -1 λ G,φ M (B, α) = c α B(B -1) α-1 µ * -1 2 -α P B + 1 2 , α -1 J φ -1 B -1 with dimension B(B-3) 2 eigenspace M. 3. λ ↑ L (B, α) := B α µ * -1 2 -α P B + 1 2 , α -1 λ G,φ L (B, α) = c α (B -1) α-1 µ * -1 2 -α P B + 1 2 , α -1 (B -2)α J φ (1) -J φ -1 B -1 + B B -1 J φ -1 B -1 with dimension B -1 eigenspace L. Proof. Item 1. The case of λ ↑ G (B, α). Since (φ • n)(µ * G) = (φ • n)(µ * G + ωG) for any ω, the operator dV (φ•n) (Σ G ) dΣ G Σ G =µ * G sends G to 0. Item 2. The case of λ ↑ M (B, α).

Assume Σ0 ∈ M, Thm F.31 gives (with denoting Hadamard product, i.e. entrywise multiplication)

dV φ (Σ) dΣ Σ=G { Σ0 } = λ G,φ M (B, α) Σ0 = c α B -1 B α-1 J φ -1 B -1 Σ0

Since tr( Σ0 ) = 0, Eq. ( 50) gives

G ⊗2 • dV (φ•n) (Σ G ) dΣ G Σ G =µ * G { Σ0 } = G ⊗2 • d dt V (φ•n) (Σ t ) t=0 { Σ0 } = B α µ * -1 2 -α P B + 1 2 , α -1 G ⊗2 • dV φ (Σ) dΣ Σ=G Σ0 = B α µ * -1 2 -α P B + 1 2 , α -1 c α B -1 B α-1 J φ -1 B -1 G ⊗2 { Σ0 } = c α B(B -1) α-1 µ * -1 2 -α P B + 1 2 , α -1 J φ -1 B -1 Σ0 So the eigenvalue for M is λ ↑ L = c α B(B -1) α-1 µ * -1 2 -α P B+1 2 , α -1 J φ -1 B-1 .

Item 3. The case of λ ↑ L (B, α). Let Σ0 = L(B -2, 1). Thm F.31 gives

G ⊗2 • dV φ (Σ) dΣ Σ=G { Σ0 } = λ G,φ L (B, α) Σ0 where λ G,φ L (B, α) := 1 B c α B -1 B α-1 (B -2)α J φ (1) -J φ -1 B -1 + B B -1 J φ -1 B -1

Since tr( Σ0 ) = 0, Eq. ( 50) gives

G ⊗2 • dV (φ•n) (Σ G ) dΣ G Σ G =µ * G { Σ0 } = G ⊗2 • d dt V (φ•n) (Σ t ) t=0 = B α µ * -1 2 -α P B + 1 2 , α -1 G ⊗2 • dV φ (Σ) dΣ Σ=G Σ0 = B α µ * -1 2 -α P B + 1 2 , α -1 λ Σ0 = (B -1) α-1 c α µ * -1 2 -α P B + 1 2 , α -1 × (B -2)α J φ (1) -J φ -1 B -1 + B B -1 J φ -1 B -1 × Σ0 So λ ↑ L = (B-1) α-1 c α µ * -1 2 -α P B+1 2 , α -1 (B -2)α J φ (1) -J φ -1 B-1 + B B-1 J φ -1 B-1

With some routine computation, we obtain Proposition F.34. With φ and α as above, as B → ∞,

λ ↑ L (B, α) ∼ α + B -1 (2α 3 -4α 2 -α + J φ (0) J φ (1) -J φ (0) ) + O(B -2 ) λ ↑ M (B, α) ∼ J φ (0) J φ (1) -J φ (0) + B -1 2α 2 -4α + 1 J φ (0) J φ (1) -J φ (0) - J φ (0) J φ (1) -J φ (0) - J φ (0) J φ (1) -J φ (0) 2 + O B -2

We can compute, by Proposition E.23,

J φ (0) J φ (1) -J φ (0) = (a + b) 2 1 √ π Γ( α 2 +1) 2 Γ(α+ 1 2 ) (a 2 + b 2 ) -(a -b) 2 1 2 √ π Γ( α 2 + 1 2 ) 2 Γ(α+ 1 2 ) ≤ 2Γ α 2 + 1 2 √ πΓ α + 1 2

where the last part is obtained by optimizing over a and b. On α ∈ (-1/2, ∞), this is greater than α iff α < 1. Thus for α ≥ 1, the maximum eigenvalue of U is always achieved by eigenspace L for large enough B.

G BACKWARD DYNAMICS:

Let's extend the definition of the V operator:

of Eq. (43), the eigenvalue of

G ⊗2 • F(Σ * ) corresponding to eigenspace RG is λ ↓ G = µ * -1 ∞ l=0 a 2 l l + B -3 B -3 lc -1 B-1,l C ( B-3 2 ) l (1) -C ( B-3 2 ) l -1 B -1 = ∞ l=0 a 2 l l+B-3 B-3 lc -1 B-1,l C ( B-3 2 ) l (1) -C ( B-3 2 ) l -1 B-1 ∞ l=0 a 2 l c -1 B-1,l C ( B-3 2 ) l (1) -C ( B-3 2 ) l -1 B-1 > 1 This is minimized (over choices of φ) iff φ is linear, in which case λ ↓ G = B-2 B-3 .

G.2 LAPLACE METHOD:

In this section suppose that φ is degree α positive-homogeneous. Set D = Diag(φ (y)) (and recall v = y/ y , r = y ). Then D is degree α -1 positive-homogeneous in x (because φ is).

Consequently we can rewrite Eq. ( 52) as follows,

EQUATION

) where Σ * is the BSB1 fixed point of Eq. ( 43), GΣ * G = µ * G, and

A = E r -2α D ⊗2 : y ∼ N (0, µ * G) B = E r -2(α+2) (Dyy T ) ⊗2 : y ∼ N (0, µ * G) C = E r -2(α+1) D ⊗ (Dyy T ) + (Dyy T ) ⊗ D : y ∼ N (0, µ * G)

Each term in the sum above is ultrasymmetric, so has the same eigenspaces RG, M, L (this will also become apparent in the computations below without resorting to Thm E.61). We can thus compute the eigenvalues for each of them in order.

In all computation below, we apply Lemma E.2 first to relate the quantity in question to V φ .

Computing A. For two matrices Σ, Λ, write Σ Λ for entrywise multiplication. We have by Lemma E.2

A{Λ} = E r -2α DΛD : y ∼ N (0, µ * G) = E r -2α Λ φ (y)φ (y) T : y ∼ N (0, µ * G) = Γ(α) -1 ∞ 0 ds s α-1 det(I + 2sµ * G) -1/2 Λ V φ µ * 1 + 2sµ * G = Γ(α) -1 ∞ 0 ds s α-1 (1 + 2sµ * ) -(B-1)/2 Λ µ * 1 + 2sµ * α-1 V φ (G) = Γ(α) -1 Λ V φ (G) ∞ 0 ds (µ * s) α-1 (1 + 2sµ * ) -(B-1)/2-α+1 = Γ(α) -1 Λ V φ (G) Beta B -3 2 , α µ * -1 2 -α = Λ c α-1 B -1 B α-1 BSB1 J φ (1), J φ -1 B -1 P B -3 2 , α -1 µ * -1 2 -α

Then Thm E.73 gives the eigendecomposition for

EQUATION

ϕ is degree α + 1 positive-homogeneous. Thus,

B{Λ} = Γ(α + 2) -1 ∞ 0 ds s α+1 det(I + 2sµ * G) -1/2 V (4) ψ µ * 1 + 2sµ * G {Λ} = Γ(α + 2) -1 ∞ 0 ds s α+1 (1 + 2sµ * ) -(B-1)/2 µ * 1 + 2sµ * α+1 V (4) ψ (G){Λ} = Γ(α + 2) -1 V (4) ψ (G){Λ} ∞ 0 ds (µ * s) α+1 (1 + 2sµ * ) -(B-1)/2-α-1 = Γ(α + 2) -1 V (4) ψ (G){Λ} Beta( B -3 2 , α + 2)2 -α-2 µ * -1 = P B -3 2 , α + 2 -1 2 -α-2 µ * -1 V (4) ψ (G){Λ}.

By Proposition G.8, V

EQUATION

has the following eigendecomposition (note that here µ * is still with respect to φ, not ψ = αφ)

1. Eigenspace RG with eigenvalue λ B G := α 2 B -3 2. Eigenspace L with eigenvalue λ B L := B α P B -3 2 , α + 2 -1 2 -α-2 µ * -1 2λ G,ψ L (B, α) = 2α 3 (B -2) (B -3)(B -1)(B -1 + 2α) + 2α 2 BJ φ -1 B-1 (B -3)(B -1) 2 (B -1 + 2α)(J φ (1) -J φ -1 B-1 ) 3. Eigenspace M with eigenvalue λ B M := B α P B -3 2 , α + 2 -1 2 -α-2 µ * -1 2λ G,ψ M (B, α) = 2α 2 BJ φ -1 B-1 (B -3)(B -1)(B -1 + 2α)(J φ (1) -J φ -1 B-1 )

Proof. The only thing to justify is the value of λ B G . We have Proof. Let Π t , t ∈ (-, ) be a smooth path in H B with Π 0 = I. Write D t = Diag(φ (Π t y)), so that D 0 = D. Then, using Πt to denote t derivative,

λ B G = B α P B -3 2 , α + 2 -1 2 -α-2 µ * -1 × (B -1)c α B -1 B α J ψ (1) -J ψ -1 B -1 + 2λ G,ψ M (B, α) = B α P B -3 2 , α + 2 -1 2 -α-2 µ * -1 (B -1 + 2α)c α B -1 B α J ψ (1) -J ψ -1 B -1 = c α (B -1) α (B -1 + 2α)P B -3 2 , α + 2 -1 2 -α-2 µ * -1 J ψ (1) -J ψ -1 B -1 = c α (B -1) α P B -3 2 , α + 1 -1 2 -α-1 µ * -1 J ψ (1) -J ψ -1 B -1 = c α (B -1) α P B -3 2 , α + 1 -1 2 -α-1 K -1 α,B α 2 (since J ψ = α 2 J φ ) = α 2 c α (B -1) α P B-3 2 , α + 1 -1 2 -α-1 c α P B-1 2 , α -1 B-1 2 α (Defn F.6) = α 2 2 -1 (B -3)/2 = α 2 B -3 Computing C. By Lemma E.2, C{Λ} = E r -2(α+1) DΛyy T D + Dyy T ΛD : y ∼ N (0, µ * G) = Γ(α + 1) -1 ds s α det(I + 2sµ * G) -1/2 E[DΛyy T D + Dyy T ΛD : y ∼ N (0, µ * 1 + 2sµ * G)] = Γ(α + 1) -1 ds s α (1 + 2sµ * ) -(B-1)/2 µ * 1 + 2sµ * α E[DΛyy T D + Dyy T ΛD : y ∼ N (0, G)] = Γ(α + 1) -1 ds (µ * s) α (1 + 2sµ * ) -(B-1)/2-α E[DΛyy T D + Dyy T ΛD : y ∼ N (0, G)] = Γ(α + 1) -1 Beta B -3 2 , α + 1 2 -α-1 µ * -1 E[DΛyy T D + Dyy T ΛD : y ∼ N (0, G)] = P B -3 2 , α + 1 -1 2 -α-1 µ * -1 E[DΛyy T D + Dyy T ΛD : y ∼ N (0, G)] Lemma G.

d dt V φ (Π t ΣΠ t ) = d dt E[φ(Π t y)φ(Π t y) T : y ∼ N (0, Σ)] = E[D t Πt yφ(Π t y) T + φ(Π t y)y T Πt D t : y ∼ N (0, Σ)] d dt V φ (Π t ΣΠ t ) t=0 = E[D Π0 yφ(y) T + φ(y)y T Π0 D : y ∼ N (0, Σ)]

Because for x ∈ R, αφ(x) = xφ (x), for y ∈ R B we can write φ(y) = α -1 Diag(φ (y))y. Then 

α dV φ (Σ) dΣ {Σ Π0 + Π0 Σ} = α dV φ (Σ) dΣ dΠ t ΣΠ t dt = α d dt V φ (Π t ΣΠ t ) t=0 = E[D Π0 yy T D + Dyy T Π0 D : y ∼ N (0, Σ)] Therefore, for Λ ∈ H G B , C{Λ} = P B -3 2 , α + 1 -1 2 -α-1 µ * -1 α dV φ (Σ) dΣ Σ=G {GΛ + ΛG} = 2αP B -3 2 , α + 1 -1 2 -α-1 µ * -1 dV φ (Σ)

C G := 2B α αP B -3 2 , α + 1 -1 2 -α-1 µ * -1 λ G,φ G (B, α) = 2α 2 B -3 2. eigenspace L with eigenvalue λ C L := 2B α αP B -3 2 , α + 1 -1 2 -α-1 µ * -1 λ G,φ L (B, α) = 2α 2 (B -2) (B -3)(B -1) + 2αBJ φ -1 B-1 (B -3)(B -1) 2 (J φ (1) -J φ -1 B-1 ) 3. eigenspace M with eigenvalue λ C M := 2B α αP B -3 2 , α + 1 -1 2 -α-1 µ * -1 λ G,φ M (B, α) = 2αBJ φ -1 B-1 (B -3)(B -1)(J φ (1) -J φ -1 B-1

Altogether, by Eq. ( 53) and Thms G.6, G.9 and G.11, this implies Theorem G.12. G ⊗2 • F(Σ * ) : H G B → H G B has eigenspaces RG, M, L respectively with the following eigenvalues

λ ↓ G = (B -3 + 2α) (2α -1)J φ -1 B-1 + α 2 (B -1)J φ (1) (2α -1)(B -3)(B -1)(J φ (1) -J φ -1 B-1 ) - α 2 B -3 = α 2 (B -3 + 2α) J φ (1) + 1 B-1 J α -1 φ -1 B-1 (2α -1)(B -3)(J φ (1) -J φ -1 B-1 ) - α 2 B -3 = α 2 (B -2)J φ (1) + (B -3 + 2α) 1 B-1 J α -1 φ -1 B-1 + (2α -1)(J φ (1) -J φ -1 B-1 (2α -1)(B -3)(J φ (1) -J φ -1 B-1 ) λ ↓ L = - 2α 2 (B -2)(B -1 + α) (B -3)(B -1)(B -1 + 2α) + 2 α 2 (3B -4) + α 3B 2 -11B + 8 + (B -3)(B -1) 2 (B -3)(B -1) 2 (B -1 + 2α) J φ -1 B-1 J φ (1) -J φ -1 B-1 + α 2 (B -2)(B -3 + 2α) (2α -1)(B -3)(B -1) J φ (1) J φ (1) -J φ -1 B-1 λ ↓ M = B(B 2 + 2(α -2)B + 2(α -3)α + 3) (B -3)(B -1)(B -1 + 2α) J φ -1 B-1 J φ (1) -J φ -1 B-1

.

H CROSS BATCH: FORWARD DYNAMICS:

In this section, we study the generalization of Eq. ( 43) to multiple batches. Definition H.1. For linear operators T i : X i → Y i , i = 1, . . . , k, we write i T i : i X i → i Y i for the operator (x 1 , . . . , x k ) → (T 1 (x 1 ), . . . , T k (x k )). We also write T ⊕n 1 := n j=1 T 1 for the direct sum of n copies of T 1 .

For k ≥ 2, now consider the extended ("k-batch") dynamics on Σ ∈ S kB defined by

Σ l = V B ⊕k φ ( Σ l-1 ) (54) = E[(B φ (h 1:B ), B φ (h B+1:2B ), . . . , B φ (h (k-1)B+1:kB )) ⊗2 : h ∼ N (0, Σ l-1

)] If we restrict the dynamics to just the upper left B × B submatrix (as well as any of the diagonal B × B blocks) of Σ l , then we recover Eq. ( 43).

H.1 LIMIT POINTS:

In general, like in the case of Eq. ( 43), it is difficult to prove global convergence behavior. Thus we manually look for fixed points and the local convergence behaviors around them. A very natural extension of the notion of BSB1 matrices is the following Definition H.2. We say a matrix Σ ∈ S kB is CBSB1 (short for "1-Step Cross-Batch Symmetry Breaking") if Σ in block form (k × k blocks, each of size B × B) has one common BSB1 block on the diagonal and one common constant block on the off-diagonal, i.e.

     BSB1(a, b) c11 T c11 T • • • c11 T BSB1(a, b) c11 T • • • c11 T c11 T BSB1(a, b) • • • . . . . . . . . . . . .     

We will study the fixed point Σ * to Eq. ( 54) of CBSB1 form. 

c * = Γ((B-1)/2) Γ((B-2)/2) π -1/2 π 0 dθ φ(- √ B -1 cos θ) sin B-3 θ. If φ( √ B -1x) has Gegenbauer expansion ∞ l=0 a l 1 c B-1,l C ( B-3 2 ) l (x), then c * = a 2 0 .

Proof. We will prove for k = 2. The general k cases follow from the same reasoning.

Let Σ = Σ c11 T c11 T Σ where Σ = BSB1(a, b). As remarked below Eq. ( 54), restricting to any diagonal blocks just recovers the dynamics of Eq. ( 43), which gives the claim about the diagonal blocks being Σ * through Thm F.5.

We now look at the off-diagonal blocks. where

Σ G := G ⊕2 ΣG ⊕2 = G 0 0 G Σ G 0 0 G = GΣG cG11 T G cG11 T G GΣG = (a -b)G 0 0

(ab)G (the last step follows from Lemma E.52). Thus y 1:B is independent from y B+1:2B , and where in the last line, (y, z) ∼ N (0, (Σ -1 + 2D 2 ) -1 ). We recover the equation in question with the following simplifications. T { Σ} ij = V{ Σ ii } + V{ Σ jj } T + W{ Σ ij }, ∀i = j, where U : H B → H B , V : H B → R B×B , and W : R B×B → R B×B are linear operators, and T { Σ} ij denotes the (i, j)th block of T { Σ}. We say that T is blockwise diagonal-off-diagonal semidirect, or BDOS for short. We write more specifically T = BDOS(U, V, W).

V B ⊕2

We have F = BDOS(U, V, W) where

EQUATION

)

Note that the factor 2 in W is due to the contribution of ( Σ G ) 21 by symmetry. We will first show that W is multiplication by a constant. φ (n(ey)) ⊗ φ (n(ey )) (y T e T Λeytr Λ)

= 1 2 µ * -1 √ 2P B -1 2 , 1 2 2 E v,v ∼S B-2 φ √ Bev ⊗ φ √ Bev v T e T Λev -(tr Λ) E v,v ∼S B-2 φ √ Bev ⊗ φ √ Bev

by Proposition E.33.

Note that (ê b,:

E v,v ∼S B-2 φ √ Bev ⊗ φ √ Bev = E v∼S B-2 φ √ Bev ⊗2 is a constant matrix,

) = a 1 B -1 1 c B-1,1 C ( B-3 2 ) 1 ( êb,: , êc,: ) = a 1 if b = c -a1 B-1

otherwise.

Thus, Therefore,

H{Λ} ab = B -1 B -a 1 B -1 2 c,d∈[B] Λ cd + a 1 - -a 1 B -1 -a 1 B -1 c=b XOR d=a Λ cd + a 2 1 - a 2 1 (B -1) 2 Λ ab = B -1 B -2 a 1 - -a 1 B -1 -a 1 B -1 Λ ab + a

W{Λ} = 2G ⊗2 • J{Λ} = 2a 2 1 µ * -1 P B -1 2 , 1 2 2 B B -1 Λ = a 2 1 1 2 P B 2 , 1 2 -2 B(B -1) ∞ l=0 a 2 l 1 c B-1,l C ( B-3 2 ) l (1) -C ( B-3 2 ) l -1 B-1 Λ

where we used the identity P B-1 2 , 1 2 2

B-1 = P B 2 , 1 2 -1 .

Next, we step back a bit to give a full treatment of the eigen-properties of BDOS operators. Then we specialize back to our case and give the result for F.

Definition H.9. Let H G kB := { Σ ∈ H kB : G ⊕k ΣG ⊕k = Σ} be the subspace of H kB stable under conjugation by G ⊕k . Definition H.10. Let M kB ∈ H G kB be the space of block matrices Σ ∈ H G kB where the diagonal blocks Σ ii are all zero. Likewise, let M kB ∈ H kB be the space of block matrices Σ ∈ H kB where the diagonal blocks Σ ii are all zero. We suppress the subscript kB when it is clear from the context. 

  :

BDOS operators of the form BDOS(U, V, w), where w stands for multiplication by a constant w, like F, have simple eigendecomposition, just like DOS operators. Theorem H.12. If V -wI B is not singular, then BDOS(V, U, w) : H kB → H kB has the following eigenspaces and eigenvectors 1. M with associated eigenvalue w.

2. For each eigenvalue λ and eigenvector Σ ∈ H B of U, the block matrix

L kB (Σ, -(λ -w) -1 V{Σ}) =    Σ (λ -w) -1 V{Σ} • • • (λ -w) -1 V{Σ} T 0 • • • . . . . . . . . .

along with any matrix obtained from permutating its block columns and rows simultaneously. The associated eigenvalue is λ.

If BDOS(V, U, w) also restricts to H G kB → H G kB , then with respect to this signature, the eigendecomposition is obtained by replacing H kB with H G kB in the above:

1. Block diagonal matrices Σ ∈ H G kB where the diagonal blocks Σ ii are all zero. The associated eigenvalue is w. -1 V{Σ}) along with any matrix obtained from permutating its block columns and rows simultaneously. The associated eigenvalue is λ.

Proof. The proofs are similar in both cases: One can easily verify that these are eigenvectors and eigenvalues. Then by dimensionality considerations, these must be all of them. By Thms H.8 and H.12, we easily obtain From this one sees that V [(B ⊕2 φ ) ] acts independently on each block, and consequently so does its adjoint. The diagonal blocks of Eq. ( 59) evolves according to Eq. ( 51) which we studied in Appendix G. In this section we will study the evolution of the off-diagonal blocks, which is given by 

Ξ l = E (x,

Ξ L-1 = V † {Ξ L } is in M G B .

Thus it suffices to study the eigendecomposition of

G ⊗2 • T † M G B : M G B → M G B .

Below, we will show that its adjoint T • G ⊗2 acting on M G B is in fact multiplication by a constant, and thus so is it.

We first compute 

r -1 (1 -v 2 i )φ ( √ Bv i ) = E rv∼N (0,µ * G) r -1 (1 -v 2 1 )φ ( √ Bv 1 ), ∀i ∈ [B]

and off-diagonal

E rv∼N (0,µ * G) r -1 (-v i v j )φ ( √ Bv i ) = E rv∼N (0,µ * G) r -1 (-v 1 v 2 )φ ( √ Bv 1 ), ∀i = j ∈ [B].

Thus by Lemma E.51, T has two eigenspaces, {x : Gx = x} and R1, with respective eigenvalues T 11 -T 12 and T 11 + (B -1)T 12 . However, because of the composition with G ⊗2 , only the former eigenspace survives with nonzero eigenvalue in G ⊗2 • T .

The rank-1 nature of the BSB2 1 B fixed point under G-projection is reminiscent of the case of B = 2 batchnorm. As we will see, gradients of a certain form will vanish similarly.

Let δ be a gradient vector. Then by Proposition G.4, for Σ * being the BSB2 

G B-1 =     0 0 0 • • • 0 0 1 -1 B-1 -1 B-1 • • • -1 B-1 0 -1 B-1 1 -1 B-1 • • • -1 B-1 0 -1 B-1 -1 B-1 • • • 1 -1 B-1     . Proof. Straightforward computation. Thus if δ = (δ 1 , δ 2 ) T ∈ R 1 × R B-1 , then dB φ dΣ † Σ=Σ * δ = √ Bφ -sgn(r) √ B -1 |r| -1 (0, G B-1 δ 2 ) T , r ∼ N (0, λ * )

If δ 2 = 0, then this is deterministically 0. Otherwise, this random variable does not have finite mean because of the |r| -1 term, and this is solely due to Σ * being rank 1. Thus as long as the gradient δ is not constant in δ 2 , then we expect severe gradient explosion; otherwise we expect severe gradient vanishing.

Remark K.7. We would like to note that, in pytorch, with either Float Tensors or Double Tensors, the gradient empirically vanishes rapidly with generic gradient δ, when φ is rapidly increasing and leads to a BSB2 1 B fixed point, e.g. φ(x) = relu(x) 30 . This is due to numerical issues rather than mathematical issues, where the computation of (G -q⊗2 )δ in backprop does not kill δ 1 completely.

L WEIGHT GRADIENTS:

Recall weight gradient and the backpropagation equation Eq. ( 16) of the main text. 

:

The activation function may be layer dependent, but for ease of exposition we assume that it is not.

Throughout the text, we assume that all elements of the batch are unique.

This definition of V absorbs the previous definitions of V and W inYang & Schoenholz (2017) for the scalar case

Figures and Tables:

2. a 1-dimensional eigenspace RG with eigenvalue λ G,T G . 3. a (B -1)-dimensional eigenspace L := {D G : D diagonal, tr D = 0}, with eigenvalue λ M := {Σ : Σ G = Σ, DiagΣ = 0}, with eigenvalue λ G,T

Figure 1: Numerical confirmation of theoretical predictions. (a,b) Comparison between theoretical prediction (dashed lines) and Monte Carlo simulations (solid lines) for the eigenvalues of the backwards Jacobian (see Thms 3.10 and G.12) as a function of batch size and the magnitude of gradients as a function of depth respectively for rectified linear networks. In each case Monte Carlo simulations are averaged over 200 sample networks of width 1000 and shaded regions denote 1 standard deviation. Dashed lines are shifted slightly for easier comparison. (c,d) Demonstration of the existence of a BSB1 to BSB2 symmetry breaking transition as a function of α for α-ReLU (i.e. the αth power of ReLU) activations. In (c) we plot the empirical variance of the diagonal and off-diagonal entries of the covariance matrix which clearly shows a jump at the transition. In (d) we plot representative covariance matrices for the two phases (BSB1 bottom, BSB2 top). Published as a conference paper at ICLR 2019

Figure A.1: Batch norm leads to a chaotic input-output map with increasing depth. A linear network with batch norm is shown acting on two minibatches of size 64 after random orthogonal initialization. The datapoints in the minibatch are chosen to form a 2d circle in input space, except for one datapoint that is perturbed separately in each minibatch (leftmost datapoint at input layer 0). Because the network is linear, for a given minibatch it performs an affine transformation on its inputs -a circle in input space remains an ellipse throughout the network. However, due to batch norm the coefficients of that affine transformation change nonlinearly as the datapoints in the minibatch are changed. (a) Each pane shows a scatterplot of activations at a given layer for all datapoints in the minibatch, projected onto the top two PCA directions. PCA directions are computed using the concatenation of the two minibatches. Due to the batch norm nonlinearity, minibatches that are nearly identical in input space grow increasingly dissimilar with depth. Intuitively, this chaotic input-output map can be understood as the source of exploding gradients when batch norm is applied to very deep networks, since very small changes in an input correspond to very large movements in network outputs. (b) The correlation between the two minibatches, as a function of layer, for the

Figure A.1: Batch norm leads to a chaotic input-output map with increasing depth. A linear network with batch norm is shown acting on two minibatches of size 64 after random orthogonal initialization. The datapoints in the minibatch are chosen to form a 2d circle in input space, except for one datapoint that is perturbed separately in each minibatch (leftmost datapoint at input layer 0). Because the network is linear, for a given minibatch it performs an affine transformation on its inputs -a circle in input space remains an ellipse throughout the network. However, due to batch norm the coefficients of that affine transformation change nonlinearly as the datapoints in the minibatch are changed. (a) Each pane shows a scatterplot of activations at a given layer for all datapoints in the minibatch, projected onto the top two PCA directions. PCA directions are computed using the concatenation of the two minibatches. Due to the batch norm nonlinearity, minibatches that are nearly identical in input space grow increasingly dissimilar with depth. Intuitively, this chaotic input-output map can be understood as the source of exploding gradients when batch norm is applied to very deep networks, since very small changes in an input correspond to very large movements in network outputs. (b) The correlation between the two minibatches, as a function of layer, for the same network. Despite having a correlation near one at the input layer, the two minibatches rapidly decorrelate with depth. See ?? for a theoretical treatment.



Figure 3: Batch normalization strongly limits the maximum trainable depth. Colors show test accuracy for rectified linear networks with batch normalization and γ = 1, β = 0, = 10 -3 , N = 384, and η = 10 -5 B. (a) trained on MNIST for 10 epochs (b) trained with fixed batch size 1000 and batch statistics computed over sub batches of size B. (c) trained using RMSProp. (d) Trained on CIFAR10 for 50 epochs. In each case, White dashed line indicates a theoretical prediction of trainable depth, as discussed in the text.

Figure 4: Gradients in networks with batch normalization quickly achieve dynamical equilibrium. Plots of the relative magnitudes of (a) the weights (b) the gradients of the loss with respect to the pre-activations and (c) the gradients of the loss with respect to the weights for rectified linear networks of varying depths during the first 10 steps of training. Colors show step number from 0 (black) to 10 (green).

Figure 5: Three techniques for counteracting gradient explosion. Test accuracy on MNIST as a function of different hyperparameters along with theoretical predictions (white dashed line) for the maximum trainable depth. (a) tanh network changing the overall scale of the pre-activations, here γ → 0 corresponds to the linear regime. (b) Rectified linear network changing the mean of the preactivations, here β → ∞ corresponds to the linear regime. (c,d) tanh and rectified linear networks respectively as a function of , here we observe a well defined phase transition near ∼ 1. Note that in the case of rectified linear activations we use β = 2 so that the function is locally linear about 0. We also find initializing β and/or setting > 0 having positive effect on VGG19 with batchnorm. See Figs.8 and 9.

Figure7: relative gradient norms of different parameters in layer order (input to output from left to right), with γ and W interleaving. From dark to light blue, each curve is separated by (a) 3, (b) 5, or (c) 10 epochs. We see that after 10 epochs, the relative gradient norms of both γ and W for all layers become approximately equal despite gradient explosion initially.

Figure9: In the same setting as Fig.8, except we don't take the max over the unseen hyperparameter but rather set it to 0 (the default value).

we will also write B = B id . Definition C.3. Define the matrix G B = I -1 B 11 T . Let S G B be the space of PSD matrices of size B × B with zero mean across rows and columns, S G B := {Σ ∈ S B : G B ΣG B = Σ} = {Σ ∈ S B : Σ1 = 0}.

Theorem E.25. Let Σ be a positive-definite D × D matrix. Let φ : R → R. Suppose for any b > a > 0, 1. b a d(r 2 )|φ(z a /r)φ(z b /r)| exists and is finite for any z a , z b ∈ R (i.e. φ(z a /r)φ(z b /r) is locally integrable in r).

(t), . . .. Proof of Thm E.47. We have

respectively with the following domains Definition E.55. Denote by H B the space of symmetric matrices of dimension B. Also write H G B for the space of symmetric matrices Σ of dimension B such that GΣG = Σ (which is equivalent to saying rows of Σ sum up to 0). As in the case of S, we omit subscript B when it's clear from context. The following subspaces of H B and H G B will turn out to comprise the eigenspaces of G ⊗2 •T •G ⊗2 .

and λ T L,1 and λ T L,2 are the roots to the quadratic x 2 -(β 11 + β 22 )x + β 11 β 22β 12 β 21 .

10. Suppose φ is degree α positive-homogeneous. Then for Λ ∈ H B , E[D ⊗ (Dyy T ) + Dyy T ⊗ D : y ∼ N (0, Σ)] = α dV φ (ΠΣΠ) dΠ Π=I E[DΛyy T D + Dyy T ΛD : y ∼ N (0, Σ)] = α dV φ (Σ) dΣ {ΣΛ + ΛΣ} where D = Diag(φ (y)).

.1 SPHERICAL INTEGRATION AND GEGENBAUER EXPANSION Theorem H.3. Let Σ ∈ S kB . If Σ is CBSB1 then V B ⊕k φ ( Σ) is CBSB1 and equals Σ * := Σ * c * c * Σ * where Σ * is the BSB1 fixed point of Thm F.5 and c * ≥ 0 with √

1:B,B+1:2B = E[B φ (z 1:B ) ⊗ B φ (z B+1:2B ) : z ∼ N (0, Σ)] = E[φ • n(y 1:B ) ⊗ φ • n(y B+1:2B ) : y ∼ N (0, Σ G )]

( Σ) 1:B,B+1:2B = E[φ • n(x) : x ∼ N (0, (ab)G)] ⊗2By symmetry,E[φ • n(x) : x ∼ N (0, (ab)G)] = √ c * 1 E[φ • n(x) : x ∼ N (0, (ab)G)] ⊗2 = c * 11 T where √ c * := E[φ(n(x) 1 ) : x ∼ N (0, (ab)G)]. We can compute √ c * = E[φ(n(x) 1 ) : x ∼ N (0, G)] because n is scale-invariant = E[φ(n(ex) 1 ) : x ∼ N (0, I B-1 )]where e is as in Eq. (33)= E[φ((en(x)) 1 ) : x ∼ N (0, I B-1 )]because e is an isometry= E[φ( √ B(ev) 1 ) : x ∼ N (0, I B-1 )] with v = xθ) sin B-3 θby Lemma E.37. At the same time, by Proposition E.33, we can obtain the Gegenbauer expansion √ c(ê 1,: ) = a 0 .Proof. If Σ is full rank, we can show Fubini-Tonelli theorem is valid in the following computation by the same arguments of the proof of Lemma E.2.E[ y -2a z -2b f (y, z) : (y, z) ∼ N (0, Σ)] = E[f (y, z) ∞ 0 ds Γ(a) -1 s a-1 e -y 2 s ∞ 0dt Γ(b) -1 t b-1 e -y 2 t : (y, z) ∼ N (0a-1 t b-1 det(Σ(Σ -1 + 2D 2 )) -1/2 E f (y, z)

(Σ -1 + 2D 2 ) -1 = Σ(I A+B + 2D 2 Σ) -1 = Σ(D -1 + 2DΣ) -1 D -1 = ΣD(I A+B + 2DΣD) -1 D -1 = D -1 Ω(I A+B + 2Ω) -1 D -1 det(Σ(Σ -1 + 2D 2 )) = det(I A+B + 2ΣD 2 ) = det(I A+B + 2DΣD) = det(I A+B + 2Ω)The case of general Σ with rank Σ > 2(a + b) is given by the same continuity arguments as in Lemma E.2.Let Σ ∈ S 2B , Σ = Σ Ξ Ξ T Σ where Σ, Σ ∈ S B and Ξ ∈ R B×B . Consider the off-diagonal block of V B ⊕2 φ ( Σ). E[B φ (z) ⊗ B φ (z ) : (z, z ) ∼ N (0, Σ)] = E[φ(n(y)) ⊗ φ(n(y )) : (y, y ) ∼ N (0, Σ G )] = B α E[ y -α y -α φ(y) ⊗ φ(y ) : (y, y ) ∼ N (0, Σ G )] = B α Γ(α/2) ) α/2-1 det(I 2B + 2Ω) -1/2 E (y,y )∼N (0,Π) φ(y) ⊗ φ(y ) ) α/2-1 det(I 2B + 2Ω) -1/2 V φ (Π)6 (Rephrasing and adding to Corollary H.4). Let Σ ∈ S kB and φ be positive homoge-neous of degree α. If Σ is CBSB1 then V B ⊕k φ ( Σ) is CBSB1 and equals Σ * = Σ * c * 11 T c * 11 TΣ * where Σ * is the BSB1 fixed point of Thm F.5 andc * = c α B -b) 2 1 4π (B -1) α Γ((B -1)/2)Γ((α + 1)/2) Γ((α + B -1)/2)Definition H.7. Let T : H kB → H kB be a linear operator such that for any block matrix Σ ∈ H kB with k × k blocks Σ ij , i, j ∈ [k], of size B, we haveT { Σ} ii = U{ Σ ii }

Then for anyΛ ∈ R B×B satisfying GΛG = Λ, G ) 12 E (h,h )∼N (0, Σ G ) φ(n(h)) ⊗ φ(n(h )) Σ G =µ * G ⊕2. Combining Lemma F.18 with projection by e, we getJ{Λ} = 1 2 E y,y ∼N (0,µ * I B-1 )φ (n(ey)) ⊗ φ (n(ey )) µ * -2 yy Tµ * -1 I B-1 , e T

because of the independence of v from v and the spherical symmetry of the integrand. Thus this term drops out in G ⊗2 • J{Λ}.Next, we analyze H{Λ}= E v,v ∼S B-2 φ √ Bev ⊗ φ √ Bevv T e T Λev . This is a linear function in Λ with coefficientsH ab|cd = B -1 B E v∼S B-2 φ( √ B -1ê b,: v)(ê c,: v) E v∼S B-2 φ( √ B -1ê a,: v)(ê d,: v) .By the usual zonal spherical harmonics argument,

Definition H.11. Define the block L-shaped matrix L kB (Σ, Λ) :=

2. For each eigenvalue λ and eigenvectorΣ ∈ H G B of U, the matrix L kB (Σ, -(λw)

space generated by L 2B (G, -(λ ↑ Lλ ↑ M ) -1 V{G}) along with any matrix obtained from permutating its block columns and rows simultaneously. The associated eigenvalue is 0.3. The space generated by L2B (L, -(λ ↑ Lλ ↑ M ) -1 V{L}), where L = L B (B -2, 1), alongwith any matrix obtained from permutating its block columns and rows simultaneously. The associated eigenvalue is λ ↑ L . 4. The space generated byL 2B (M, -(λ ↑ Mλ ↑ M ) -1 V{M }), for any M ∈ M, along with any matrix obtained from permutating its block columns and rows simultaneously. The associated eigenvalue is λ ↑ M . Here, V is as in Eq. (57). The local convergence dynamics of Eq. (54) to the CBSB1 fixed point Σ * has the same eigenvalues λ ↑ M , λ ↑ L , λ ↑ M , and 0.Its dynamics is then given by the eigendecomposition of V B ⊕k φ ( Σ * ), in parallel to the single batch scenario.WLOG, we only need to consider the case k = 2. y)∼N (0, Σ * ) B φ (x)ΣB φ (x) T B φ (x)ΞB φ (y) T B φ (y)Ξ T B φ (x) T B φ (y)ΛB φ (y) T

y)∼N (0, Σ * ) [B φ (x) ⊗ B φ (y)] † {Ξ l+1 } = E (x,y)∼N (0, Σ * ) [B φ (y) T ⊗ B φ (x) T ]{Ξ l+1 } (60) = E (x,y)∼N (0, Σ * ) [B φ (y) T Ξ l+1 B φ (x) T ]Definition I.1. Let M B := R B×B be the set of all B × B real matrices, and letM G B := {Ξ ∈ M B : GΞG = Ξ}. Define T := E (x,y)∼N (0, Σ * ) [B φ (x)⊗B φ (y)].As in the single batch case, after one step of backprop,

since Σ * is CBSB1 with diagonal blocks Σ * (Corollary H.4). Then ξ is independent from η, so this is justE (x,y)∼N (0, Σ * ) [B φ (x) ⊗ B φ (y)] = E ξ∼N (0,µ * G) = ξ/ ξ . But notice that T := E ξ∼N (0,µ * G) Dr -1 (Ivv T ) is actually BSB1, with diagonal E rv∼N (0,µ * G)

G(Ivv T )Dδ where D = Diag(φ (n(y))), r = y , v = y/ y , and y ∼ N (0, GΣ * G). Because GΣ * G = λ * q⊗2 by Thm K.5, n(y) = √ Bv = ± √ B q with equal probability and r ∼ |N (0, λ * )| (the absolute value of a Gaussian). Furthermore, G(I-q⊗2 ) = G -q⊗2 because G q = q. Thus dB φ dΣ † Σ=Σ * δ = √ B|r| -1 (G -q⊗2 )(φ ( √ B sgn(r) q) δ), r ∼ N (0, λ * )Now notice the following Lemma K.6. G B -q⊗2 = 0 0 0

used x l βi to denote B φ (h l β ) i , δ l i = ∇ h l i L, and subscript i denotes slice over sample i in the batch, and subscript α denotes slice over neuron α in the layer. Then E

Peer review (iclr2019):
 This paper develops a mean field theory for batch normalization (BN) in fully-connected networks with randomly initialized weights. There are a number of interesting predictions made in this paper on the basis of this analysis. The main technical results of the paper are Theorems 5-8 which compute the statistics of the covariance of the activations and the gradients. Comments: 1. The observation that gradients explode in spite of BN is quite counter-intuitive. Can you give an intuitive explanation of why this occurs? 2. In a similar vein, there a number of highly technical results in the paper and it would be great if the authors provide an intuitive explanation of their theorems. 3. Can the statistics of activations be controlled using activation functions or operations which break the symmetry? For instance, are BSB1 fixed points good for training neural networks? 4. Mean field analysis, although it lends an insight into the statistics of the activations, needs to connected with empirical observations. For instance, when the authors observe that the structure of the fixed point is such that activations are of identical norm equally spread apart in terms of angle, this is quite far from practice. It would be good to mention this in the introduction or the conclusions.