A Mean Field Theory of Batch Normalization

Motivation

• The paper targets batch normalization, an important architectural component that had remained elusive to existing theoretical frameworks.Prior work analyzed initialization distributions and gradients in architectures including skip-connected, convolutional, and gated recurrent networks (Yang & Schoenholz, 2017; 2018; Xiao et al., 2018; Chen et al., 2018; Gilboa et al., 2019).Source: §1 · pp. 1–2

Theory and validation

• The authors develop a theory for randomly weighted and biased fully connected batch-normalized networks by reformulating random-network analysis for data batches and validating it with Monte Carlo simulations.Because batch-normalization statistics depend non-locally on the entire batch, the framework also predicts valid hyperparameter configurations.Source: §1 · p. 2

Gradient instability

• For most practical nonlinearities, batch normalization induces strong symmetry among batch-sample embeddings, causing gradients to explode exponentially with depth for every nonlinearity.This symmetry imposes strong limits on the maximum trainable depth of randomly initialized batch-normalized networks.Source: §1 · p. 2
• Making activations more linear at initialization reduces gradient explosion and partially increases trainable depth, but cannot eliminate the instability completely.The limitation arises even when the nonlinearity is changed, so tuning activation functions alone cannot fully resolve it.Source: §1 · p. 2

Learning dynamics

• After one gradient-descent step, batch-normalized networks naturally equilibrate gradient scales, enabling adequate training in shallower networks when initial gradients avoid numerical instability.Thus, the initially unfavorable gradient dynamics can settle into a regime with substantially more stable gradient scales.Source: §1 · p. 2

Bayesian inference

• In the wide-network limit, the induced prior over functions for batch-normalized networks can be computed exactly, enabling exact Bayesian inference.This extends prior exact-inference work for random neural-network function priors (Lee et al., 2017; de G. Matthews et al., 2018; Novak et al., 2019; Garriga-Alonso et al., 2019; Yang, 2019).Source: §1 · p. 2

Prior theoretical perspectives

• Earlier studies attributed batch normalization’s benefits to loss-landscape smoothing, favorable gradient dynamics, larger learning rates, or provable acceleration, while comparisons among normalization schemes remained inconclusive.These works include Santurkar et al. (2018), Bjorck et al. (2018), Luo et al. (2018), and Kohler et al. (2018).Source: §2 · p. 2

Scope and distinction

• Unlike most prior mathematical analyses, this work handles arbitrarily many batch-normalization layers, very general activation functions, and gradients flowing through batch variance.The paper argues that this deep stacking of batch normalization leads to gradient explosion.Source: §2 · pp. 2–3

Empirical validation

• The theory is an infinite-width initialization analysis whose resulting insights predict both training- and test-time behavior experimentally, strengthening earlier empirical observations by Philipp & Carbonell (2018).Philipp & Carbonell (2018) empirically noted gradient explosion in deep batch-normalized networks with various nonlinearities; this work provides a precise theoretical characterization of gradient statistics.Source: §2 · p. 3

3 THEORY

• Mean field theory replaces wide-network pre-activations with Gaussian variables, yielding an efficient covariance-matrix dynamical system whose fixed points govern deep signal propagation.The approximation is exact in the infinite-width limit, and covariance statistics can be computed using O(B^2) two-dimensional integrals.Source: §3 · p. 3
• Asymptotically deep vanilla networks converge to symmetric fixed points where all inputs have identical pre-activation norms and pairwise angles, a regime previously shown to be untrainable.Fixed-point stability is determined by Jacobian eigenvalues: all eigenvalues below one imply exponential convergence, whereas any eigenvalue above one implies instability.Source: §3 · p. 4

3.1 BATCH NORMALIZATION

• Batch normalization is modeled as a batchwise nonlinearity, with Laplace, Fourier, spherical-integration, and Gegenbauer methods producing analytic covariance recurrences across activation functions.The Gegenbauer method also establishes that gradient explosion occurs regardless of nonlinearity under general conditions.Source: §3.1 · pp. 5–6
• Batch-normalized covariance dynamics have permutation-symmetric BSB1 fixed points, but rapidly growing activations can cause winner-take-all behavior and symmetry-breaking BSB2 fixed points.Most practical nonlinearities, including ReLU, leaky ReLU, tanh, and sigmoid, are observed to converge to BSB1 fixed points under the stated assumption.Source: §3.1 · p. 6

3.1.1 LINEARIZED DYNAMICS

• Linearized batch-normalization dynamics decompose into three orthogonal eigenspaces capturing average norms, norm fluctuations, and inter-sample covariances.The average-norm component vanishes after one step, while increasing batch size makes deep ReLU-batchnorm networks require more layers to reach a BSB1 fixed point; more-linear activations reduce the dominant instability.Source: §3.1.1 · pp. 7–8

3.1.2 GRADIENT BACKPROPAGATION

• For batch-normalized fully connected networks whose covariances converge to BSB1 fixed points, gradient norms squared grow as β^L+o(L) for every permitted nonlinearity under gradient independence.The backpropagation recurrence includes an additional sum over batch elements because batch-normalization activations are not pointwise.Source: §3.1.2 · p. 8
• Batch normalization forces gradient signals and layer-Jacobian singular values away from 1, causing exponential gradient explosion that cannot be removed by tuning initialization variances or nonlinearities.The eigendecomposition of the linearized gradient dynamics establishes the explosion; this contrasts with vanilla networks, where variances can control mean-field dynamics.Source: §3.1.2 · p. 9
• For ReLU batch-normalized networks, gradient norms explode exponentially, with the rate decreasing to π/(π−1) ≈1.467 as batch size B grows; linear networks explode at rate B^-2.Figure 1 confirms the theoretical eigenvalue and depth predictions using Monte Carlo simulations.Source: §3.1.2 · p. 9, Figure 1 · p. 10
• Weight-gradient norms explode at the same rate as hidden-preactivation gradients, because they are determined by the corresponding gradient covariance.The weight gradient at layer l is expressed through the trace pairing ⟨Π_l, μ*_G⟩ = μ* tr Π_l.Source: §3.1.2 · p. 9
• Increasing epsilon can ameliorate but not eliminate gradient explosion, while the linear activation maximally suppresses explosion in the ideal infinite-batch limit.The epsilon result concerns values somewhat larger than those typically used; the linear-activation conclusion is stated alongside cross-batch information-loss results.Source: §3.1.2 · p. 9, §3.2 · p. 11

3.2 CROSS-BATCH DYNAMICS

• For ReLU and identity activations, the cross-batch covariance dynamics converge to a BSB1 fixed point whose off-diagonal blocks share a common constant.After mean-centering, the limiting covariance is a multiple of the identity, indicating loss of mutual information between batches.Source: §3.2 · pp. 10–11
• Batch normalization causes cross-batch covariance and cross-batch gradient covariance to decay at the same rate, so similar input batches decorrelate exponentially with depth.The decay rate is always below 1, is maximized by the identity activation for fixed B, and approaches 1 from below as B increases.Source: §3.2 · p. 11

• The depth scale ξ = 1/log λ↓ closely agrees with the maximum trainable depth across experiments, linking trainability to gradient explosion.Here, λ↓ is associated with the largest eigenvalue G of the Jacobian through the theory’s gradient-growth analysis.Source: §4 · p. 11
• Batch normalization strongly limits trainable depth: networks deeper than about 50 layers remain untrainable across batch sizes, optimizers, and MNIST/CIFAR10 settings.The observed trainable regions are nearly identical across the tested configurations, matching theoretical predictions overlaid in Fig. 3.Source: §4 · pp. 11–12, Figure 3 · p. 12
• After one training step, exploding initialization gradients drive exponential weight growth, then batch normalization rescales gradients into a stable equilibrium that can vanish exponentially in depth.In sufficiently deep networks, lower-layer gradients become nearly zero, causing those layers to remain mostly unchanged and impairing training.Source: §4 · p. 12, Figure 4 · p. 12
• Trainable depth increases when nonlinearities approach the linear regime, while sufficiently large ϵ suppresses gradient explosion and makes very deep networks trainable.For tanh, this occurs as γ decreases toward 0; for rectified linear activations, as β increases; both trends and the ϵ transition agree with theory.Source: §4 · pp. 12–13, Figure 5 · p. 13

5 CONCLUSION

• Batch normalization unavoidably causes gradients to explode with depth at initialization.The work characterizes this as a counterintuitive property of batch-normalized networks.Source: §5 · p. 13
• Several methods reduce gradient explosion, enabling training of significantly deeper networks with batch normalization.Source: §5 · p. 13
• The theory provides a foundation for future research on more advanced network architectures and topologies.Source: §5 · p. 13

Gradient equilibration

• After 10 epochs, relative gradient norms for W and γ become approximately equal across all layers despite initial gradient explosion.The norms ∥∇θL∥/∥θ∥ equilibrate to about the same magnitude for weight and BN-scale parameters.Source: A VGG19 WITH BATCHNORM ON CIFAR100 · p. 16, Figure 7 · p. 17
• The experiment tracks relative gradient norms across layers and parameter types at 3, 5, and 10 epochs.γ and W parameters are interleaved in layer order from input to output.Source: Figure 7 · p. 17

Training acceleration

• Setting ϵ > 0 and/or initializing β > 0 produces acceleration effects, especially during initial training.The study sweeps learning rate, β initialization, and ϵ across eight random seeds per combination, aggregating median training and validation accuracy.Source: A VGG19 WITH BATCHNORM ON CIFAR100 · p. 17, Figure 8 · p. 17

B GRADIENT INDEPENDENCE ASSUMPTION

• During backpropagation, every multiplication by W^T is replaced by multiplication with an iid copy of W, following prior mean-field analyses.This is the paper’s gradient independence assumption, attributed to Schoenholz et al. (2016), Yang & Schoenholz (2017), and Xiao et al. (2018).Source: B GRADIENT INDEPENDENCE ASSUMPTION · p. 17
• Computations under the gradient independence assumption agree closely with simulations, as shown in Fig. 1.Yang (2019) rigorously justified the assumption for moment computations in several architectures without batch normalization.Source: B GRADIENT INDEPENDENCE ASSUMPTION · p. 18
• The existing proof does not extend directly to batch normalization because its Jacobian is singular at zero, although the authors expect an extension to be possible.The anticipated extension is based on Yang & Schoenholz (2017)’s work.Source: B GRADIENT INDEPENDENCE ASSUMPTION · p. 18

C NOTATIONS

• For measurable Φ:R^B→R^A, the transform V_Φ maps a PSD covariance matrix Σ to E[Φ(h)⊗2] for h∼N(0,Σ).For scalar φ, V_φ denotes the same transform applied coordinatewise.Source: C NOTATIONS · p. 18
• The notation B_φ represents batch normalization followed by coordinatewise application of φ to a batch of neuronal activations.The paper also writes B=B_id when the activation is the identity.Source: C NOTATIONS · p. 18
• G denotes the projection that removes a batch’s sample mean, while S^G_B contains PSD matrices with zero row and column means.The zero-mean covariance space is characterized by GΣG=Σ, equivalently Σ1=0.Source: C NOTATIONS · pp. 18–19
• The paper defines Gamma, Pochhammer, and Beta functions, uses dot and trace inner products, and treats vectors as columns in its Jacobian convention.These conventions support the subsequent matrix and differential calculations.Source: C NOTATIONS · p. 19

E TECHNIQUES

• The appendix compares Laplace, Fourier, spherical-integration, and Gegenbauer methods, with ultrasymmetric-operator eigendecompositions supporting the dynamics’ asymptotic analysis.Laplace yields closed forms for positive-homogeneous nonlinearities; spherical integration supports general measurable nonlinearities numerically; Gegenbauer methods expose eigenvalue structure but may require many integrals.Source: D A GUIDE TO THE REST OF THE APPENDIX · pp. 19–20

D A GUIDE TO THE REST OF THE APPENDIX

• The appendix analyzes forward, backward, and cross-batch dynamics, then extends the results beyond BSB1 fixed points and computes weight-gradient moments.Corollary F.3 proves global convergence for φ = id, while other nonlinearities receive local analyses; gradient explosion and cross-batch-gradient decorrelation are also characterized.Source: D A GUIDE TO THE REST OF THE APPENDIX · p. 20

E.1 LAPLACE METHOD

• Schwinger parametrization produces the Laplace master equation, which is well-defined and continuous for rank Σ > 2k and simplifies positive-homogeneous transforms.The master equation rewrites Gaussian expectations involving ∥y∥^-2k as an integral over s, while matrix simplifications exploit the projection structure of G.Source: §E.1 · pp. 21–23

E.1.1 α-RELU

• For Gaussian settings, the cited results give J2(θ) = 3 sin θ cos θ + (π − θ)(1 + 2 cos2 θ), while for α ∈ [1/2, 1), Jα(c) = c has two solutions.The fixed points are unstable at 1 and stable at c∗ ∈ (0, 1).Source: §E.1.1 · p. 24

E.1.2 POSITIVE-HOMOGENEOUS FUNCTIONS IN 1 DIMENSION

• Every one-dimensional degree-α positive-homogeneous function with φ(0) = 0 equals aρα(x) − bρα(−x), enabling its V and W transforms to be expressed through α-ReLU transforms.The resulting kernel is Jφ(c) = (a^2 + b^2)Jα(c) − 2abJα(−c), with Jφ(1) = a^2 + b^2 and Jφ(0) = (a − b)^2Jα(0).Source: §E.1.2 · pp. 25–26
• For φ(c) = aρα(c) − bρα(−c) with α > 1/2 on R \ {0}, Proposition E.23 derives identities connecting φ′ and J_α-related expressions.The proof uses φ′(c) = aαρα−1(c) + bαρα−1(−c) and Proposition E.14 to establish the first identity, followed by two further equations.Source: §E.1.2 · pp. 26–27
• For any Σ ∈SB and diagonal matrix D, Proposition E.24 states that diagonal matrices can be factored out of V_φ when φ is positive-homogeneous of degree α.This generalizes the subsection’s preceding observation that diagonal matrices can be factored out of V_φ.Source: §E.1.2 · p. 27

E.2 FOURIER METHOD

• The Fourier derivation is formally invalid at intermediate steps, but Theorem E.25 rigorously justifies its final identity using mollifiers under explicit integrability assumptions satisfied by common nonlinearities.The problematic steps involve Fourier-expanding a delta function and swapping a non-absolutely-integrable integration order; mollifier limits and dominated convergence repair the argument.Source: §E.2 · pp. 28–30

E.3 SPHERICAL INTEGRATION

• Spherical coordinates express batch-normalization integrals as expectations over the uniform sphere, while symmetry and coordinate-change lemmas reduce many high-dimensional integrals to only a few dimensions.The reductions use sparse realizations of the centering basis and integrate out irrelevant angular coordinates when the integrand depends on only a small subset.Source: §E.3 · pp. 33, 35, 36

E.4 GEGENBAUER EXPANSION

• A new Gegenbauer identity diagonalizes a derivative-dependent quadratic form expressed through a spherical integral, enabling the proof of necessary gradient explosion under the BSB1 fixed-point assumption.The construction uses the relationship between Gegenbauer polynomials and zonal spherical harmonics, including their reproducing property.Source: §E.4 · p. 37
• Theorem E.47 shows that the Dirichlet-like form of φ diagonalizes over the Gegenbauer basis.The proof expands the relevant expressions in Gegenbauer polynomials and establishes cancellation of higher-order coefficients.Source: §E.4 · pp. 38–39

E.5 SYMMETRY AND ULTRASYMMETRY

• Common nonlinearities converge to BSB1 fixed points, whose covariance matrices have one diagonal and one off-diagonal value and are isotropic in a codimension 1 subspace.A BSB1(a, b) matrix has eigenvalues a − b and (B − 1)b + a, while G BSB1(a, b) G = (a − b)G.Source: §E.5 · p. 40
• Under ultrasymmetry, T decomposes into BSB1, L, and M invariant eigenspaces, including an M eigenspace of dimension B(B − 3)/2 with eigenvalue λT M = [12|12] + [12|21] − 4[12|13] + 2[12|34].The BSB1 and L components each yield two eigenspaces, while simultaneous row-column permutations preserve the corresponding eigenvectors and eigenvalues.Source: §E.5 · pp. 42, 44
• Although these eigenspaces are generally nonorthogonal under the trace inner product, projection by G⊗2 makes the operator self-adjoint.Thus T is typically not self-adjoint before projection, but gains self-adjointness after the batch-normalization projection.Source: §E.5 · p. 44
• Projected ultrasymmetric operators decompose into orthogonal eigenspaces RG and L with explicit eigenvalues determined by operator coefficients.For DOS operators, the eigenspaces and eigenvalues simplify through direct two-dimensional linear maps.Source: §E.5 · pp. 44–46

F.1 GLOBAL CONVERGENCE

• For φ = id, the forward covariance dynamics converge exponentially to a rank-C limit whose nonzero eigenvalues become equal, with gap O((1 − 2 C+2)^l).The result applies to rank-C initial covariances and identifies the limiting matrix in the initial eigenspace basis.Source: §F.1 · p. 50
• The same exponential convergence and asymptotic gap bound extend to projected dynamics when GΣG has rank C < B.For general nonlinearities, the basin of attraction is not characterized and local fixed-point analysis is used instead.Source: §F.1 · p. 50

F.2 LIMIT POINTS

• Permutation-equivariance motivates studying fixed points of the form μI + ν11^T, and symmetry ensures the covariance map preserves this form.For any such input covariance, the mapped covariance is the same, yielding the unique fixed point within this class.Source: §F.2 · p. 50, §F.2.1 · p. 50

F.2.1 SPHERICAL INTEGRATION

• The unique BSB1 fixed point has diagonal and off-diagonal entries expressed through one- and two-dimensional integrals, enabling numerical computation.Theorem F.5 establishes uniqueness and supplies the integral characterization.Source: §F.2.1 · pp. 50–51

F.2.2 LAPLACE METHOD

• For degree-α positive-homogeneous φ, VBφ maps BSB1 matrices to BSB1 matrices, yielding a closed-form fixed point that exists and is unique.The diagonal and off-diagonal entries are expressed using Kα,B, Jφ, and cα.Source: §F.2.2 · p. 51
• At the fixed point, GΣ*G is proportional to the identity on the B−1-dimensional centered subspace, with proportionality μ*.This follows from the fixed-point characterization for BSB1 covariance matrices.Source: §F.2.2 · p. 52
• For any BSB1 covariance, the relative batch-normalization map remains BSB1; for φ(x)=relu(x)−relu(−x), its centered second-order map scales G by B/(B−1).The latter specializes the positive-homogeneous result to α=1 and φ=id.Source: §F.2.2 · p. 52

F.2.3 GEGENBAUER EXPANSION

• The BSB1 fixed point also admits a clean description through the Gegenbauer coefficients of φ.The gradient-explosion rate later decomposes similarly, whereas forward-convergence eigenvalues do not generally share this decomposition.Source: §F.2.3 · pp. 52–53

F.3 LOCAL CONVERGENCE

• Local convergence is studied by linearizing Eq. (43) on PSD operators and analyzing the Jacobian at the unique BSB1 fixed point.Batch normalization is decomposed as mean-centering, standard-deviation normalization, and coordinate-wise nonlinearity.Source: §F.3 · p. 53
• The auxiliary operator U captures the relevant Jacobian structure, and its nonzero eigenvalues determine those of the related composed operators through a general AB–BA eigenspace correspondence.The operators can be diagonalized, and the Jacobian spectrum can be recovered from U.Source: §F.3 · p. 54

F.3.1 SPHERICAL INTEGRATION

• Spherical integration reduces the eigenvalue calculation to ultrasymmetric expectations Wij|kl involving φ(n(ex)), coordinates of ex, and Gaussian x.Ultrasymmetry permits restricting indices to {1,2,3,4}, while further symmetry and trigonometric expansion integrate out θ3 and θ4.Source: §F.3.1 · p. 56

F.3.2 GEGENBAUER EXPANSION

• The Gegenbauer expansion expresses the local convergence eigenvalue on L as a ratio of quadratic forms in the coefficients {a_l}ᵢ.The analysis is obtained by differentiating the fixed-point relation along paths in the batch-normalized covariance space.Source: §F.3.2 · pp. 57, 59
• The eigenvalue expressions are independent of a_0 because batch normalization is invariant under additive shifts of φ.The corresponding Gegenbauer coefficients v_{B−1,0}, w_{B−1,0}, and u_{B−1,0} vanish.Source: §F.3.2 · p. 59
• The L eigenvalue is typically ≥1, while more explosive φ makes convergence to a BSB1 fixed point less likely.Higher-degree Gegenbauer components can increase endpoint explosiveness, consistent with slower or failed convergence.Source: §F.3.2 · p. 59
• For B ≥10, Fig. 10 suggests w_{B−1,l}−v_{B−1,l} ≥0 for every l ≥1, supporting the eigenvalue analysis.The comparison is based on plots across B for several l values.Source: §F.3.2 · p. 59, Figure 10 · p. 60
• The M eigenvalue is typically much smaller than the L eigenvalue, although counterexamples such as sin exist.The M expansion is derived using analogous Gegenbauer and symmetry calculations.Source: §F.3.2 · pp. 60–62

F.3.3 LAPLACE METHOD

• The Laplace method derives eigendecompositions for G⊗2 ◦ dVφ and identifies the operator’s eigenvalue structure for positive-homogeneous activations.The analysis includes general Gaussian-moment conditions and specialized positive-homogeneous cases.Source: §F.3.3 · p. 64
• The Laplace analysis also obtains an asymptotic characterization of the relevant operator as B → ∞.This result is stated as Proposition F.34 after routine computation.Source: §F.3.3 · p. 67
• For α ≥ 1, the maximum eigenvalue of U is attained by eigenspace L for sufficiently large B.For α ∈ (−1/2, ∞), the optimized comparison exceeds α precisely when α < 1.Source: §F.3.3 · p. 67

G BACKWARD DYNAMICS

• Backward dynamics are formulated through the V operator and analyzed by linearizing around the BSB1 fixed point, whose large-depth behavior is governed by an eigendecomposition.The adjoint G⊗2 ◦ V_B′φ(Σ∗) is studied because it has the same eigenvalues and is self-adjoint.Source: G BACKWARD DYNAMICS · p. 68

G.1.1 SPHERICAL INTEGRATION

• At the BSB1 fixed point, the derivative operator F(Σ∗) is ultrasymmetric, enabling its eigenvalues to be obtained from structured spherical expectations.The spherical method reduces entries to two-dimensional integrals after integrating out θ3 and θ4, then applies Theorem E.62.Source: §G.1.1 · pp. 69–70

G.1.2 GEGENBAUER EXPANSION

• The Gegenbauer expansion simplifies F(Σ∗){G} because G acts as the identity on v, with separate formulas for diagonal and off-diagonal entries.The diagonal calculation lifts a spherical integral to dimension B, whereas the off-diagonal calculation uses Theorem E.47 instead.Source: §G.1.2 · pp. 70–72
• The Gegenbauer calculation gives the eigenvalue associated with eigenspace R_G, and this eigenvalue is minimized over φ if and only if φ is linear.The statement assumes Σ∗ is the unique BSB1 fixed point.Source: §G.1.2 · pp. 72–73

G.2 LAPLACE METHOD

• The Laplace method rewrites the relevant operator for degree-α positive-homogeneous activations and exploits shared ultrasymmetric eigenspaces R_G, M, and L.The decomposition uses D = Diag(φ′(y)), with D degree α − 1 homogeneous, and evaluates terms through V_φ-related identities.Source: §G.2 · pp. 73–74

3. Eigenspace M with eigenvalue

• The analysis computes the eigenvalues associated with eigenspaces R_G, L, and M, then combines them to obtain the eigendecomposition of B.The final result follows from Eq. (53) and Theorems G.6, G.9, and G.11.Source: §3 · pp. 77–78

H CROSS BATCH: FORWARD DYNAMICS

• For multiple batches, the extended covariance dynamics act on kB-dimensional covariances while restricting any diagonal B × B block recovers the single-batch dynamics.The construction uses direct-sum operator notation and applies batch normalization separately to each batch block.Source: H CROSS BATCH: FORWARD DYNAMICS · p. 78

H.1 LIMIT POINTS

• The cross-batch analysis searches for fixed points with CBSB1 structure: a common BSB1 diagonal block and a constant off-diagonal block.Global convergence is not established; the study instead examines fixed points and local behavior around them.Source: §H.1 · p. 78

H.1.1 SPHERICAL INTEGRATION AND GEGENBAUER EXPANSION

• The proposed cross-batch fixed point has diagonal blocks equal to the single-batch fixed point Σ* and a nonnegative off-diagonal parameter c*.The diagonal claim follows because restricting to diagonal blocks recovers Eq. (43) and invokes Theorem F.5.Source: §H.1.1 · p. 79
• The off-diagonal calculation uses Gaussian independence, isometry, spherical integration, and a Gegenbauer expansion to characterize the cross-batch coupling.The proof is presented for k = 2, with the general k case stated to follow by the same reasoning.Source: §H.1.1 · p. 79
• For φ(x) = aρ_α(x) − bρ_α(−x), the resulting off-diagonal parameter c* is given explicitly using gamma functions and depends on a − b.The expression is obtained by expanding the beta function and combining it with Theorem H.3.Source: §H.1.1 · p. 80

H.1.2 LAPLACE METHOD

• The Laplace cross-batch master equation makes the weighted Gaussian expectation well-defined and continuous whenever rank Σ > 2(a + b), with the full-rank proof extended by continuity.The derivation uses Fubini–Tonelli for full-rank Σ and recovers the general-rank case through continuity.Source: §H.1.2 · pp. 80–81

H.2 LOCAL CONVERGENCE

• The linearized dynamics around CBSB1 fixed points are organized as blockwise diagonal-off-diagonal semidirect operators, BDOS(U, V, W), separating diagonal and off-diagonal block effects.The off-diagonal components of Λ do not affect the diagonal blocks, and symmetry contributes the factor 2 in W.Source: §H.2 · pp. 82–83

H.2.1 SPHERICAL INTEGRATION AND GEGENBAUER EXPANSION √

• Spherical integration and Gegenbauer expansions exploit Gaussian independence and spherical symmetry to simplify the operator acting on cross-batch perturbations.A symmetry argument removes one term from G⊗2 ◦ J{Λ}, enabling the subsequent eigendecomposition analysis.Source: §H.2.1 · pp. 83–84
• For nonsingular V − wI_B, BDOS(V, U, w) has eigenspaces consisting of zero-diagonal block matrices with eigenvalue w and lifted eigenvectors of U with their corresponding eigenvalues.The lifted eigenvectors are formed from Σ and (λ − w)^−1V{Σ}, together with simultaneous block-row and block-column permutations.Source: §H.2.1 · p. 85
• Specializing the BDOS eigendecomposition to ˜F gives the local convergence dynamics toward the CBSB1 fixed point, with the same framework extending from 2 batches to k batches.The eigenvalue associated with eM is also the eigenvalue governing cross-batch correlations in gradient propagation.Source: §H.2.1 · pp. 85–86

H.2.2 LAPLACE METHOD

• For positive homogeneous activations, the new eM eigenvalue admits a direct expression in terms of the activation’s J function via the Laplace method.The result applies to degree α positive homogeneous φ and is derived around the CBSB1 fixed point using a smooth perturbation path.Source: §H.2.2 · p. 86
• Along the perturbation path, (I + 2Ω_0)^−1 ˙Ω_0 has zero diagonal blocks and zero trace, so only the second term in the governing sum remains.This cancellation isolates the contribution used to derive the eM eigenvalue at the CBSB1 fixed point.Source: §H.2.2 · pp. 86–87

Off-diagonal eigendecomposition

• Only the eigenspace {x : Gx = x}⊗2 = MG has a nonzero eigenvalue after composition with G⊗2, while the other eigenspace vanishes.The underlying operator has two eigenspaces, but composition with G⊗2 preserves only the former.Source: I CROSS BATCH: BACKWARD DYNAMICS · pp. 88–89

Positive-homogeneous activations

• For positive-homogeneous activations, the off-diagonal operator is explicitly characterized through the activation’s degree and coefficients, with an equivalent Gegenbauer-coefficient expression.The paper derives the eigenvalue for φ(c) = aρα(c) − bρα(−c) and then rewrites it using Gegenbauer coefficients.Source: I CROSS BATCH: BACKWARD DYNAMICS · pp. 89–90

Cross-batch backward dynamics

• The backward cross-batch dynamics is separable across B × B blocks, with diagonal blocks following Eq. (51) and off-diagonal blocks evolving by a scalar eigenvalue.The operator acts independently on blocks, and its off-diagonal dynamics is multiplication by λ↓.Source: I CROSS BATCH: BACKWARD DYNAMICS · p. 91

Activation choice

• Among fixed-batch-size activations, φ = id maximizes the backward eigenvalue, with λ↓eM = 1 and convergence to this limit from below.Thus the linear activation maximally mitigates the exponential loss of cross-batch information.Source: I CROSS BATCH: BACKWARD DYNAMICS · p. 91

Implications for deep networks

• For any activation inducing convergence to a CBSB1 fixed point, cross-batch correlation decreases exponentially with depth in both forward and backward propagation.The result identifies this as chaotic behavior caused by stacking batch normalization.Source: I CROSS BATCH: BACKWARD DYNAMICS · p. 91

J UNCOMMON REGIMES

• For B = 2, batch normalization gradients vanish immediately, whereas B ≥ 4 produces gradient explosion; B = 3 empirically appears similarly explosive.B = 1 is not well-defined, and the authors conjecture that Theorem 3.9 extends to B = 3.Source: J UNCOMMON REGIMES · p. 91
• When the fixed point is not BSB1, experiments indicate convergence to BSB2, whose dominant block can emerge in different diagonal locations through symmetry breaking.Finite-network weight randomness breaks the symmetry lost in the mean-field limit, while small input fluctuations can select among BSB1 fixed points.Source: J UNCOMMON REGIMES · p. 91
• At BSB2 fixed points, the dominant block is expected to follow smaller-batch BSB1-like dynamics and therefore exhibit similar gradient explosion.A specific BSB2 form with a 1 × 1 dominant block is analyzed further in Appendix K.1.Source: J UNCOMMON REGIMES · p. 91, K TOWARD UNDERSTANDING BSB2 FIXED POINTS · p. 93
• Finite width can change the favored fixed point: for B = 10, width 100 favors BSB2, whereas widths 1000 and above favor BSB1.This behavior is attributed to finite-width effects for certain nonlinearities.Source: J UNCOMMON REGIMES · p. 92

K TOWARD UNDERSTANDING BSB2 FIXED POINTS

• No rigorous BSB2 result is currently established, and the forward and backward dynamics are left for future work.The authors nevertheless expect the main block to undergo BSB1-like dynamics and produce similar gradient explosion.Source: K TOWARD UNDERSTANDING BSB2 FIXED POINTS · pp. 92–93
• The BSB2 analysis characterizes these matrices through eigenspaces and eigenvalues, including block-specific zero-sum subspaces and a one-dimensional zero-eigenvalue space.Theorem K.2 gives the eigendecomposition for general BSB2 matrices, and the section specializes it to the structured BSB2 form.Source: K TOWARD UNDERSTANDING BSB2 FIXED POINTS · p. 92
• For fixed B′ < B, forward BSB2 dynamics reduce to a generally nondegenerate two-dimensional system on the Z1 and Rq eigenvalues.The Gegenbauer method does not immediately generalize because these eigenspaces generally have different eigenvalues, preventing a spherical-integral reduction.Source: K TOWARD UNDERSTANDING BSB2 FIXED POINTS · p. 93

K.1 \ BSB2

• For rapidly growing activations such as φ(x) = relu(x)30 and not-too-large width, the analysis sometimes observes pathological gradient vanishing.These situations resemble gradient vanishing for B = 2 batch normalization.Source: §K.1 · p. 93
• The BSB2 fixed point is rank 1, with B(d, c, b)G equal to λq̂⊗2 for λ = B^-1.The analysis extends Corollary K.3 to B′ = 1 and establishes a unique fixed point.Source: §K.1 · p. 93
• At the rank-1 fixed point, backpropagation reduces the normalized activation structure to q̂ with equal probability and radial scale r distributed as |N(0, λ*)|.This follows from GΣ*G = λ*q̂⊗2 and explains the singular |r|^-1 factor in the gradient expression.Source: §K.1 · p. 94
• At the BSB2 fixed point, gradient behavior is dichotomous: nonconstant δ produces severe explosion, whereas δ2 = 0 yields deterministic vanishing.The explosion arises from a |r|^-1 factor, while the vanishing case is exactly zero.Source: §K.1 · p. 94
• In PyTorch Float and Double tensors, gradients empirically vanish rapidly for generic δ with φ(x) = relu(x)30 near the BSB2 fixed point.The passage attributes this to numerical issues because finite-precision computation does not completely eliminate δ1.Source: §K.1 · p. 94

L WEIGHT GRADIENTS

• The weight-gradient analysis assumes large depth and well-behaved φ so that Σl−1 approximately equals the unique BSB1 fixed point of Eq. (43).It uses the gradient-independence assumption together with the weight-gradient backpropagation equation.Source: L WEIGHT GRADIENTS · pp. 94–95
• In general, the magnitude of weight gradients is expected to follow the magnitude of the hidden-unit gradient δl.This links weight-gradient scale directly to hidden-unit gradient propagation under the stated fixed-point approximation.Source: L WEIGHT GRADIENTS · p. 95