arXiv CSopen access

Inter-Layer Hessian Analysis of Neural Networks with DAG Architectures

Unknown · 2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

machine learning, deep learning, neural networks

Inter-Layer Hessian Analysis of Neural Networks with DAG Architectures Maxim Bolshima,∗, Alexander Kugaevskikha a ITMO University, Saint Petersburg, Russia

arXiv:2604.11639v1 [cs.LG] 13 Apr 2026

Abstract Modern automatic differentiation frameworks (JAX, PyTorch) return the Hessian of the loss function as a monolithic tensor, without exposing the internal structure of inter-layer interactions. This paper presents an analytical formalism that explicitly decomposes the full Hessian into blocks indexed by the DAG of an arbitrary architecture. The canonical decomposition H = H GN + H T separates the Gauss–Newton component (convex part) from the tensor component (residual curvature responsible for saddle points). T For piecewise-linear activations (ReLU), the tensor component of the input Hessian vanishes (Hv,w ≡ 0 a.e., f GN Hv,w = Hv,w ⪰ 0); the full parametric Hessian contains residual terms that do not reduce to the GGN. Building on this decomposition, we introduce diagnostic metrics (inter-layer resonance R, geometric coupling C, stable rank D, GN-Gap) that are estimated stochastically in O(P ) time and reveal structural curvature interactions between layers. The theoretical analysis explains exponential decay of resonance in vanilla networks and its preservation under skip connections; empirical validation spans fully connected MLPs (Exp. 1–5) and convolutional architectures (ResNet-18, ∼11M parameters, Exp. 6). When the architecture reduces to a single node, all definitions collapse to the standard Hessian ∇2θ L(θ) ∈ Rp×p . Keywords: Hessian analysis, DAG architecture, Gauss–Newton decomposition, curvature diagnostics, inter-layer interaction, neural network

1. Introduction The second-order Hessian ∇2 L plays a fundamental role in curvature analysis of the loss landscape and in the design of optimization methods for neural networks (Martens, 2020). Second-order methods such as Newton, trust-region, and their variants require accurate curvature information for efficient optimization (Nocedal and Wright, 2006). However, for deep neural networks computing and storing the full Hessian is computationally prohibitive, requiring various approximations. The most common approach—the Gauss–Newton (GN) approximation—captures only a subset of the second derivatives, disregarding significant curvature contributions (Schraudolph, 2002; Martens, 2010). Block-diagonal approximations (K-FAC; Martens and Grosse, 2015) are restricted to diagonal blocks and do not capture inter-layer interactions in arbitrary DAG architectures. Hessian–vector product (HVP) methods (Pearlmutter, 1994) yield exact Hessian–vector products but do not reveal the structure of interlayer curvature. Moreover, the GN approximation is inherently blind to negative curvature: the block matrix GN [Hv,w ] ⪰ 0, rendering it insensitive to saddle points (Dauphin et al., 2014). The full Hessian H = H GN +H T recovers this information through the tensor component H T . Non-smooth activations (ReLU, max-pooling) additionally require a generalized second-derivative apparatus that has not been developed for full inter-layer blocks in the existing literature. ∗ Corresponding author

Email addresses: [email protected] (Maxim Bolshim), [email protected] (Alexander Kugaevskikh) URL: https://orcid.org/0009-0009-3037-9852 (Maxim Bolshim), https://orcid.org/0000-0002-6676-0518 (Alexander Kugaevskikh)

Preprint submitted to Neural Networks

April 14, 2026

Recent studies of per-layer Hessian spectra (Bolshim and Kugaevskikh, 2025; Ghorbani et al., 2019) suggest that diagonal blocks carry diagnostic information about overfitting and expressivity that is invisible in the global Hessian; however, existing analyses remain confined to ∂ 2 L/∂θv2 and do not address inter-layer interactions. The present work addresses this gap: we construct an analytical formalism for structural analysis of interlayer Hessian blocks ∂ 2 L/∂fv ∂fw and ∂ 2 L/∂θv ∂θw in neural networks with arbitrary DAG architectures. Standard APIs of modern frameworks (JAX, PyTorch) return the Hessian as a monolithic tensor, without exposing its internal block structure. We propose an explicit analytical decomposition: recurrence relations over the graph, the canonical splitting H GN + H T , path-wise contribution analysis, and diagnostic metrics (GN-Gap, resonance, coupling). This enables one to prove structural theorems about architectures (effects of skip connections, rank constraints, curvature decay) that do not follow from the computation of individual Hessian components. Key distinction from the GGN literature: our contribution is not a more efficient computation of the GN approximation (K-FAC and its variants already achieve this), but rather an explicit analytical decomposition of the full Hessian H = H GN + H T at the level of arbitrary DAG blocks, together with structural theorems that cannot be derived from H GN alone. Contributions. f C1. Full inter-layer Hessian for DAG architectures. A recurrence relation for blocks Hv,w accounting for all pure, mixed, and cross-block second derivatives with respect to inputs and parameters, including weight sharing across nodes (Eq. equation 2, Eq. equation 5, Appendix F).

C2. Canonical decomposition H = H GN + H T . Splitting into the Gauss–Newton component (block GN ] ⪰ 0) and the tensor component, formally explaining the blindness of the GN approxmatrix [Hv,w imation to saddle points; the quantitative deviation measure is the GN-Gap (Theorem 7, Definition 11). C3. Non-smooth case: input Hessian for ReLU networks. A rigorous proof that for piecewiseT ≡ 0 a.e., and linear activations (ReLU) the tensor component of the input Hessian vanishes: Hv,w f GN Hv,w = Hv,w in activation space. The parametric residual equation 5 is determined by mixed input–parameter derivatives that do not involve σ ′′ ; this residual does not vanish in general. The GN-Gap metric quantitatively measures the deviation of H f from H GN ; correctness of AD computations is established via conservative set-valued fields (CSVF) theory (Bolte and Pauwels, 2021) (Appendix D). C4. Curvature routing and the HVP operator. Recurrence equation 2 formalizes curvature routing through the DAG: decomposition of H GN +H T along edges and path-wise contributions. Collapsing the graph recursion over all nodes analytically reduces it to an O(P ) HVP operator whose cost matches that of a backward pass (Definition 20, Corollary 21, Appendix I). Building on this framework, we define diagnostic metrics (inter-layer resonance R, geometric coupling C, stable rank D, GN-Gap) that reveal structural curvature interactions between layers; details are in Appendix G. The metrics are estimated stochastically in O(P ) time and O(P ) memory (storing a probe vector of dimension P ) via Hessian–vector products. The theoretical analysis establishes exponential decay of resonance in vanilla networks, preservation of long-range coupling under skip connections, and rank constraints imposed by bottleneck layers. Practical applicability of the framework is confirmed on a convolutional architecture (ResNet-18, ∼11M parameters) through stochastic curvature estimation. Special case: if the architecture reduces to a single node with no parents or children, all proposed definitions naturally collapse to the standard Hessian ∇2θ L(θ) ∈ Rp×p . Scope. The formalism is designed for structural analysis, not for novel optimization algorithms. Experiments 1–5 verify theoretical predictions on fully connected MLPs and toy Attention; Exp. 6 extends the validation to a convolutional architecture (ResNet-18, ∼11M parameters). Scaling to ImageNet-scale models and Transformers is discussed in Section 6. 2

2. Related Work Second-order curvature approximations. Methods based on the GN approximation (Martens, 2010; Martens and Sutskever, 2012; Botev et al., 2017) and its Kronecker variants (K-FAC (Martens and Grosse, 2015), KFRA (George et al., 2018), Shampoo (Gupta et al., 2018)) operate on block-diagonal structures and ignore both the tensor component H T and inter-layer cross-blocks. AdaHessian (Yao et al., 2021) estimates only the Hessian diagonal; Sophia (Liu et al., 2024) uses a stochastic diagonal estimate for large language models. Pearlmutter HVP (Pearlmutter, 1994) yields exact Hessian–vector products but does not reveal block structure. Abreu et al. (2025) studied the full GN for LLMs, confirming the practical significance of the decomposition. Kunstner et al. (2019) demonstrated the Fisher/GGN/Hessian divergence; our decomposition H = H GN +H T explicitly separates these contributions and quantifies the gap via the GN-Gap. Dangel et al. (2022) exploit the low-rank structure of the GGN for curvature access; our path decomposition reveals the architectural origin of this structure. Modular curvature backpropagation. Botev et al. (2017) proposed recursive computation of GGN blocks via backpropagation (modular curvature backpropagation); Dangel et al. (2020) scaled this approach in the BackPACK framework, implementing modular computation of layer-wise GGN blocks and their Kronecker approximations for sequential architectures. Our formalism extends the modular paradigm: (i) interf layer cross-blocks Hv,w (v ̸= w), absent in layer-wise schemes; (ii) arbitrary DAGs with fan-in > 1; (iii) an explicit tensor component H T with quantitative diagnostics (GN-Gap, R, C, D). A key consequence: the structural theorems of this work—exponential resonance decay (Theorem 25), path decomposition (Theof , v ̸= w, and rem 33), bottleneck rank constraint (Proposition 31)—are formulated through cross-blocks Hv,w have no analog in block-diagonal schemes operating only on Hv,v . Spectral analysis and loss landscape. Papyan (2019), Pennington and Worah (2018), Ghorbani et al. (2019), and Sagun et al. (2017) studied the global Hessian spectrum; Tang et al. (2025) and Singh (2025) analyzed Hessian structure in CNNs and LLMs. Per-layer Hessian spectra have been linked to local curvature properties and generalization (Bolshim and Kugaevskikh, 2025; Sagun et al., 2017); the present work generalizes the per-layer approach from diagonal blocks to the full inter-layer block structure, enabling diagnostics of specific layer pairs. Schoenholz et al. (2017) and Hayou et al. (2019) established edge of chaos conditions for signals; our condition sρ < 1 (Corollary 26) generalizes this idea to curvature in DAG architectures. Cohen et al. (2021) empirically showed that GD operates at the edge of the maximum Hessian eigenvalue (edge of stability); our analysis provides a formal framework for studying this phenomenon at the level of inter-layer blocks (Theorem 25). In the NTK limit (Jacot et al., 2018) H T ≡ 0; our analysis operates at finite width where H T ̸= 0. Non-smooth activations. Clarke (1990) and Bolte and Pauwels (2021) provided the formal foundation of generalized derivatives and CSVF for neural networks. Zhang et al. (2018) proposed computing local Hessian blocks during backpropagation, but without explicit inter-layer cross-blocks for DAGs. f for arbitrary Positioning. Unlike GN/K-FAC/HVP, we construct full inter-layer cross-blocks Hv,w GN T DAGs with the decomposition H +H , path analysis, and quantitative curvature metrics; unlike NTK/ spectral analysis, we operate at finite width for specific layer pairs. 3. Methodology 3.1. Analysis framework Remark 1 (Scope and paradigm of analysis). The object of analysis is the AD-Hessian—the numerical object returned by modern AD frameworks when computing second derivatives. For smooth activations it coincides with the classical Hessian; for piecewise-linear activations (ReLU) AD assigns σ ′′ (0) := 0, annihilating the tensor component H T a.e. Correctness is justified via CSVF theory (Bolte and Pauwels, 2021). The goal of the formalism is an explicit decomposition of the Hessian into inter-layer cross-blocks, pathwise contributions, and tensor components over the DAG (curvature routing), analogous to backpropagation routing gradients. The full matrix (O(P 2 )) is intended for analysis of small architectures; for scaling, the framework collapses to an O(P ) HVP operator (Corollary 21). 3

3.2. Notation and network model Core notation is summarized in Table 1; detailed index conventions and tensor contraction rules are deferred to Appendix A.

Symbol

Table 1: Core notation. Definition

Graph and derivatives G = (V, E) DAG of the neural network Pa(v), Ch(v) Parents / children of node v fv ∈ Rdv Output of node v θv ∈ Rpv Parameters of node v δv Gradient ∇fv L Du←v Jacobian ∂fu /∂fv Tu;v , Tu;v,w Hessians / mixed second derivatives w.r.t. inputs f Hv,w Input Hessian block 2f ∂C Clarke Hessian v Decomposition and metrics (Def. 23, 24) H GN , H T GN and tensor components f ∥F R(v, w) Inter-layer resonance ∥Hv,w C(v, w) Geometric coupling (normalized R) D(v, w) Stable rank equation 7 T ∥ /(∥H GN ∥ + ϵ) equation 4 GN-Gap ∥Hv,w F v,w F

We distinguish two regularity cases (detailed definitions of function spaces are in Appendix A). Assumption 2 (Regularity of node functions). For every node v ∈ V : (A) in the smooth case gv ∈ C 2 ; (B) in the non-smooth case gv ∈ P C 2 and is locally Lipschitz. Architecture. A neural network is specified by a DAG G = (V, E). For each node v ∈ V : output fv = gv (fPa(v) , θv ) ∈ Rdv , parameters θv ∈ Rpv , loss L : Rdout → R. Computation proceeds in reverse topological order over G. Gradients δv , Jacobians Du←v , Dv , and second-derivative tensors Tu;v , Tu;v,w are defined in Table 1; element-wise definitions, contraction rules, and notational conventions are in Appendix A (notation follows Magnus and Neudecker, 2019). Remark 3 (Element-wise form of the tensor terms). For a node u with a single parent v and node function fu = σ(Wu fv + bu ), the second-derivative tensor has entries [Tu;v ]i,j,k =

∂ 2 fu,i = σ ′′ (zu,i ) Wu,ij Wu,ik , ∂fv,j ∂fv,k

(1)

where zu = Wu fv + bu . For a node u with two parents v, w ∈ Pa(u) (fan-in ≥ 2), the mixed tensor [Tu;v,w ]i,j,k = ∂ 2 fu,i /(∂fv,j ∂fw,k ) captures cross-input curvature; it vanishes for linear merge (fu = fv + fw ) ′′ and P is nonzero when the merge involves a nonlinearity with σ ̸= 0 (cf. Example 6). The contraction i [Tu;v ]i,•,• δu,i appearing in equation 2 is a dv ×dv matrix weighted by the gradient at node u. 3.3. Non-smooth case: AD-Hessian for piecewise-linear networks For piecewise-linear activations (ReLU, Leaky ReLU, max-pooling) σ ′′ (z) = δ(z); AD frameworks assign ′′ σ (0) := 0, so that the tensor component Tu;v ≡ 0 for all such nodes. By the canonical decomposition T GN (Theorem 7), Hv,w ≡ 0 a.e.: the input Hessian coincides with Hv,w ⪰ 0 in activation space; the parametric residual does not vanish in general (Remark 16). For smooth activations (GELU, Swish) H T ̸= 0, and the GN-Gap (Definition 11) quantitatively measures the curvature lost under the GN approximation. Correctness of AD computations for non-smooth networks is justified via CSVF theory (Bolte and Pauwels, 2021): the AD-Hessian coincides with the classical one a.e., is bounded at non-smooth points, and preserves convergence of stochastic estimates (formal definitions and proofs are in Appendix D). 4

3.4. Full input Hessian f Definition 4 (Input Hessian). The full input Hessian is the block matrix {Hv,w }v,w∈V , where each block f dv ×dw Hv,w ∈ R is defined recursively: f The input Hessian Hv,w is the sum of four components: f Hv,w =

Du⊤1 ←v Huf1 ,u2 Du2 ←w

u1 ∈Ch(v) u2 ∈Ch(w) du X

+ 1v̸=w

[Tu;v,w ]i,•,• δu,i

u∈Ch(v)∩Ch(w) i=1

+ 1v=w

du X X [Tu;v ]i,•,• δu,i u∈Ch(v) i=1

∂ L ∂fv ∂fw

(2)

The terms correspond to: (1) the full double sum over Ch(v) × Ch(w), accounting for all cross-blocks; (2) mixed tensor inputs (when v ̸= w); (3) pure tensor terms for a single input (when v = w); (4) direct loss dependence. Example 5 (Two-layer network). Consider v1 → v2 → out with no skip connections: Ch(v1 ) = {v2 }, Ch(v2 ) = {out}. For v = w = v1 the double sum in equation 2 reduces to a single term u1 = u2 = v2 (the unique pair in Ch(v1 ) × Ch(v1 )): Hvf1 ,v1 = Dv⊤2 ←v1 Hvf2 ,v2 Dv2 ←v1 +

dv2 X [Tv2 ;v1 ]i,•,• δv2 ,i . i=1

The first term propagates curvature from the descendant v2 ; its GN component is positive semi-definite (PSD) and is given via Dout←v1 (Theorem 7), but the full Hvf2 ,v2 also contains a tensor part. The second term is the tensor part (may be indefinite). 2 L ∂2L = Cross-block Hvf1 ,v2 : since fv1 affects L only through v2 , by the chain rule ∂fv∂ ∂f = Dv⊤2 ←v1 ∂f 2 v 1

Dv⊤2 ←v1 Hvf2 ,v2 . This is a special case of the one-directional path (Appendix H, Remark 76). With a skip connection v1 → out, term (1) of equation 2 for the pair (v1 , v1 ) would include all pairs (v2 , v2 ), (v2 , out), (out, v2 ), (out, out) from Ch(v1 ) × Ch(v1 ), accounting for the cross-block Hvf2 ,out . Example 6 (Diamond graph). Consider a DAG with branching: v1 → v2 , v1 → v3 , v2 → v4 , v3 → v4 —two branches converging at node v4 . For the cross-branch pair (v2 , v3 ): Ch(v2 ) ∩ Ch(v3 ) = {v4 } (common child with fan-in ≥ 2). Equation equation 2 gives: Hvf2 ,v3 = Dv⊤4 ←v2 Hvf4 ,v4 Dv4 ←v3 +

dv4 X [Tv4 ; v2 ,v3 ]i,•,• δv4 ,i . i=1

The first term is the GN component, propagating curvature through v4 (always present). The second is the mixed tensor term (2) of equation 2: it vanishes for linear merging (fv4 = fv2 + fv3 , since Tv4 ;v2 ,v3 ≡ 0 by linearity), but does not vanish for nonlinear merging with σ ′′ ̸= 0. Unlike Example 5, here v2 and v3 are not connected by a directed path—the cross-block Hvf2 ,v3 arises exclusively through the common descendant v4 . This is the minimal topology in which term (2) can be nonzero.

with base cases: f Hout,out

= ∇2 L(fout ), X f ⊤ = Du←v Hu,out

f Hv,out

(∀v ̸= out),

u∈Ch(v) f Hout,v

f = Hv,out

⊤

Since node out has no children (Ch(out) = ∅), term (1) of equation 2 does not apply to pairs involving w = f out. Instead, blocks Hv,out are computed by a one-sided recursion (Remark 76 in Appendix), which closes from out in reverse topological order. 3.5. Canonical decomposition of the Hessian Equation equation 2 naturally splits into two structurally distinct components, each with a clear geometric interpretation. f Theorem 7 (Canonical decomposition). The input Hessian Hv,w admits a canonical decomposition:

(3)

f GN T Hv,w = Hv,w + Hv,w ,

where the components are defined as follows: GN (self-contained recursion): 1. Gauss–Newton component Hv,w GN Hv,w =

Du⊤1 ←v HuGN Du2 ←w 1 ,u2

u1 ∈Ch(v) u2 ∈Ch(w)

∂2L . ∂fv ∂fw

GN ⊤ GN Unrolling the recursion to the output node: Hv,w = Dout←v ∇2L Dout←w . The block matrix [Hv,w ]v,w∈V is positive semi-definite (see Remark 8 below). GN f T satisfies the recursion with base case − Hv,w := Hv,w 2. Tensor component (residual curvature) Hv,w T = 0: Hout,out T Hv,w =

Du⊤1 ←v HuT1 ,u2 Du2 ←w

u1 ∈Ch(v) u2 ∈Ch(w)

+ 1v̸=w

du X

[Tu;v,w ]i,•,• δu,i

u∈Ch(v)∩Ch(w) i=1

+ 1v=w

du X X [Tu;v ]i,•,• δu,i . u∈Ch(v) i=1

This component encodes node function curvature and can be either positive or negative semi-definite. Proof sketch. The recursion for H GN is self-contained (uses HuGN , not full H f ); H T := H f − H GN inherits 1 ,u2 T GN ⊤ the recursion with base Hout,out = 0. Unrolling to the output: Hv,w = Dout←v ∇2L Dout←w ; PSD follows from Remark 8. Full proof: Appendix F. L GN Remark 8 (PSD of the block GN matrix). When Hout ⪰ 0 the block matrix HGN := [Hv,w ]v,w∈V is PSD: P ⊤ GN 2 L r H r = ∥ v Dout←v rv ∥H L ≥ 0. PSD is inherited from the output Hessian Hout ; individual blocks out

f Hu,u (which include H T ) need not be PSD.

Remark 9 (H GN and the GGN). In unrolled form: GN ⊤ Hv,w = Dout←v ∇2L Dout←w .

Let J := [Dout←v1 , . . . , Dout←vn ]⊤ be the full network Jacobian. Then HGN = J ⊤∇2L J, which is precisely the Generalized Gauss–Newton (GGN) matrix (Schraudolph, 2002; Martens and Grosse, 2015). One should distinguish GGN from the Fisher matrix: GGN is a deterministic construction requiring no probabilistic model; under log-likelihood loss GGN coincides with the Fisher, but in general they differ (Martens, 2020). The canonical decomposition represents the full Hessian as H = GGN + H T , explicitly isolating the part discarded by optimizers such as K-FAC. Proposition 10 (Input Hessian for piecewise-linear networks). For networks with piecewise-linear activaT tions (ReLU, Leaky ReLU, max-pool) the tensor component of the input Hessian Hv,w ≡ 0 a.e., and in f GN activation space Hv,w = Hv,w for all pairs (v, w). The parametric Hessian equation 5 contains a residual term determined by node function structure but independent of activation curvature σ ′′ ; in parameter space the full Hessian does not reduce to GGN in general (Remark 16). Proof sketch. For piecewise-linear activations, σ ′′ (z) = 0 a.e. By Remark 3, [Tu;v ]i,j,k = σ ′′ (zu,i ) Wu,ij Wu,ik = 0 a.e. for every activation node u; likewise Tu; v,w = 0. In the recursion for H T (Theorem 7, part 2) all source T terms (tensor summands) vanish, and the base case is Hout,out = 0. By structural induction in reverse T topological order, Hv,w = 0 for all (v, w). Full proof: Appendix D. Definition 11 (GN-Gap). The deviation of the full Hessian from the Gauss–Newton approximation is measured by the GN-Gap: T ∥Hv,w ∥F GapGN (v, w) := , (4) GN ∥Hv,w ∥F + ϵ where ϵ > 0 is a small constant for numerical stability. Corollary 12 (Interpretation of GN-Gap). GN-Gap quantifies the quality of the Gauss–Newton approximation for a given layer pair: • GapGN (v, w) ≈ 0: the GN approximation is accurate, and optimizers such as K-FAC are expected to be close to exact second order. • GapGN (v, w) ≫ 0: a substantial contribution of the tensor component H T is present; the full Hessian is required for precise optimization. Since GN-Gap uses the Frobenius norm, it measures the magnitude of the ignored tensor component, not the sign of curvature. To assess negative curvature specifically, one should use the negative-curvature mass (Appendix H). Scope of applicability. For piecewise-linear activations (ReLU) the tensor component vanishes (Proposition 10), and GapGN ≈ 0 in activation space; the non-convexity of the landscape concentrates on linearregion boundaries (as delta contributions) unobservable by AD. GN-Gap is critically informative for networks with smooth activations (GELU, Swish, Softmax in Attention), where H T ̸= 0 everywhere and the GN approximation systematically loses information about negative curvature. Remark 13 (Relation to practical optimizers). Within decomposition equation 3, approximations such as GN K-FAC (Martens and Grosse, 2015) correspond to dropping the tensor component H T . Since [Hv,w ]⪰0 (Remark 8), such approximations contain no information about negative curvature, which may limit their effectiveness near saddle points (Dauphin et al., 2014). f Proposition 14 (Structural sparsity and routing). The block Hv,w is nonzero if and only if there exists a path from v and w to a common descendant u, or a direct dependence ∂ 2 L/∂fv ∂fw ̸= 0. case PFor the special f ⊤ f of a one-directional path (v →∗ w) the computation reduces to the recursion Hv,w = u∈Ch(v) Du←v Hu,w , which excludes tensor terms.

Proof. Follows from the topological structure of the graph and the chain rule (see Appendix F). f f Symmetry property. In the smooth case (Case A) Hv,w = (Hw,v )⊤ for all v, w ∈ V , following from equality of mixed partial derivatives for twice continuously differentiable functions. In the non-smooth case (Case B), at points of non-differentiability, symmetry of the AD-Hessian may f f f + (Hw,v )⊤ ). not hold. In such cases one can symmetrize: Ĥv,w = 12 (Hv,w f f f Remark 15 (Symmetrization in the non-smooth case). Symmetrization Ĥv,w = 12 (Hv,w + (Hw,v )⊤ ) modifies the spectrum and is recommended when positive definiteness or Cholesky factorization is required; it should be avoided if asymmetry carries curvature information at non-smooth points.

3.6. Full parametric Hessian f Given the input block Hv,w , the parametric Hessian ∇2θ L = {Hθv ,θw } is obtained by the standard chain rule: X X f (i) Hθv ,θw = D⊤ δv,i Tv(i) + δu,i D⊤ (5) v Hv,w Dw + 1v=w v Tu;v Dw , i

u,i

(i) (i) where Tv = ∂ 2 fv,i /∂θv2 and Tu;v = ∂ 2 fu,i /(∂fv ∂θv ) are the pure-parameter and mixed input–parameter

second-derivative tensors, respectively. With weight sharing an additional cross-tensor term appears (full expansion: Appendix P F.1, Eq. equation Appendix F.1). Weight tying (Pascanu et al., 2013) generalizes trivially: Hθ,θ = a,b Hθva ,θvb (Appendix A). Remark 16 (Activation vs. parameter space). The canonical decomposition H f = H GN +H T (Theorem 7) is T ≡ 0 (piecewise-linear case) it does not follow that Hθv ,θw = HθGN : defined for the input Hessian. From Hv,w v ,θw (i)

(i)

the second and third terms in equation 5 involve Tv and Tu;v , which depend on the node function structure and are nonzero for linear layers (∂ 2 (W fv )/(∂fv ∂θv ) ̸= 0). Contribution C3 refers to the decomposition of the input Hessian; full expansion of the parametric terms is in Appendix F.1. Corollary 17 (Input-to-parameter Hessian bridge). From equation 5 one obtains a norm bound for the parametric cross-block via the input one: f ∥Hθv ,θw ∥F ≤ ∥Dv ∥2 ∥Dw ∥2 ∥Hv,w ∥F + Rv,w ,

where Rv,w is the contribution of the pure-parameter and mixed tensors (second and third terms) in equation 5. Consequently, the metrics R(v, w) and C(v, w), defined for activations, upper-bound the corresponding parametric cross-blocks up to the Jacobian scale factors ∥Dv ∥2 , ∥Dw ∥2 . High inter-layer resonance R(v, w) indicates potentially strong coordination of weight updates for θv and θw . 3.7. Curvature routing and HVP operator Formulas equation 2 and equation 5 define the structure of the inter-layer Hessian. Just as classical backpropagation routes gradient signals δv = ∇fv L over the DAG, the recursion equation 2 extends this f paradigm to matrices of curvature signals—blocks of the inter-layer Hessian Hv,w . The Pearlmutter trick (HVP) (Pearlmutter, 1994) treats the Hessian as a monolithic operator; the graph recursion decomposes it into cross-blocks, path-wise contributions, and tensor components. When collapsed over all nodes (Corollary 21), it reduces to an O(P ) HVP. Definition 18 (Curvature routing recursion). Let G = (V, E) be the DAG of a neural network with loss f L. Curvature routing is defined as the scheme of propagating curvature matrices {Hv,w }v,w∈V over G in f 2 reverse topological order, starting from the base condition Hout,out = ∇ L(fout ), computing boundary blocks f Hv,out by one-sided recursion (Remark 76), and applying the recurrence equation 2 for pairs of internal nodes (Figure 1). Theorem 19 (Recursion completeness). For a neural network with DAG G = (V, E) and L ∈ C 2 , the recursion equation 2, equation 5 exactly recovers ∇2θ L = [Hθvi ,θvj ]ni,j=1 . Proof. Correctness is guaranteed by induction over the reverse topological order of graph nodes using the generalized chain rule. Full formal proof: Appendix I. 8

(a) DAG & forward paths

(b) Backward pass: curvature routing Hvf ,out = D1⊤ Hvf ,out + I ⊤ Hvf ,out

I (skip)

out

D1⊤

v1 I

D2⊤

D3⊤

out

⊤ f Hvf ,out = D3⊤ ∇2L Hout,out = ∇2L 2

Dout←v0 = (I + D2 D1 ) D3

Hvf ,out = D2⊤ Hvf ,out 2 1

Figure 1: (a) DAG of a neural network with skip connection (v0 → v2 ): forward edges (solid) and skip edge (dashed). The full f are propagated Jacobian is a sum of path-wise contributions. (b) Backward pass of curvature routing: curvature matrices Hv,w f from out over the graph via transposed Jacobians; with a skip connection, the block Hv0 ,out receives contributions from both paths.

Operator formulation: Hessian–vector products The recursion equation 2 admits an operator formulation that avoids storing the full P ×P matrix. Definition 20 (HVP operator). A tangent vector rv is propagated by a forward pass: X rv = Dv←w rw , v ∈ / Vin ,

(6)

w∈Pa(v)

and the HVP operator action is defined as HVPv (r) :=

f w∈V Hv,w rw .

When collapsed over nodes, the recursion reduces to an O(P ) HVP (Pearlmutter, 1994) with explicit decomposition H GN +H T (algorithm and proof: Appendix I; AD implementation: Appendix E). Corollary 21 (HVP complexity). The computational cost of HVP(r) for a tangent vector r obtained by the forward pass equation 6 coincides with the cost of one backward pass (up to a multiplicative constant). For ⊤ r costs O(du dv ), the total complexity is O(P ), where P is fully connected layers, where multiplication Du←v the total number of parameters. ⊤ by a vector of dimension du , costing Proof. At each edge (v, u) the mainP operation is multiplication of Du←v O(du dv ). Summing over all edges: (v,u)∈E du dv , which coincides with the analogous sum for the backward pass.

Remark 22 (Relation to AD and practical significance). The recursive HVP formula is implemented via double application of AD (Pearlmutter, 1994): a forward pass, backward for δv , then backward-over-forward with JVP in direction r. The graph recursion provides a micro-decomposition of each HVP summand via equation 3, enabling isolation of H GN and H T contributions and underpinning second-order optimization algorithms (Newton-CG, L-BFGS, Trust Region) for arbitrary DAG architectures. 4. Diagnostics and structural analysis f The block structure {Hv,w } and the canonical decomposition induce quantitative metrics and structural theorems that connect the geometry of L with the network architecture.

Definition 23 (Inter-layer resonance and geometric coupling). For a pair of nodes v, w ∈ V we define f the inter-layer resonance R(v, w) := ∥Hv,w ∥F —a measure of joint curvature influence, and the geometric p coupling C(v, w) := R(v, w) R(v, v) · R(w, w)—a normalized connectivity measure. In the PSD regime (H T = 0) the metric is strictly bounded: C ≤ 1. A value C > 1 is mathematically possible only due to the tensor component H T and serves as an indicator of an indefinite contribution of H T to the curvature 9

between layers v and w. The metric is invariant to rescaling of adjacent weights Wv 7→ αWv , Ww 7→ α−1 Ww for v ∈ Pa(w) (Theorem 44 in Appendix). Since the absolute value R(v, w) naturally depends on layer dimensions dv × dw , the normalized metric C(v, w) is used for correct comparison across layers of different width. Definition 24 (Stable rank of inter-layer block). For a pair of nodes v, w ∈ V we define the stable rank of the inter-layer block: f ∥Hv,w ∥2F . (7) D(v, w) := f ∥Hv,w ∥22 The metric measures the effective dimensionality of curvature interaction between layers: f 1 ≤ D ≤ rank(Hv,w ), with D = 1 only if the block has rank 1. Unlike deff := ∥H∥∗ /∥H∥2 (Definition 57), the stable rank admits stochastic estimation via HVP: ∥H∥2F by the Hutchinson estimator (Avron and Toledo, 2011), ∥H∥22 by power iteration (Nocedal and Wright, 2006); total cost is O (m + 2T ) · Backprop , memory √ O(P ). The two-sided bound D ≤ deff ≤ r D makes D a practical replacement for deff in large-scale analysis (Proposition 62, Algorithm 4 in Appendix G.2.2). Theorem 25 (Lyapunov decay of resonance). Let v0 → v1 → · · · → vL be a sequential chain in the DAG with |Ch(vi )| = 1 for 0 ≤ i < L, Πj←i = Dvj ←vj−1 · · · Dvi+1 ←vi the product of inter-layer Jacobians, and (m)

λ1

1 log Πi+m←i 2 m

(8)

the finite-depth Lyapunov exponent. Then for 0 ≤ i < j ≤ L: (j−i)

R(vi , vj ) ≤ Cj Πj←i 2 = Cj e (j−i) λ1

(9)

f ∥F ; for piecewise-linear activations this coincides where Cj = ∥Hvfj ,vj ∥F . At j = L we have CL = ∥Hout,out GN 2 with C = ∥∇ L∥F .

Proof. Remark 76 with |Ch(vk )| = 1 gives Hvfk ,vj = Dv⊤k+1 ←vk Hvfk+1 ,vj . Applying (j−i) times from k = i to k = j−1: f Hvfi ,vj = Dv⊤i+1 ←viDv⊤i+2 ←vi+1 · · · Dv⊤j ←vj−1 Hvfj ,vj = Π⊤ j←i Hvj ,vj , whence R(vi , vj ) ≤ ∥Πj←i ∥2 Cj . Corollary 26 (Exponential curvature sensitivity for arbitrary DAGs). Let s = maxv |Ch(v)| be the maximum out-degree and ∥Du←v ∥2 ≤ ρ for all edges (v, u) ∈ E, with sρ < 1. Then for GN-resonance GN RGN (v, w) := ∥Hv,w ∥F : RGN (v, w) ≤ C GN (sρ)dist(v,w) , (10) GN where C GN = ∥Hout,out ∥F = ∥∇2L∥F . For s = 1 (sequential chains) the bound follows from Theorem 25 by submultiplicativity ∥Πj←i ∥2 ≤ ρj−i ; for arbitrary DAGs (s > 1) the proof is in Appendix H.

Remark 27 (Multiplicative ergodic theorem). Bound equation 9 is strictly no weaker than equation 10 for s = 1; the gap can be exponential (at ρ = 2.6, L = 8: ρ8 ≈ 780 vs. ∥Π∥2 ≈ 5). If the Jacobians {Di } form a (m) stationary ergodic sequence, then by the Oseledets theorem (Oseledets, 1968) λ1 = limm→∞ λ1 exists a.s. (λ1 +o(1))L and R ≤ C e ; the condition λ1 < 0 is necessary and sufficient for decay. Corollary 28 (Block-banded structure). When sρ < 1 the GN-Hessian is approximately block-banded: GN blocks Hv,w with dist(v, w) > ⌈log ε/ log(sρ)⌉ are negligibly small, justifying block-diagonal approximations (Martens and Grosse, 2015). For piecewise-linear networks (Proposition 10) this extends to the full f input Hessian [Hv,w ].

Corollary 29 (K-FAC approximation error). K-FAC (Martens and Grosse, 2015) approximates Hθ ≈ blkdiag(Hθv ,θv ), discarding all cross-blocks Hθv ,θw , v ̸= w. By Corollary 28, when sρ < 1 the approximation error is exponentially small for distant pairs and concentrates in blocks with dist(v, w) ≤ ⌈log ε/ log(sρ)⌉. Thus, the framework quantitatively characterizes the accuracy region of block-diagonal optimizers. Remark 30 (Tightness of the sρ < 1 condition). The condition sρ < 1 is sufficient but not necessary. When sρ ≥ 1, the finer bound of Theorem 25 via the Lyapunov exponent λ1 remains applicable: λ1 < 0 is necessary and sufficient for exponential decay (Remark 27). In practice (Table 2), ρmax ≫ 1 at He initialization, but λ1 < 0 due to decorrelation of singular bases. Beyond norm decay, narrow layers constrain the dimensionality of coordination between distant blocks. If all paths between v and w pass through a bottleneck of width du , then the Jacobian factors through Rdu , and the rank of the GN component inherits this compression. Proposition 31 (Rank bottleneck). If all paths from v and w to out pass through a node u with du ≪ GN min(dv , dw ), then the Gauss–Newton component is strictly bounded: rank(Hv,w ) ≤ du . Consequently, the f T full Hessian Hv,w is low-rank up to the tensor perturbation Hv,w —training coordination is limited to du directions in the Gauss–Newton regime. (Follows directly from submultiplicativity of rank under matrix multiplication (Horn and Johnson, 2012).) Proof sketch. Since every directed path from v (resp. w) to out traverses u, the end-to-end Jacobian factors as Dout←v = Dout←u Du←v , where Du←v ∈ Rdu ×dv . Substituting into the Gauss–Newton block gives GN ⊤ 2 Hv,w = D⊤ u←v Dout←u ∇ L Dout←u Du←w . {z } | ∈ Rdu ×du

By submultiplicativity of rank under matrix products (Horn and Johnson, 2012), the outer factors cannot GN ) ≤ du . Full details are in Appendix H. increase rank beyond the inner du ×du core, hence rank(Hv,w Remark 32 (Skip connections and the rank constraint). Proposition 31 applies when all paths pass through node u. With P skip connections bypassing u, additional terms appear in the path decomposition equation 11, and rank( p Dp⊤ (· · · )Dp ) may exceed du . Thus, skip connections not only restore decaying resonance (Theorem 34) but also break “information bottlenecks,” preserving full-rank inter-layer Hessians. Corollary 26 and Proposition 31 show that resonance decays exponentially and narrow layers constrain coordination. The stable rank D(v, w) (Definition 24) quantitatively measures the effective number of coordination directions. The decay mechanism is revealed by the path decomposition. Theorem 33 (Path decomposition Q of the Hessian). Let P(v → c) be the set of directed paths from v to a common descendant c, and Dp = (ui ,ui+1 )∈p Dui+1 ←ui the path-wise Jacobian. Then f Hv,w =

T Dp⊤v HfLc Dpw + Hv,w ,

(11)

c∈Desc(v)∩Desc(w) pv ∈P(v→c) pw ∈P(w→c) T where Hv,w is the tensor component of Theorem 7; its path-wise expansion is given in Appendix G. Each pair of paths (pv , pw ) to a common descendant c contributes additively to resonance R(v, w).

Proof: Appendix G. Equation equation 11 explains why skip connections restore decaying resonance: adding an edge creates new paths and thereby additional summands in R(v, w). Theorem 34 (Resonance decay: vanilla network vs. Pre-Activation ResNet). Consider an L-layer network with ∥Dvi+1 ←vi ∥2 ≤ ρ < 1. Parts (b), (b′ ) apply to the Pre-Activation ResNet architecture (He et al., 2016: xl+1 = xl + F (ReLU(xl ))) with din = dout for each residual block, where the skip edge has Jacobian I: L→∞ (a) Vanilla network (upper bound): R(v0 , vL ) ≤ CρL −−−−→ 0. 11

(b) ResNet with identity skips every k layers (upper bound): R(v0 , vL ) ≤ Cρk (independent of L). (b′ ) ResNet (lower bound, PSD regime): If H T = 0 and ∇2L ⪰ 0, then R(v0 , vL ) ≥ (1−ρk )L/k ∥∇2L∥F > 0. k

The decay rate ρk is exponentially slower than the vanilla rate ln ρ1 when ρk ≪ 1. With residual branch scaling α = c/L: (1−(αρ)k )L/k → 1, i.e. the lower bound is independent of L. Proof sketch. (a) Submultiplicativity gives ∥ΠL←0 ∥2 ≤ ρL ; applying Theorem 25 yields R ≤ CρL → 0. (b) With identity skips every k layers, the path decomposition equation 11 includes a direct path whose QL/k Jacobian contribution is the product j=1 (I + αk Mj ) where ∥Mj ∥2 ≤ ρk . The upper bound becomes R ≤ C(1 + ρk )L/k ; when α = O(1/L), this converges to (1 + (αρ)k )L/k → e0 = 1. (b′ ) In the PSD regime (H T = 0, ∇2L ⪰ 0), the skip-path summands are individually PSD; by Weyl’s inequality each contributes a positive-semidefinite lower bound via σmin , yielding R ≥ (1 − ρk )L/k ∥∇2L∥F . Full proof: Appendix G. Remark 35 (Edge of chaos for curvature). Corollary 26 shows that vanilla networks require strict balancing ρ = 1 (He/Xavier init, BatchNorm) to prevent exponential decay/explosion of the Hessian—an analog of “edge of chaos” (Balduzzi et al., 2017). Theorem 34 shows that skip connections remove this requirement: the lower bound R ≥ (1−ρk )L/k ∥∇2L∥F at α = O(1/L) is independent of L. Remark 36 (Practical implications). For architectures with dense resonance (ResNet, DenseNet), crossblocks Hθv ,θw with dist(v, w) ≤ k are not negligibly small, so block-diagonal optimizers (K-FAC) lose coupling information; banded preconditioners are optimal (George et al., 2018). The path decomposition predicts resonance growth in DenseNet and high C in U-Net (skip encoder→decoder); applicability to Attention (Example 37), GN insensitivity, BatchNorm—Appendix M, Appendix H. Example 37 (Attention block as DAG: density of cross-blocks). A√single-head Attention defines a DAG with nodes {vin , vQ , vK , vV , vout }, where vout computes O = σ(QK ⊤/ dk ) V (Softmax×Value) with Pa(vout ) = {vQ , vK , vV } (fan-in = 3). Softmax is a smooth nonlinearity with a dense Jacobian Sij = σi (δij − σj ) ̸= 0 f for all i, j (since σi > 0) and nonzero Hessian Tσ; z . By equation 2, the cross-block HQ,K contains the mixed T ̸= 0. tensor term (2) through common descendant vout : Tout; Q,K ̸= 0, hence HQ,K Structural implications: (i) Softmax creates dense cross-blocks across all sequence positions—the blockbanded approximation (Corollary 28) is inapplicable inside an Attention block; (ii) GN-Gap for the pair (Q, K) is fundamentally nonzero, qualitatively distinguishing Attention from piecewise-linear layers and rendering the Gauss–Newton approximation inexact; (iii) multi-head analysis reduces to block diagonalization over heads and requires no new formalism. Full Jacobian/tensor expansions and proofs: Appendix M. 5. Experimental validation of diagnostic metrics The purpose of the experiments is falsification-based verification: we test whether empirical data can refute the theoretical predictions, rather than propose a new algorithm. Six experiments are conducted. Experiments 1–4 use CIFAR-10 (CIFAR-100 for Exp. 2) with fully connected MLPs on the flattened input (din = 3072; channel normalization to [−1, 1], no augmentation), enabling exact Hessian block computation. Experiment 5 uses synthetic data and a minimal Attention model (§ 5.5). Experiment 6 extends the validation to a convolutional architecture (ResNet-18, ∼11M parameters) on CIFAR-10 via stochastic curvature estimation (§ 5.6). Each configuration is trained with 5 random seeds {42, . . . , 46}; metrics are computed at checkpoints: init (before training), mid (half epochs), final. Common protocol. Training: SGD (η = 0.01, momentum 0.9, weight decay 10−4 ), cosine schedule, 50 epochs, batch size 128, gradient clipping (∥∇∥2 ≤ 1). Plain MLP blocks: Linear–LayerNorm–σ (LayerNorm disabled in Exp. 3, where piecewise linearity is required—§ 5.3; and in Exp. 1b with spectral normalization—Appendix N); Residual MLP blocks: Linear–σ + identity skip (no LayerNorm, α = L−1/2 ).1 1 LayerNorm is a smooth nonlinearity (∇2 ̸= 0), contributing to H T even for ReLU. In Exp. 1/2/4 we measure R(d), whose decay rate is determined by the spectral radius of the Jacobian ρ and is independent of the magnitude of H T . Validation of Proposition 10 for piecewise-linear activations: Exp. 3, where LN is disabled.

Stochastic estimates: Hutchinson with 30 Rademacher probes and a common probe vector for correct estimation of ∥Havg ∥F ; subsample of 32 examples; power iteration: 20 steps for ∥H∥2 . 5.1. Exp. 1: Resonance and coupling decay Protocol. For each depth L ∈ {8, 10, 12}, Plain MLP and Residual MLP (skip at every block, α = L−1/2 ) of width 64 are trained with SGD and cosine schedule, 50 epochs (additional depths L ∈ {16, 32}: Appendix). ¯ At each checkpoint we compute R̄(d) and C(d)—averages over pairs at distance d. Results (Figure 3). (a) In Plain MLP, R̄(d) decays exponentially with d, confirming Corollary 26. Loglinear regression log R̄ = a − bd yields R2 > 0.91 for all depths L ∈ {8, 10, 12} at init and final (Table 2); for L = 12 norms decrease by 3–4 orders of magnitude (to 10−14 at L = 32, see Appendix). (b) In Residual MLP, ¯ R̄(d) stabilizes (b ≈ 0.02, two orders slower than Plain), consistent with Theorem 34.2 (c) Coupling C(d) (Figure 2) qualitatively reproduces the same trends: decay in Plain, preservation in ResNet. Condition sρ < 1 and the Lyapunov exponent. For He initialization ρmax ≈ 2.6 (Table 2), so the sufficient condition sρ < 1 of Corollary 26 is formally not satisfied, but exponential decay of R̄(d) is robustly observed. The reason is decorrelation of singular bases: at L = 8, ∥J8 · · · J1 ∥2 ≈ 5 (vs. ρ8 ≈ 780), λ1 ≈ −0.02 (Theorem 25). Exp. 1b (Appendix N) verifies the theorem in strict mode: spectral normalization ensures ρ ≤ 1, s = 1 ⇒ sρ ≤ 1. For Residual MLP, ρmax > 1 by construction (J = αDW +I); stability of R̄ is provided by skip branches (Theorem 34). Remark 38 (Batch-averaged decay beyond the Lyapunov bound). Theorem 25 bounds the per-sample resλ1 L onance: when grows with depth. The P λ1 > 0 the worst-case R(v, w; xb ) ≤ C e Q batch-averaged quantity 1 f (xb ) R̄(d) = B b R(v, w; xb ) can nevertheless decay, because the Jacobian products i Di (xb ) entering Hv,w have sample-dependent singular subspaces; the off-diagonal contributions partially cancel upon averaging, √ analogous to the O(1/ B) variance reduction that arises from summing random-sign terms. A tight bound on R̄ would require distributional assumptions on the singular subspaces of {Di (xb )}B b=1 ; characterizing the precise rate of cancellation remains an open question. U-shaped artifact. At L ∈ {16, 32} the unnormalized R̄(d) exhibits a rise for d > L/2 due to non¯ uniformity of R(v, v) near the output; the normalized C(d) decays monotonically (R2 ≥ 0.90, Table Appendix N.3 in Appendix). Table 2: Spectral norms of Jacobians ρmax , Lyapunov exponent λ1 , and test accuracy for Exp. 1 configurations (CIFAR-10, mean over 5 seeds). L

Arch

ρinit max

ρmid max

ρfinal max

λinit 1

λmid 1

λfinal 1

Acc (%)

8 8 10 10 12 12

Plain Res. Plain Res. Plain Res.

2.59 1.25 2.69 1.22 2.69 1.20

2.50 1.31 2.47 1.27 2.48 1.23

2.78 1.33 2.59 1.29 2.55 1.25

−0.02 0.02 −0.03 0.02 −0.03 0.01

0.05 0.04 0.03 0.03 0.03 0.03

0.08 0.04 0.04 0.04 0.04 0.03

49.8 51.3 49.7 51.3 49.7 51.3

5.2. Exp. 2: Bottleneck ablation Protocol. BottleneckMLP (depths L ∈ {6, 8}, base width dbase = 256, CIFAR-100, K = 100) with a single narrow layer of width du ∈ {4, 8, 16, 32, 64, 128, 256} at position L/2. Control run († ): dbase = du = 512 (uniform network without bottleneck; stochastic estimate, 100 probes, subsample 64). Metrics Dfar and Cfar : averages over cross-bottleneck pairs (i<L/2<j). 2 The experimental scaling α = L−1/2 directly verifies the upper bound (b), which is independent of both L and α. The lower bound (b′ ) guarantees depth-independent curvature preservation under the stronger condition α = O(1/L); with α = L−1/2 the √ lower bound decays as e−ρ L —still exponentially slower than the vanilla rate ρL —but is not formally L-independent. Note that Theorem 34 assumes contracting activations (ρ < 1); the experimental regime (ρmax > 1) falls outside this formal scope, so the observed stabilization is a qualitative, not quantitative, confirmation of the theorem.

L=8

¯ C(d)

0.5 Plain init Res. init

Plain final Res. final

3 4 Distance d L = 12

¯ C(d)

0.5

−1

5 6 Distance d

¯ Figure 2: Exp. 1: Geometric coupling C(d) vs. distance d (error bars: ±1σ over 5 seeds). (a) L = 8, (b) L = 12. Plain MLP: C decays monotonically from 1 to 0.24 (L = 12, init), reflecting loss of geometric coherence between distant layers. Residual MLP: C > 0.93 at all distances; skip connections preserve coupling. After training both architectures shift upward (curvature becomes more uniform).

Plain init

Plain final

Res. init

Res. final

L=8

L = 12

100

R̄(d)

100

10−1 10−1

Distance d

Figure 3: Exp. 1: Mean resonance R̄(d) vs. distance d between layers (log scale on y; error bars ±1σ over 5 seeds). (a) L = 8, (b) L = 12. Plain MLP exhibits exponential decay (straight lines on log scale, R2 > 0.91); Residual MLP shows stabilization (b ≈ 0.02). Scale differs: init ∼ 10−1 , final ∼ 100 (reflects overall curvature growth during training).

Table 3: Exp. 2: bottleneck ablation (L = 6, CIFAR-100). Stable rank Dfar and coupling Cfar vs. narrow layer width du (mean ±1σ over 5 seeds). † Control: dbase = 512 (stochastic estimate, 100 probes, subsample 64). ‡ Acc at uninformative model level (du ≪ K); metrics for du ≥ 8 are preferred for conclusions. du

init Dfar

final Dfar

init Cfar

final Cfar

Acc (%)

4 8 16 32 64 128 256

1.65 ± 0.23 2.58 ± 0.11 4.04 ± 0.65 5.67 ± 0.99 6.95 ± 0.52 9.16 ± 0.91 10.73 ± 0.84

2.03 ± 0.20 3.39 ± 0.36 4.90 ± 0.24 5.01 ± 0.41 5.15 ± 0.34 5.45 ± 0.42 5.18 ± 0.49

0.21 ± 0.06 0.35 ± 0.05 0.43 ± 0.03 0.52 ± 0.03 0.61 ± 0.03 0.66 ± 0.04 0.73 ± 0.02

0.37 ± 0.02 0.58 ± 0.03 0.68 ± 0.01 0.73 ± 0.01 0.77 ± 0.02 0.84 ± 0.01 0.86 ± 0.02

8.9‡ 14.2 17.8 19.8 21.1 22.0 22.7

512†

11.17 ± 0.77

5.18 ± 0.36

0.77 ± 0.01

0.88 ± 0.01

23.3

Results. Stable rank Dfar increases monotonically with du at initialization (from 1.7 at du = 4 to ≈ 10.7 at du = 256; suppression ≈ 6.5×) and satisfies Dfar ≤ min(du , K−1) for all configurations, verifying the rank constraint of Proposition 31. At du = 4 the model operates in severe compression (Acc ≈ 9 % at K = 100); main conclusions are robust for du ≥ 8. The control run († , Table 3) with dbase = du = 512 (uniform init = 11.2: without a narrow layer, stable rank is determined by network network, no bottleneck) yields Dfar width and class count (see § 6). Coupling Cfar is less sensitive (≈ 2×): the bottleneck constrains the rank GN of both the cross-block Hv,w and the diagonals Hv,v , Hw,w , and normalization in C partially compensates. Unnormalized resonance Rfar decreases with growing du after training, whereas at initialization Rfar ≈ const for all du ; normalized metrics C and D isolate the structural bottleneck effect from training-induced scale changes. Summary for L = 6: Table 3 and Figure 4; L = 8 data are qualitatively analogous (Dfar lower for du ≥ 64, up to 20 % at init; see Appendix N).

Dfar

100

L = 6, init L = 6, final L = 8, init L = 8, final min(du , K−1)

0 4

128

256

512

Bottleneck width du Figure 4: Exp. 2: Stable rank Dfar vs. bottleneck width du (CIFAR-100). Solid/dashed: init/final; bright: L = 6, faded: L = 8. init is lower for d ≥ 64 (∼8–20 %), consistent with additional rank compression at greater distance. Dotted: theoretical L = 8: Dfar u bound min(du , K−1). du = 512 († , control: dbase = 512, no bottleneck): Table 3.

5.3. Exp. 3: GN-Gap by activation type Motivation. The canonical decomposition (Theorem 7) predicts that the GN-Gap is governed by the activation curvature σ ′′ : GapGN ≈ 0 for piecewise-linear activations (σ ′′ = 0 a.e.) and GapGN > 0 for smooth ones (σ ′′ ̸= 0). We test this prediction across five activations under two complementary protocols. Protocol. Plain MLP (depth 6, width 64) with five activations: ReLU, LeakyReLU (piecewise-linear), Softplus, SiLU, GELU (smooth); training for 20 epochs; subsample for the Hessian: 16 examples (exact computation at width = 64, memory-bound). For each pair (v, w) the exact decomposition H = H GN +H T and GN-Gap equation 4 are computed; additionally, E[σ ′′ (z)2 ] is measured over pre-activation values. Two configurations are evaluated: (a) standard conditions with LayerNorm, where all activations reach Acc > 15

50 %; (b) isolation protocol without LayerNorm, removing the smooth LayerNorm contribution3 to isolate the pure effect of σ ′′ . Gap stabilizes between mid and final checkpoints (column ∆ in Tables 4, 5). Standard training conditions (Exp. 3b, with LayerNorm).. Under standard training all activations converge (Acc > 50 %, Table 4). For ReLU + LN, GapGN > 0 is expected: LayerNorm is a smooth nonlinearity contributing σ ′′ ̸= 0. Results (Table 4) confirm: (i) the baseline GapGN ≈ 0.2 for ReLU/LeakyReLU + LN is entirely due to the smooth LayerNorm component; (ii) smooth activations additionally increase Gap, preserving the monotonic ranking by E[σ ′′ 2 ] (ρs = 1.0). Table 4: Exp. 3b: GN-Gap with LayerNorm (L = 6, width = 64, CIFAR-10, mean ± 1σ over 5 seeds). ∆: relative change init → final. Activation Gapinit Gapfinal ∆ (%) E[σ ′′ 2 ] Acc (%) ReLU LeakyReLU Softplus SiLU GELU

0.618 ± 0.160 0.626 ± 0.161 0.459 ± 0.098 1.967 ± 0.317 2.296 ± 0.415

0.199 ± 0.030 0.212 ± 0.013 0.272 ± 0.024 0.359 ± 0.051 0.516 ± 0.117

−67.8 −66.1 −40.9 −81.7 −77.5

0 0 0.045 0.142 0.271

51.3 ± 0.3 51.5 ± 0.3 52.7 ± 0.2 52.6 ± 0.3 52.0 ± 0.4

Isolation of activation curvature (without LayerNorm).. Disabling LayerNorm yields a strict test: GapGN ≈ 0 for piecewise-linear activations is now an exact null hypothesis (σ ′′ ≡ 0 a.e. throughout the network). Without LayerNorm, smooth activations do not reach meaningful accuracy (Acc ≲ 12 %); however, the GNGap is a structural property of the computation graph (Corollary 12): it tests whether σ ′′ = 0, independently of convergence quality. (a) For ReLU and LeakyReLU, GapGN ≈ 0 at all checkpoints and pairs, including maximum distance, strictly confirming H T ≡ 0 (Corollary 12). For ReLU/LeakyReLU the tensor component is zero at all recursion levels (σ ′′ ≡ 0 a.e.), verified as ∥H T ∥F < 10−7 . (b) For smooth activations, the correlation GapGN vs. E[σ ′′ (z)2 ] manifests already at initialization (R2 > 0.9, ρs = 0.97, p < 0.01; Figure 5) and persists at all checkpoints, consistent with the scaling H T ∼ σ ′′ 2 . The sample size n = 5 is inherently limited by the number of qualitatively distinct activation classes (two piecewise-linear plus three smooth); the near-perfect rank ordering (ρs = 0.97, p < 0.01) provides strong evidence for monotonicity despite the small n. Decomposition of ∥H T ∥F and ∥H GN ∥F by distance d (Table Appendix N.5 in Appendix) shows that both norms decay T exponentially with d, and at d = L−1 the tensor component vanishes (Hd=L−1 ≡ 0), confirming the structural routing prediction. (c) Stochastic estimates (30 probes) agree with exact values within 5–10 %. Summary: Table 5. Table 5: Exp. 3: GN-Gap, isolation protocol (no LayerNorm, L = 6, width = 64, CIFAR-10, mean ± 1σ over 5 seeds). Activation

Gapinit

Gapfinal

∆ (%)

E[σ ′′ 2 ]

Acc (%)

ReLU LeakyReLU Softplus SiLU GELU

< 10−7 ≲ 10−7 0.122 ± 0.020 0.225 ± 0.005 0.358 ± 0.008

< 10−7 ≲ 10−7 0.111 ± 0.003 0.217 ± 0.008 0.345 ± 0.012

— — −8.9 −3.5 −3.4

0 0 0.058 0.246 0.618

45.1 ± 0.9 45.0 ± 0.8 10.0 ± 0.0 11.7 ± 2.1 11.9 ± 2.4

5.4. Exp. 4: Diamond MLP — activation of tensor term (2) Motivation. Experiments 1–3 use sequential architectures in which Ch(v) ∩ Ch(w) = ∅ for v ̸= w at different depths, so that term (2) of formula equation 2 (the mixed tensor Tu; v,w for fan-in ≥ 2) is identically 3 LayerNorm is itself a smooth nonlinearity with σ ′′ ̸= 0, contributing to H T even for piecewise-linear activations. The isolation protocol disables it so that GapGN ≈ 0 for ReLU/LeakyReLU serves as a strict null hypothesis.

GapGN (init)

0.4

GELU

R2 = 0.908, ρs = 0.97 SiLU

0.2 Softplus

ReLU

0 −5 · 10−2 0

5 · 10−2 0.1

0.15

0.2

0.25

0.3

0.35 ′′

0.4

0.45

0.5

0.55

0.6

0.65

0.7

E[σ (z) ] Figure 5: Exp. 3: GN-Gap (init) vs. E[σ ′′ (z)2 ] for 5 activations (isolation protocol, LayerNorm off). Linear fit R2 = 0.908; Spearman rank correlation ρs = 0.97 (p < 0.01). With n = 5, R2 has limited power, but monotonicity of ranking (ρs ≈ 1) robustly confirms the scaling H T ∼ σ ′′ 2 .

zero. To demonstrate the non-triviality of formula equation 2 on DAG topologies, we introduce the minimal Diamond MLP construction (Example 6). Protocol. Diamond MLP: stem → two parallel branches (depth k ∈ {1, 2, 3}, width 32) → merge node → head; CIFAR-10. Sweep: merge type (sum vs. cat + σ + Linear) × activation (ReLU vs. SiLU). Metrics: R̄(d) intra-branch, R̄(dgraph ) cross-branch, ρmax = maxi ∥Ji ∥2 (spectral norm of Jacobians), GNGap equation 4 for the pair (Ak−1 , Bk−1 ) through the merge node. Hypotheses. (H4.1) Intra-branch resonance R̄AA (d) decays monotonically with distance. (H4.2a) At initialization, cross-branch R̄AB (dgraph ) decays exponentially. (H4.2b) After training the decay is violated: R grows with distance due to increasing ρ and weakened sign cancellation (Remark 81). (H4.3) Linear merge: Tmerge; A,B = 0 ⇒ GapGN ≈ 0 regardless of σ. Nonlinear merge + SiLU: Tmerge; A,B ̸= 0 ⇒ GapGN > 0. Results. (H4.1) Confirmed: of 48 (config×distance) pairs, 45 exhibit monotone decay of R̄AA (d); 3 marginal deviations are within noise (< 5 %). (H4.2a) Confirmed: at initialization, cross-branch R̄AB (dgraph ) decays exponentially; R(d+2)/R(d) ≈ 0.28–0.39 across configs. Here ρmax > 1 (LayerNorm increases the spectral norm of the Jacobian), and the elementwise upper bound of Corollary 26 grows; the observed batchaveraged decay is explained by mutual sign cancellation when averaging over x. (H4.2b) Confirmed: after training, ρmax increases and cross-branch R grows with dgraph for cat_relu, sum_relu at k ≥ 2, consistent with weakened sign cancellation as ρ grows. (H4.3) Strongly confirmed: for sum-merge, ∥Tmerge ∥F ≈ 0 and Gap ≈ 0 (merge is linear ⇒ T = 0); for cat + ReLU, ∥T ∥ > 0 but σ ′′ = 0 a.e., so Gap ≈ 0 (∼ 10−7 , machine precision); for cat + SiLU, σ ′′ ̸= 0 yields Gap ∼7 orders of magnitude above baseline (separation ≈ 2·107 ). Model accuracy ≈50 % (width = 32, CIFAR-10), sufficient for curvature diagnostics. Quantitative Gap values for k = 2: Table 6 and Figure 6; varying branch depth (k = 1, 2, 3) does not change qualitative conclusions (Table Appendix N.4 in the Appendix); cross-branch R̄AB dynamics are shown in Figure Appendix N.2, Appendix N. Table 6: Exp. 4: GN-Gap at the merge node of Diamond MLP (k = 2, width = 32, CIFAR-10, Acc ≈ 50 %; mean ±1σ over 5 seeds). ∥Tmerge ∥F : Frobenius norm of the tensor component. Merge sum sum cat cat

Gapinit

Gapfinal

∥T ∥init F

∥T ∥final F

ReLU SiLU ReLU SiLU

6.7·10−8

9.3·10−8

4.1·10−9

7.0·10−8 8.4·10−8 1.33 ± 0.17

9.6·10−8 1.2·10−7 0.082 ± 0.007

4.1·10−9 6.6·10−10 6.8·10−3 ± 5.0·10−4

2.0·10−8 2.0·10−8 3.3·10−8 2.5·10−2 ± 1.1·10−3

T 5.5. Exp. 5: Toy-Attention — verification of HQ,K ̸= 0 T Motivation. Experiments 1–4 are restricted to MLP topologies; Example 37 predicts HQ,K ̸= 0 for Softmax self-attention due to its nonzero second derivative, but this has not been empirically verified. We

GapGN

100

init final

10−3

10−6

10−9

eLU

L +Si

sum

U ReL

cat+

sum

SiL

cat+

Figure 6: Exp. 4: GN-Gap at the merge node of Diamond MLP (k = 2). Logarithmic scale highlights the ∼7 order-of-magnitude separation between cat + SiLU (σ ′′ ̸= 0, nonlinear merge) and the remaining configurations (linear merge or σ ′′ = 0 a.e.). Summerge yields T = 0 by construction — Gap ≈ 0 regardless of σ.

design a minimal experiment contrasting a single-head Attention model with a per-position ReLU-MLP ′′ T ̸= 0 is scale-invariant: σSoftmax (z) ̸= 0 for all z ∈ Rd control of comparable size. The prediction HQ,K regardless of dimension or parameter values; hence a minimal architecture suffices for falsification, and the qualitative distinction (Gap ≫ 0 for Softmax vs. Gap ≈ 0 for piecewise-linear activations) is a structural invariant of the computation graph. Protocol. Single-head self-attention: d = 16, S = 8, projections WQ , WK , WV ∈ Rd×d (no bias), meanpool, linear head; 785 parameters. Control: per-position ReLU-MLP, 3×[Linear(d, d)+ReLU] shared across positions, mean-pool, linear head; 833 parameters. Data: synthetic regression y = σ̄(XWt ) · 1/S · wr + ε, where σ̄ = tanh, Wt ∈ Rd×d , wr ∈ Rd are fixed random teacher weights, ε ∼ N (0, 0.01); targets standardized to zero mean and unit variance. Training: Adam (η = 10−3 ), cosine schedule, 30 epochs, batch 128, 5 seeds. f GN , GN-Gap for pairs (Q, K) (attention) and (block0 , block1 ) (MLP); computed via exact Metrics: Hv,w , Hv,w autograd on 64 validation samples at checkpoints init, mid, final. Hypotheses. (H5.1) For Attention, GN-Gap for the pair (Q, K) is substantially nonzero at all check′′ T ̸= 0 due to σSoftmax ̸= 0). (H5.2) For per-position ReLU-MLP, GN-Gap ≈ 0 at machine precision points (HQ,K ′′ (Proposition 10: σ = 0 a.e.). (H5.3) Both models converge on the synthetic target (monotone loss decrease), ensuring that the GN-Gap comparison is not confounded by training failure. Results. All hypotheses confirmed (Table 7). (H5.1) Attention yields GN-Gap = 69.3 ± 20.1 at conf GN ∥F accounts for less than 2 % of ∥HQ,K ∥F , confirming that the Gauss–Newton approxvergence; ∥HQ,K imation is structurally inadequate for the (Q, K) cross-block (Example 37). (H5.2) ReLU-MLP yields Gap = 7.4·10−8 ±2.0·10−8 (numerical zero); ∥H T ∥F /∥H f ∥F < 10−7 at all checkpoints. The ratio of gaps exceeds 109 , surpassing the seven-order-of-magnitude separation of Exp. 4 (cat + SiLU vs. σ ′′ = 0). (H5.3) Both models converge (Attention: MSE 0.99 → 0.57; ReLU-MLP: 0.92 → 0.03); the GN-Gap for Attention remains large and stable across checkpoints (57.5 → 79.9 → 69.3), confirming that it reflects architectural structure rather than training state. Table 7: Exp. 5: GN-Gap for Toy-Attention vs. per-position ReLU-MLP (d = 16, S = 8, synthetic regression, mean ± 1σ over 5 seeds). ∥H T ∥F : Frobenius norm of the tensor component; Gap = ∥H T ∥F /∥H GN ∥F . Model

Checkpoint

∥H f ∥F

∥H GN ∥F

∥H T ∥F

Gap

Attention Attention Attention

init mid final

9.2·10−3 5.2·10−3 5.2·10−3

1.7·10−4 7.4·10−5 8.7·10−5

9.2·10−3 5.2·10−3 5.2·10−3

57.5 ± 16.4 79.9 ± 23.3 69.3 ± 20.1

ReLU-MLP ReLU-MLP ReLU-MLP

init mid final

3.5·10−2 3.5·10−1 3.6·10−1

3.4·10−9 2.8·10−8 2.6·10−8

9.6·10−8 8.8·10−8 7.4·10−8

5.6. Exp. 6: Convolutional architecture (ResNet-18, CIFAR-10) Motivation. Experiments 1–5 operate on fully connected or toy architectures (P ≤ 105 ) where exact Hessian blocks are computable. To verify the predictions of the formalism at practical scale we conduct an experiment on ResNet-18 (He et al., 2016) (P ≈ 11·106 parameters) with CIFAR-10. Protocol. For each activation σ ∈ {ReLU, SiLU} two architectures are trained: (i) SegmentedResNet18—a standard torchvision ResNet-18 with 5 measurement points (stem, layer1–layer4), and (ii) SegmentedPlainResNet18—the same parameterization with identity shortcuts removed (conv–BN–act chain; stridebased downsampling is preserved). Training: SGD (η = 0.1, momentum 0.9, weight decay 5·10−4 ), cosine schedule, 100 epochs, batch size 128, standard CIFAR-10 augmentation (random crop, horizontal flip). Checkpoints: init, mid (epoch 50), final (epoch 100). Metrics R, C, D are computed via the Hutchinson estimator (30 Rademacher probes, subsample of 32 examples); GN-Gap via StochasticGNGapEstimator (30 probes, common probe vector). Five seeds {42, . . . , 46}. Segmentation and linear classifier. The model is partitioned into n = 6 segments: seg0 (conv1–BN– σ–maxpool), seg1–4 (ResNet layer1–layer4), seg5 (avgpool–flatten–fc). Measurement points are the outputs fstem , . . . , flayer4 ; maximum distance dmax = 4. In particular, seg5 (the linear classifier) contains no nonlinear activations: logits = W · avgpool(flayer4 ) + b. As a consequence, for any pair (v, w) with d = dmax , the path from fw to the logits passes exclusively through the linear map seg5 , which has a principled implication for GN-Gap at d = dmax (Remark 39). Hypotheses. (H6.1) For piecewise-linear activation (ReLU): GapGN ≈ 0 at all distances (H T ≡ 0 a.e.). (H6.2) For smooth activation (SiLU): GapGN > 0 at d < dmax (σ ′′ ̸= 0). (H6.3) Skip connections stabilize R̄(d) (ResNet vs. Plain). Results. (H6.1) Strongly confirmed. For ReLU, GN-Gap < 2·10−6 at all distances and both architectures at the final checkpoint (Table 8); across all three checkpoints the bound GapGN < 10−5 holds uniformly (Table Appendix N.8). This is consistent with H T ≡ 0 a.e. for piecewise-linear activations (Proposition 10) and reproduces the Exp. 3 result on a convolutional architecture with 11M parameters. (H6.2) Confirmed with refinement. For SiLU, GN-Gap at the final checkpoint is 0.43 (ResNet, d = 0) and decreases with distance to 0.15 at d = 3 (Table 8), confirming the presence of the tensor component H T ̸= 0. However, at d = dmax = 4, GN-Gap collapses to ∼10−6 for all configurations, including SiLU (Remark 39). (H6.3) Confirmed. At initialization, ReLU ResNet yields stable R̄(d) (R(0)/R(4) ≈ 1), whereas Plain exhibits ∼627× decay; for SiLU the contrast is stronger: 25× (ResNet) vs. 47 000× (Plain) (Figure 7). ¯ Coupling C(d) decreases monotonically in all configurations; skip connections slow the decay (Table 9). Test accuracy: 88.3 % (ReLU ResNet), 87.4 % (ReLU Plain), 88.2 % (SiLU ResNet), 87.3 % (SiLU Plain)— sufficient for meaningful curvature diagnostics. Remark 39 (Linear head and GN-Gap at d = dmax ). For a pair (v, w) with d = dmax (here: stem, layer4), the path from fw to the logits passes only through the linear map seg5 : fw 7→ W avgpool(fw ) + b. Since ∂ 2logits/∂fw2 = 0, the Jacobian Jw = ∂logits/∂fw does not depend on fw , and the block Hessian becomes f GN Hv,w = Jv⊤ ∇2z ℓ Jw = Hv,w , T i.e. Hv,w = 0 exactly, regardless of the activation type in intermediate layers. The residual value

∥H T ∥F ∥H GN ∥F ∼ 10−6 is a float32 arithmetic artifact. This effect is specific to architectures with a linear classifier head (standard fc head) and does not arise with nonlinear heads (e.g., MLP head in ViT).

ReLU

SiLU 102 ResNet init ResNet final Plain init Plain final

10−1

ResNet init ResNet final Plain init Plain final

10−1

R̄(d)

101

10−4

10−3

10−5 0

10−7

Distance d

Figure 7: Exp. 6: Mean resonance R̄(d) for ResNet-18 (CIFAR-10, log-scale on y; ±1σ over 5 seeds). (a) ReLU: ResNet preserves R̄ at init (R(0)/R(4) ≈ 1); Plain decays by ∼627×. (b) SiLU: decay is stronger—Plain init 47 000×; ResNet 25×, slower than Plain (skip connections stabilize curvature).

Table 8: Exp. 6: GN-Gap by distance d for ResNet-18 (CIFAR-10, final checkpoint, mean ± 1σ over 5 seeds). Row d = 4: linear classifier (Remark 39).

ReLU

0 1 2 3 4†

SiLU

ResNet

Plain

−6

< 2·10 < 2·10−6 < 2·10−6 < 2·10−6 < 2·10−6

< 10 < 10−6 < 2·10−6 < 10−6 < 2·10−6

ResNet

Plain

0.43 ± 0.12 0.24 ± 0.07 0.22 ± 0.06 0.15 ± 0.04 < 2·10−6

0.21 ± 0.11 0.12 ± 0.04 0.06 ± 0.02 0.05 ± 0.02 < 2·10−6

† H T = 0 exactly (linear head); see Remark 39.

¯ Table 9: Exp. 6: mean resonance R̄(d) and coupling C(d) at initialization (ResNet-18, CIFAR-10, mean over 5 seeds).

¯ C(d)

R̄(d) d

ReLU Res.

0 1 2 3 4

−2

4.5·10 3.5·10−2 3.3·10−2 3.6·10−2 4.8·10−2

SiLU Plain −2

2.1·10 3.5·10−4 6.0·10−5 2.9·10−5 3.3·10−5

Res. −2

3.0·10 1.0·10−2 4.2·10−3 2.1·10−3 1.2·10−3

ReLU Plain −2

2.1·10 1.1·10−4 6.9·10−6 1.2·10−6 4.4·10−7

SiLU

Res.

Plain

Res.

Plain

1.00 0.95 0.92 0.89 0.82

1.00 0.93 0.83 0.71 0.58

1.00 0.56 0.31 0.17 0.08

1.00 0.24 0.06 0.01 < 0.01

5.7. Ablation and reproducibility Comparison of exact and stochastic (30 Rademacher probes) GN-Gap estimates shows discrepancy < 1 % for all smooth activations (Table Appendix N.6, AppendixpN), confirming the adequacy of the Hutchinson estimator. At fixed m = 30 the relative error scales as O( rank(A)/m) (Avron and Toledo, 2011); probecount scaling guidelines are given in Appendix N. Implementation: PyTorch 2.x; all experiments are feasible on a single GPU (CIFAR-10/100, width ≤ 256). Seeds, hyperparameters, and configurations are fixed; code is available at https://github.com/comiam/dag-hesse. 20

6. Discussion and limitations Scaling. Explicit materialization of H is practical for P ≲ 105 ; for larger networks the framework reduces to an O(P ) HVP operator whose cost is 2–4× that of a backward pass. Experiments have been verified on MLPs up to L = 32 (Table Appendix N.3), on a toy single-head Attention model (Exp. 5), and on ResNet-18 (∼11M parameters, Exp. 6). Exp. 6 confirms the key predictions of the formalism under stochastic curvature estimation at practical scale. Scaling to ImageNet-scale models, extension to Transformer architectures, and integration of the proposed metrics into online training diagnostics remain key directions for future work. Role of experiments. The experiments serve a falsification role: they verify theoretical predictions (decay, rank constraints, vanishing H T ) on exact Hessian blocks. Experiments 1–4 use fully-connected f MLPs, where exact computation of Hv,w is feasible at moderate P . Experiment 5 extends the empirical T scope to self-attention, confirming the prediction of Example 37 (HQ,K ̸= 0). Experiment 6 extends the validation to a convolutional architecture (ResNet-18) via stochastic estimation, reproducing the resonance stabilization and GN-Gap predictions at 11M-parameter scale. The formalism itself is architecture-agnostic (linear, convolutional, Attention — Examples 5–37). ReLU and GN-Gap. For piecewise-linear activations the tensor component vanishes in activation space (Proposition 10), and the GN-Gap metric is uninformative; non-convexity is concentrated on the boundaries of linear regions. GN-Gap is critically important for smooth activations (GELU, Swish, Softmax), where H T ̸= 0 everywhere. Static estimates. All theorems are stated for fixed θ; during re-training the metrics reflect both architecture and the geometry induced by learning. To isolate the architectural effect, measurements at initialization are preferred. Potential applications. The formalism yields several consequences that are actionable without additional theoretical development. (i) Banded preconditioner design: Corollary 29 quantifies the approximation error of block-diagonal optimizers; cross-blocks beyond distance ⌈log ε/ log(sρ)⌉ are negligible, providing a principled bandwidth selection for banded preconditioners (George et al., 2018). (ii) Bottleneck diagnostics: Proposition 31 bounds the coordination rank of inter-layer blocks by the bottleneck width du ; this constraint is invisible from the scalar loss and can inform width allocation in architecture design. (iii) Activation/optimizer compatibility: the GN-Gap quantifies how much curvature information the GGN discards for a given activation type, serving as a criterion for choosing between GGN-based and full-Hessian methods. (iv) Online curvature monitoring: all metrics are estimated stochastically in O(P ) via HVP (Corollary 21); abrupt changes in R, C, or D during training may indicate phase transitions in the loss landscape geometry. 7. Conclusion This paper presents an analytical formalism for block-wise curvature analysis in neural networks of arbitrary DAG architecture, comprising a recursive formula for inter-layer blocks (C1), the canonical decomposition H GN +H T (C2), and diagnostic metrics: inter-layer resonance R, geometric coupling C, stable rank D, and GN-Gap. The theoretical analysis establishes exponential resonance decay in vanilla networks, the stabilizing effect of skip connections on curvature, rank constraints of bottleneck layers, and vanishing of the tensor component of the input Hessian for piecewise-linear activations (C3). Curvature routing (C4) formalizes the propagation of curvature through the graph and reduces to an O(P ) HVP operator. Beyond these theoretical contributions, the empirical validation (Section 5) confirmed the key predictions of the theory: exponential resonance decay without skip connections, rank constraints of bottlenecks on T coupling, selective vanishing of GN-Gap for piecewise-linear activations, nonzero HQ,K in Softmax selfattention (Exp. 5), and reproducibility of these patterns on a convolutional architecture (ResNet-18, ∼11M parameters, Exp. 6). At initialization, the decay is guaranteed by Theorem 25 (λ1 < 0); after training, λ1 > 0 for all configurations (Table 2), yet batch-averaged resonance still decays—an empirical finding whose formal characterization remains open (Remark 38). Key directions for future work include scaling the diagnostic metrics to ImageNet-scale models and integrating them with second-order optimizers. 21

Declaration of competing interests The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. CRediT authorship contribution statement Maxim Bolshim: Conceptualization, Methodology, Software, Formal analysis, Writing – original draft, Visualization. Alexander Kugaevskikh: Validation, Writing – review & editing. Acknowledgements The authors thank the anonymous reviewers for their constructive feedback. Funding This work was supported by the state assignment (project FSER-2025-0004). Declaration of generative AI and AI-assisted technologies in the manuscript preparation process During the preparation of this work the authors used Claude Opus (Anthropic) for partial editing and typesetting of the manuscript, and Claude Sonnet (Anthropic) for partial development of the accompanying software code. After using these tools, the authors reviewed and edited the content as needed and take full responsibility for the content of the published article. Data availability The source code and experimental results supporting the findings of this study are openly available at https://github.com/comiam/dag-hesse (archived: https://doi.org/10.5281/zenodo.19545553). References Abreu, N., Vyas, N., Kakade, S., Morwani, D., 2025. The potential of second-order optimization for LLMs: A study with full Gauss–Newton. arXiv preprint arXiv:2510.09378 doi:10.48550/arXiv.2510.09378. preprint. Avron, H., Toledo, S., 2011. Randomized algorithms for estimating the trace of an implicit symmetric positive semi-definite matrix. Journal of the ACM 58, 1–34. doi:10.1145/1944345.1944349. Balduzzi, D., Frean, M., Leary, L., Lewis, J.P., Ma, K.W.D., McWilliams, B., 2017. The shattered gradients problem: If resnets are the answer, then what is the question?, in: Proceedings of the 34th International Conference on Machine Learning (ICML), pp. 342–350. doi:10.48550/arXiv.1702.08591. Ben-Israel, A., Greville, T.N.E., 2003. Generalized Inverses: Theory and Applications. 2nd ed., Springer, New York, NY, USA. doi:10.1007/b97366. Bolshim, M., Kugaevskikh, A., 2025. Local properties of neural networks through the lens of layer-wise Hessians. arXiv preprint arXiv:2510.17486 doi:10.48550/arXiv.2510.17486. preprint. Bolte, J., Pauwels, E., 2021. Conservative set-valued fields, automatic differentiation, stochastic gradient methods and deep learning. Mathematical Programming 188, 19–51. doi:10.1007/s10107-020-01501-5. Botev, A., Ritter, H., Barber, D., 2017. Practical Gauss–Newton optimisation for deep learning, in: Proceedings of the 34th International Conference on Machine Learning (ICML), pp. 557–565. doi:10.48550/arXiv.1706.03662. Clarke, F.H., 1990. Optimization and Nonsmooth Analysis. SIAM, Philadelphia, PA, USA. doi:10.1137/1.9781611971309. Cohen, J.M., Kaur, S., Li, Y., Kolter, J.Z., Talwalkar, A., 2021. Gradient descent on neural networks typically occurs at the edge of stability, in: Proceedings of the 9th International Conference on Learning Representations (ICLR). doi:10.48550/ arXiv.2103.00065. Dangel, F., Kunstner, F., Hennig, P., 2020. BackPACK: Packing more into backprop, in: Proceedings of the 8th International Conference on Learning Representations (ICLR). doi:10.48550/arXiv.1912.10985. Dangel, F., Tatzel, L., Hennig, P., 2022. ViViT: Curvature access through the generalized Gauss–Newton’s low-rank structure. Transactions on Machine Learning Research (TMLR) doi:10.48550/arXiv.2106.02624.

Dauphin, Y.N., Pascanu, R., Gulcehre, C., Cho, K., Ganguli, S., Bengio, Y., 2014. Identifying and attacking the saddle point problem in high-dimensional non-convex optimization, in: Advances in Neural Information Processing Systems (NeurIPS), pp. 2933–2941. doi:10.48550/arXiv.1406.2572. George, T., Laurent, C., Bouthillier, X., Ballas, N., Vincent, P., 2018. Fast approximate natural gradient descent in a Kroneckerfactored eigenbasis, in: Advances in Neural Information Processing Systems (NeurIPS), pp. 9550–9560. doi:10.48550/arXiv. 1806.03884. Ghorbani, B., Krishnan, S., Xiao, Y., 2019. An investigation into neural net optimization via Hessian eigenvalue density, in: Proceedings of the 36th International Conference on Machine Learning (ICML), pp. 2232–2241. doi:10.48550/arXiv.1901. 10159. Gupta, V., Koren, T., Singer, Y., 2018. Shampoo: Preconditioned stochastic tensor optimization, in: Proceedings of the 35th International Conference on Machine Learning (ICML), pp. 1842–1850. doi:10.48550/arXiv.1802.09568. Hanin, B., Rolnick, D., 2019. Complexity of linear regions in deep networks, in: Proceedings of the 36th International Conference on Machine Learning (ICML), pp. 2596–2604. doi:10.48550/arXiv.1901.09021. Hayou, S., Doucet, A., Rousseau, J., 2019. On the impact of the activation function on deep neural networks training, in: Proceedings of the 36th International Conference on Machine Learning (ICML), pp. 2672–2680. doi:10.48550/arXiv.1902. 06853. He, K., Zhang, X., Ren, S., Sun, J., 2016. Deep residual learning for image recognition, in: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pp. 770–778. doi:10.1109/CVPR.2016.90. Horn, R.A., Johnson, C.R., 2012. Matrix Analysis. 2nd ed., Cambridge University Press, Cambridge, UK. doi:10.1017/ CBO9781139020411. Jacot, A., Gabriel, F., Hongler, C., 2018. Neural tangent kernel: Convergence and generalization in neural networks, in: Advances in Neural Information Processing Systems (NeurIPS), pp. 8571–8580. doi:10.48550/arXiv.1806.07572. Jeyakumar, V., Luc, D.T., 1999. Nonsmooth calculus, minimality, and monotonicity of convexificators. Journal of Optimization Theory and Applications 101, 599–621. doi:10.1023/A:1021790120780. Kunstner, F., Hennig, P., Balles, L., 2019. Limitations of the empirical Fisher approximation for natural gradient descent, in: Advances in Neural Information Processing Systems (NeurIPS), pp. 4156–4167. doi:10.48550/arXiv.1905.12558. Liu, H., Li, Z., Hall, D., Liang, P., Ma, T., 2024. Sophia: A scalable stochastic second-order optimizer for language model pre-training, in: Proceedings of the 12th International Conference on Learning Representations (ICLR). doi:10.48550/ arXiv.2305.14342. Magnus, J.R., Neudecker, H., 2019. Matrix Differential Calculus with Applications in Statistics and Econometrics. 3rd ed., John Wiley & Sons, Chichester, UK. doi:10.1002/9781119541219. Martens, J., 2010. Deep learning via Hessian-free optimization, in: Proceedings of the 27th International Conference on Machine Learning (ICML), pp. 735–742. doi:10.48550/arXiv.1003.0358. Martens, J., 2020. New insights and perspectives on the natural gradient method. Journal of Machine Learning Research 21, 1–76. doi:10.48550/arXiv.1412.1193. Martens, J., Grosse, R., 2015. Optimizing neural networks with Kronecker-factored approximate curvature, in: Proceedings of the 32nd International Conference on Machine Learning (ICML), pp. 2408–2417. doi:10.48550/arXiv.1503.05671. Martens, J., Sutskever, I., 2012. Training deep and recurrent networks with Hessian-free optimization, in: Neural Networks: Tricks of the Trade. 2nd ed.. Springer, Berlin, Germany. volume 7700 of Lecture Notes in Computer Science, pp. 479–535. doi:10.1007/978-3-642-35289-8_27. Miyato, T., Kataoka, T., Koyama, M., Yoshida, Y., 2018. Spectral normalization for generative adversarial networks, in: Proceedings of the 6th International Conference on Learning Representations (ICLR). doi:10.48550/arXiv.1802.05957. Mordukhovich, B.S., 2006. Variational Analysis and Generalized Differentiation I: Basic Theory. Springer, Berlin, Germany. doi:10.1007/3-540-31247-1. Nocedal, J., Wright, S.J., 2006. Numerical Optimization. 2nd ed., Springer, New York, NY, USA. doi:10.1007/ 978-0-387-40065-5. Oseledets, V.I., 1968. A multiplicative ergodic theorem. Characteristic Ljapunov exponents of dynamical systems. Trudy Moskov. Mat. Obshch. 19, 179–210. Papyan, V., 2019. Measurements of three-level hierarchical structure in the outliers in the spectrum of deepnet Hessians, in: Proceedings of the 36th International Conference on Machine Learning (ICML), pp. 5012–5021. doi:10.48550/arXiv.1901. 08244. Pascanu, R., Mikolov, T., Bengio, Y., 2013. On the difficulty of training recurrent neural networks, in: Proceedings of the 30th International Conference on Machine Learning (ICML), pp. 1310–1318. doi:10.48550/arXiv.1211.5063. Pearlmutter, B.A., 1994. Fast exact multiplication by the Hessian. Neural Computation 6, 147–160. doi:10.1162/neco.1994. 6.1.147. Pennington, J., Worah, P., 2018. The spectrum of the Fisher information matrix of a single-hidden-layer neural network, in: Advances in Neural Information Processing Systems (NeurIPS), pp. 5410–5419. doi:10.48550/arXiv.1810.01075. Poole, B., Lahiri, S., Raghu, M., Sohl-Dickstein, J., Ganguli, S., 2016. Exponential expressivity in deep neural networks through transient chaos, in: Advances in Neural Information Processing Systems (NeurIPS), pp. 3360–3368. doi:10.48550/arXiv. 1606.05340. Rockafellar, R.T., Wets, R.J.B., 1998. Variational Analysis. Springer, Berlin, Germany. doi:10.1007/978-3-642-02431-3. Sagun, L., Evci, U., Güney, V.U., Dauphin, Y., Bottou, L., 2017. Empirical analysis of the Hessian of over-parametrized neural networks. arXiv preprint arXiv:1706.04454 doi:10.48550/arXiv.1706.04454. preprint. Schoenholz, S.S., Gilmer, J., Ganguli, S., Sohl-Dickstein, J., 2017. Deep information propagation, in: Proceedings of the 5th International Conference on Learning Representations (ICLR). doi:10.48550/arXiv.1611.01232.

Schraudolph, N.N., 2002. Fast curvature matrix–vector products for second-order gradient descent. Neural Computation 14, 1723–1738. doi:10.1162/08997660260028683. Serra, T., Tjandraatmadja, C., Ramalingam, S., 2018. Bounding and counting linear regions of deep neural networks, in: Proceedings of the 35th International Conference on Machine Learning (ICML), pp. 4558–4566. doi:10.48550/arXiv.1711. 02114. Singh, S.P., 2025. The geometry of learning via loss landscape curvature. Ph.D. thesis. ETH Zurich. Tang, Q.Y., Gu, Y., Cai, Y., Sun, M., Li, P., Zhou, X., Xie, Z., 2025. Investigating the overlooked Hessian structure: From CNNs to LLMs, in: Proceedings of the 42nd International Conference on Machine Learning (ICML), pp. 58805–58831. URL: https://proceedings.mlr.press/v267/tang25d.html. Yao, Z., Gholami, A., Shen, S., Mustafa, M., Keutzer, K., Mahoney, M.W., 2021. AdaHessian: An adaptive second order optimizer for machine learning, in: Proceedings of the AAAI Conference on Artificial Intelligence, pp. 10665–10673. doi:10. 1609/aaai.v35i12.17275. Zhang, H., Chen, W., Liu, T.Y., 2018. On the local Hessian in back-propagation, in: Advances in Neural Information Processing Systems (NeurIPS). URL: https://proceedings.neurips.cc/paper/2018/hash/ b6d67a24906e8a8541291882f81d31ca-Abstract.html.

Appendix A. Notation details and function spaces Remark 40 (Index convention). i — component index of the node output; j, k — component indices of inputs; k, ℓ (in parameter context) — component indices of θv ; v, w, u — node indices. Standard function spaces are employed: C 2 (twice continuously differentiable, smooth case), C 1,1 (functions with Lipschitz derivatives), P C 2 (piecewise C 2 with piece boundaries of measure zero). In the nonsmooth case we rely on the apparatus of generalized differentiation (Mordukhovich, 2006) and, in particular, Clarke subdifferentials (Clarke, 1990): for a locally Lipschitz function, ∂C f (x) is non-empty, convex, and 2 F (x) exists when ∇Fi is Lipschitz. compact; the Clarke Hessian ∂C Remark 41 (Tensor notation and contraction rules). Index i refers to the output component of a node (fu or fv ); j, m to inputs from parent nodes; α, β to parameters θv . Notation [T ]i,•,• denotes a tensor slice Pdu at fixed i. Contraction [Tu;v ]i,j,k δu,i yields a dv ×dv matrix with entries i=1 [Tu;v ]i,j,k δu,i . Under matrix ⊤ f multiplication Du←v Hu,u Du←w , the dimensions chain as (dv ×du )(du ×du )(du ×dw ) = dv ×dw . Rules for differentiating matrix expressions follow the standard formalism (Magnus and Neudecker, 2019). Element-wise definitions of second-derivative tensors (v, w ∈ Pa(u)): [Tu;v ]i,j,k =

∂ 2 (fu )i ∈ Rdu ×dv ×dv , ∂(fv )j ∂(fv )k

[Tu;v,w ]i,j,k =

∂ 2 (fu )i ∈ Rdu ×dv ×dw , ∂(fv )j ∂(fw )k

[Tv;w,θ ]i,j,k =

∂ 2 (fv )i ∈ Rdv ×dw ×pv , ∂(fw )j ∂(θv )k

[Tvθ ]i,k,ℓ =

∂ 2 (fv )i ∈ Rdv ×pv ×pv . ∂(θv )k ∂(θv )ℓ

Remark 42 (Symmetry of the tensor component). The mixed tensor Tu;v,w ∈ Rdu ×dv ×dw for v ̸= w produces rectangular dv × dw slices that cannot be symmetrized. Symmetry of the full Hessian is ensured by the f f relation Hv,w = (Hw,v )⊤ between blocks. For diagonal blocks (v = w) the slices are square; in the smooth ⊤ case [Tu;v ]i,:,: = [Tu;v ]i,:,: by Schwarz’s theorem. The double sum over (u1 , u2 ) ∈ Ch(v) × Ch(w) in term 1 entails no double-counting: each pair contributes independently. Assembling blocks into the global Hessian.. The global Hessian ∇2θ L ∈ RP ×P is assembled from the computed parametric blocks Hθvi ,θvj as a block matrix [Hθvi ,θvj ]ni,j=1 with symmetry by Schwarz’s theorem.

Appendix B. Stability and invariance of the curvature recursion Theorem 43 (Error accumulation in the curvature recursion). Under local Jacobian errors ∥∆D∥2 ≤ ϵD ∥D∥2 and tensor errors ∥∆T ∥F ≤ ϵT ∥T ∥F : εv,w ≤ dist(v, w)(2ϵD + ϵT ) + O(ϵ2 ). The protocol is stable whenever L(2ϵD + ϵT ) < 1. The Hutchinson estimator (Avron and Toledo, 2011) is used for trace and norm computation, converging p as ∥T̂ − T ∥F ≤ C∥T ∥F log(1/δ)/m with probability ≥ 1−δ using m probe vectors. Theorem 44 (Invariance of geometric coupling). Let v ∈ Pa(w) and consider the rescaling Wv 7→ αWv , f Ww 7→ α−1 Ww that preserves the output fw (and hence all downstream activations and L). Then: Hv,v 7→ −2 f f −1 f f α Hv,v , Hv,w 7→ α Hv,w , Hw,w is invariant, so that ˜ w) = q C(v,

f α−1 ∥Hv,w ∥F

= C(v, w).

f f α−2 ∥Hv,v ∥F · ∥Hw,w ∥F

In a DAG with branching (v has several children), rescaling Wv 7→ αWv affects all children: preserving outputs requires compensating transformations on every child, but the ratio C(v, w) remains invariant provided fw and all fu , u ∈ Ch(v) are preserved. Appendix C. Optimization-theoretic remarks Appendix C.1. Convergence guarantees Using the exact Hessian (or its AD analog for non-smooth networks) as Bk in second-order methods yields classical convergence guarantees: quadratic for Newton’s method and superlinear when ∥Bk −∇2 L∥ → 0 (Nocedal and Wright, 2006). In practice, regularization Bk + λk I and symmetrization 12 (Bk + Bk⊤ ) ensure numerical stability. Appendix C.2. Pseudoinverse and Newton direction In Definition 47 the AD-Hessian is defined as a valid element of a conservative field (via CSVF (Bolte and Pauwels, 2021)). Here we clarify its connection with the Moore–Penrose pseudoinverse. 2 Remark 45 (Relation to the pseudoinverse). For an arbitrary H ∈ ∂C L, the Newton direction d = −H † g † (g = ∇L, H is the Moore–Penrose pseudoinverse) is the minimum-norm solution of Hd = −g (Ben-Israel and Greville, 2003). At full rank (rank(H) = p), H † = H −1 and the direction coincides with the standard Newton step.

Appendix D. Non-smooth case: Clarke Hessian and CSVF This section contains the full theory of the non-smooth case omitted from the main text. Methodologically we rely on the non-smooth calculus of convexificators (Jeyakumar and Luc, 1999) and work within the Conservative Set-Valued Fields (CSVF) framework (Bolte and Pauwels, 2021): for a piecewise-analytic 2 network F , the mapping x 7→ ∂C L(x) is a conservative field whose potential is ∇L (a.e.). This circumvents the well-known intransitivity of the chain rule for the Clarke subdifferential (Rockafellar and Wets, 1998): instead of an element-wise chain rule, the global conservativity property guarantees correctness of AD computations for the entire computation graph. Theorem 46 (AD-Hessian of a ReLU DAG network). Let F be a ReLU DAG network and ℓ ∈ C 2 . Then: 1. Each fv is locally Lipschitz. 25

2. Under a transversality condition on the activation maps (the mapping x 7→ [fu1 (x), . . . ]⊤ intersects non-smoothness hyperplanes transversally), the set of non-smooth points has Lebesgue measure zero. 3. At points of differentiability, the AD-Hessian coincides with the classical one. 4. At non-smooth points, AD frameworks return an AD-Hessian (with T ≡ 0), which amounts to selecting the zero element from the generalized derivative and preserves algorithmic convergence by CSVF theory (Bolte and Pauwels, 2021). Definition 47 (AD-Hessian at a non-smooth point). For ReLU networks, AD frameworks set σ ′′ (0) := 0, so that Tu;v ≡ 0. The input Hessian block at a non-smooth point is defined through Jacobian products Du←v (computed via the convention σ ′ (0) := 0). This choice is valid under CSVF theory (Bolte and Pauwels, 2021): the AD-Hessian is a valid element of a conservative field, guaranteeing convergence of second-order optimization algorithms. Theorem 48 (Stability of the AD-Hessian). Let xn → x0 with the network twice differentiable at every xn . 2 Then ∇2 L(xn ) → M implies M ∈ ∂C L(x0 ). At smooth points the AD-Hessian is continuous in x. Proof of Theorem 46. Step 1 (local Lipschitz continuity). Consider the topological ordering v1 , . . . , v|V | . For the input node, fv1 (x) = x is 1-Lipschitz. If node v receives outputs fu1 , . . . , fuk and applies Wv ( · ) + bv followed by ReLU: fv (x) = ReLU Wv [fu1 (x), . . . , fuk (x)]⊤ + bv . Qk The linear map has Lipschitz constant ∥Wv ∥2 ; ReLU has constant 1. Hence fv is i=1 Lui ∥Wv ∥2 -Lipschitz, where Lui is the Lipschitz constant of fui . Induction over the topological order yields local Lipschitz continuity for all fv . Step 2 (measure of the smooth set). L is non-smooth on submanifolds corresponding to boundaries of ReLU linear regions. For each neuron the set of points where the pre-activation equals zero is given by the equation Wv [fu1 (x), . . . , fuk (x)]⊤ +bv = 0. Under the rank condition (the mapping x 7→ [fu1 (x), . . . , fuk (x)]⊤ has locally full rank), this is a codimension-1 hypersurface of measure zero. Geometrically, this is a transver sality condition: the activation sign map x 7→ sign(preactj (x)) j intersects the non-smoothness hyperplanes {zj = 0} transversally, which guarantees codimension ≥ 1 for each neuron (Hanin and Rolnick, 2019). When transversality is violated (degenerate Jacobians, dead ReLU, bottleneck layers), the rank condition may fail; for F (x) = max{0, x}, L(x) = F (x)2 /2 the function is twice non-differentiable on (−∞, 0]. When the rank condition holds, the non-smoothness set is covered by at most 2N codimension-1 subspaces (Hanin and Rolnick, 2019; Serra et al., 2018). 2 Step 3 (coincidence at smooth points). At a twice-differentiable point, ∂C L(x) = ∇2 L(x) . Inside a linear region F is affine (F (x) = Ax + b), so ∇2 L(x) = A⊤ ∇2 ℓ(F (x))A. Step 4 (AD-Hessian at non-smooth points). For piecewise-linear activations (ReLU), σ ′′ (z) = δ(z) is a distribution, not an element of L∞ ; the classical Clarke apparatus is not directly applicable. However, AD frameworks assign σ ′′ (0) := 0, zeroing out Tu;v ≡ 0 for all ReLU nodes. The AD-Hessian is therefore defined solely through Jacobian products Du←v . Under CSVF theory (Bolte and Pauwels, 2021), the AD gradient g(x) := AD[∇L](x) is an element of a conservative field for piecewise-analytic functions. Key observation: with σ ′′ (0):=0, the mapping g is piecewise-affine (off a measure-zero set it equals x 7→ A⊤ r ∇ℓ(Ar x+br ) for linear region r). Since g is piecewiseanalytic, its Jacobian Jg (x) = AD[∇2 L](x) is defined a.e. and is itself an element of a conservative field ∂C g (closure by the CSVF chain rule, Proposition 3 in (Bolte and Pauwels, 2021)). Thus the AD-Hessian satisfies the conservative-field axioms and is valid for second-order optimization.

Appendix E. Implementation via autodiff frameworks Algorithm 1 Computing an input Hessian block 1: function ComputeInputHess(v, w, {fu }, L, {δu }, {H f }, C) f 2: if (v, w) ∈ C then return Hv,w 3: end if f 4: Hv,w ←0 f 5: if v, w influence L then Hv,w ← ∂ 2 L/∂fv ∂fw 6: end if 7: for u1 ∈ Ch(v) do 8: Dv ← jac(fu1 , fv ) 9: for u2 ∈ Ch(w) do 10: Dw ← jac(fu2 , fw ) 11: if (u1 , u2 ) ∈ / C then Huf1 ,u2 ← ComputeInputHess(u1 , u2 , . . . ) 12: end if f f 13: Hv,w ← Hv,w + Dv⊤ Huf1 ,u2 Dw 14: end for 15: end for 16: if v ̸= w then 17: for u ∈ Ch(v) ∩ Ch(w) do 18: for i ∈ 1..du do 19: T ← MixedHess(fu,i , fv , fw ) f f 20: Hv,w ← Hv,w + T · δu,i 21: end for 22: end for 23: else 24: for u ∈ Ch(v) do f f 25: for i do Hv,v ← Hv,v + hess(fu,i , fv ) · δu,i 26: end for 27: end for 28: end if f 29: C ← C ∪ {(v, w)}; return Hv,w 30: end function

Algorithm 2 Computing a parametric Hessian block 1: function ComputeParamHess(v, w, {fu }, {θ}, {H f }, {δ}) 2: Dv ← jac(fv , θv ); Dw ← jac(fw , θw ) f 3: Hθv ,θw ← Dv⊤ Hv,w Dw 4: if v = w then 5: for i do Hθv ,θv ← Hθv ,θv + hess(fv,i , θv ) · δv,i 6: end for 7: end if 8: for u ∈ Pa(v) ∩ Ch(w) do 9: for i, j, α do 10: T ← MixedDeriv(fv,i , fu,j , θv,α ) 11: Hθv ,θw ← Hθv ,θw + T · jac(fu , fw ) · δv,i 12: end for 13: end for 14: return Hθv ,θw 15: end function

▷ Term (1): double sum over Ch(v) × Ch(w)

▷ Term (2): mixed tensor (v ̸= w)

▷ Term (3): pure tensor (v = w)

Algorithm 3 Full Hessian computation 1: function FullHessian(G, {fv }, {θv }, L) f 2: δout ← ∇L; Hout ← ∇2 L 3: topo ← TopoSort(G).rev(); C ← ∅ f ← hess(L, fv ); C ← C ∪ {(v, v)} 4: for v ∈ InputDepNodes(L) do Hv,v 5: end for f 6: for v, w ∈ InputDep, v ̸= w do Hv,w ← mixed_hess(L, fv , fw ); C ← C ∪ {(v, w)} 7: end for 8: for v ∈ topo do 9: BackpropGrad(v) 10: for w : Ch(v) ̸= ∅ OR Ch(w) ̸= ∅ OR (v, w) influences L do f 11: Hv,w ← ComputeInputHess(v, w, . . . , C) 12: end for 13: end for 14: Init H of size P × P 15: for v, w ∈ V : ∃ path v →∗ u ←∗ w OR (v, w) influences L do 16: Hθv ,θw ← ComputeParamHess(v, w, . . . ) 17: Update H 18: end for 19: return H 20: end function

Appendix F. Proof of the canonical decomposition We give the full proof of Theorem 7. GN GN = ∇2 L and rule by a recursion with base Hout,out Proof. Define Hv,w GN Hv,w =

Du⊤1 ←v HuGN Du2 ←w + 1 ,u2

u1 ∈Ch(v) u2 ∈Ch(w)

∂2L . ∂fv ∂fw

The recursion is closed: at each step only HuGN are used, excluding tensor terms. 1 ,u2 ⊤ GN ∇2 L Dout←w . This coincides = Dout←v Unrolling to the output node and applying induction: Hv,w ⊤ 2 exactly with the GGN matrix J ∇ L J (cf. Remark 9). T The tensor component H T := H f − H GN inherits the recursion with base Hout,out = 0; substituting H f = H GN + H T into equation 2 and subtracting the H GN recursion yields only the tensor terms Tu;v , Tu;v,w weighted by δu,i . Appendix F.1. Full parametric Hessian formula Let Dv =

∂fv , ∂θv

(i) Tu;v,α =

∂ 2 fu,i , ∂fv ∂θv,α

(i)

Tv,αβ =

∂ 2 fv,i . ∂θv,α ∂θv,β

The full parametric block obtained from equation 2: f Hθv ,θw = Dv⊤ Hv,w Dw | {z } GN-like X X (i) + δv,i Tv(i) + δu,i Dv⊤ Tu;v Dw i

u,i

}

(i) δu,i Tu;v,w ,

tensor over θv (v=w)

u,i

}

mixed tensor (v̸=w, u∈Pa(v)∩Ch(w))

}

cross-tensor (shared params)

(Appendix F.1)

where the fourth term appears only with weight sharing (θv ≡ θw across nodes); for standard architectures without weight sharing it vanishes. Remark 49 (Weight sharing between nodes). If parameter P θ is shared among nodes Vθ = {v1 , . . . , vm }, the m corresponding diagonal block of the full Hessian is Hθ,θ = i,j=1 Hθvi ,θvj (direct consequence of additivity of second derivatives over shared parameter copies). Appendix G. Diagnostic metrics and architectural analysis Appendix G.1. Inter-layer resonance f The block structure of the inter-layer Hessian {Hv,w } allows introducing a quantitative measure of “geometric connectivity” between nodes. The metrics R, C, and D are defined in the main text (Definitions 23 and 24); below we develop additional properties and proofs. From Definition 23 and formula equation 2: R(v, w) = 0 iff the nodes share no path to the output; R(v, w) = R(w, v) (Hessian symmetry); R(v, w) ≥ 0 (norm property).

Remark 50 (Meaning of inter-layer resonance). The curvature recursion induces a weighted interaction graph on the DAG with edge weights R(v, w). Unlike Adam and K-FAC, which treat layers in isolation, resonance allows visualizing curvature “highways,” explaining the effectiveness of skip connections (high R between distant layers), and diagnosing architectural bottlenecks. Appendix G.2. Geometric coupling and layer coherence Geometric coupling C(v, w) (Definition 23) normalizes resonance by the diagonal blocks, yielding a scaleinvariant connectivity measure. f ]v,w ⪰ 0, then Properties: (i) C ≥ 0; C = 0 iff R = 0. (ii) PSD bound: if H T = 0 and [Hv,w C ≤ 1 (Cauchy–Schwarz for block PSD matrices). (iii) Scale invariance under rescaling Wv 7→ αWv , Ww 7→ α−1 Ww (Theorem 44). In general, C > 1 is possible only due to the tensor component H T and indicates departure from the GN PSD regime. GN ]v,w ⪰ 0 admits a decomposition Proof of the PSD bound. With H T = 0, the block matrix M := [Hv,w P M = LL⊤ . Writing L = Q (P ∈ Rdv ×r , Q ∈ Rdw ×r ): Mv,v = P P ⊤ , Mv,w = P Q⊤ , Mw,w = QQ⊤ . Set S := P ⊤P , T := Q⊤Q (S, T ⪰ 0). Then ∥Mv,w ∥2F = tr(ST ), ∥Mv,v ∥2F = tr(S 2 ), ∥Mw,w ∥2F = tr(T 2 ). By the Cauchy–Schwarz inequality for the Hilbert–Schmidt inner product: (tr(ST ))2 ≤ tr(S 2 ) tr(T 2 ), whence C(v, w) ≤ 1.

Remark 51 (Stochastic estimation and Cˆ > 1 artifacts). Under stochastic estimation the norms ∥Hv,w ∥F , ∥Hv,v ∥F , ∥Hw,w ∥F are computed by independent Hutchinson probe series. The ratio of three independent ˆ 1). estimates need not respect the theoretical bound C ≤ 1; in practice the result is projected: Cˆ := min(C, With increasing m the violation probability decreases as O(1/m). Remark 52 (Interpretation of coupling). C ≈ 1: layers function as a single module. C ∈ (0.3, 0.7): moderate connectivity, typical for adjacent layers. C ≈ 0: geometric isolation, layers can train independently. Remark 53 (Self-compensation under sequential narrowing). In an MLP, a bottleneck with du ≪ dv constrains GN GN the rank of both Hv,w and Hv,v . The normalization in C admits partial compensation: numerator and denominator are bounded by the same narrow layer. Hence C diagnoses relative geometric alignment but is insensitive to absolute rank compression. f Remark 54 (Computing R(v, w) without the full matrix). The block Hv,w ∈ Rdv ×dw is a compact matrix f (typically dv , dw ≪ P ), so R(v, w) = ∥Hv,w ∥F is computed in O(dv dw ) without building the full P × P Hessian. When even the block cannot be stored, the Hutchinson estimator (Avron and Toledo, 2011) is used: ∥A∥2F = tr(A⊤A) = Ez [z ⊤A⊤Az], where each Az is one HVP (O(P ) time).

Appendix G.2.1. Rank bottleneck Proposition 55 (Rank bound on path contributions). Let v, w ∈ V and let c be a common descendant. Consider the path contribution to the GN component along paths pv : v → c and pw : w → c: Γ(pv , pw ) = Dp⊤v ∇2Lfc Dpw . Then P rank(Γ) ≤Pminu∈pv ∪pw rank(Ju←pred(u) ). For sums, only subadditivity holds: rank( i Γi ) ≤ i rank(Γi ). Proof. By sub-multiplicativity of rank (Horn and Johnson, 2012): rank(Dpv ) ≤ minu∈pv rank(Ju←pred(u) ). Since Γ = Dp⊤v A Dpw : rank(Γ) ≤ min{rank(Dpv ), rank(Dpw )}. Corollary 56 (Architectural bottleneck). If all paths from v and w to out pass through a narrow layer u GN f with du ≪ min(dv , dw ), then rank(Hv,w ) ≤ du . Consequently, the full Hessian Hv,w is low-rank up to the T tensor perturbation Hv,w : layers v and w can coordinate training in at most du independent directions in the Gauss–Newton regime. GN Proof. All paths factor through u: Dout←v = Dout←u Mv , Dout←w = Dout←u Mw , so Hv,w = Mv⊤ (· · · )Mw du ×du with the inner matrix in R .

Definition 57 (Effective interaction dimension). deff (v, w) :=

f ∥∗ ∥Hv,w f ∥2 ∥Hv,w

(Appendix G.1)

where ∥ · ∥∗ is the nuclear norm, ∥ · ∥2 the spectral norm. f ). (ii) deff = 1 iff rank = 1. (iii) Equality Proposition 58 (Properties of deff ). (i) 1 ≤ deff ≤ rank(Hv,w with rank iff all nonzero singular values are equal. GN Remark 59 (Rank vs. norm). Proposition 31 bounds the rank of Hv,w , not the Frobenius norm. Even at rank ≤ du , R(v, w) = ∥H∥F can remain large (singular values concentrate on du directions). The correct bottleneck test is the stable rank D, not R.

Remark 60 (Double rank bound). For K-class softmax-CE loss, rank(∇2L) ≤ K −1. Together with Proposition 31: D(v, w) ≤ min(du , K −1). Remark 61 (Architectural consequences). The rank-propagation theorem explains the success of skip connections (bypassing bottlenecks via high-rank alternative paths), rank degradation in deep MLPs without skip connections, and the minimum hidden width required for training coordination. With paths bypassing u, P rank( p Dp⊤ (· · · )Dp ) can exceed du ; skip connections thus destroy “information bottlenecks” preserving full-rank Hessians. Appendix G.2.2. Stochastically computable stable rank P Computing deff requires the nuclear norm (∥H∥∗ = i σi ), equivalent to a full SVD. For scalable analysis we use the stable rank D(v, w) (Definition 24)—an alternative effective-dimension measure admitting stochastic estimation via HVP without materializing the full matrix. f Proposition 62 (Properties of stable rank). (i) 1 ≤ D ≤ rank(Hv,w ). (ii) D = 1 iff rank 1. (iii) D = rank √ iff all nonzero singular values are equal. (iv) D ≤ deff ≤ rank · D. P 2 Proof. Let and ai := σi /σ1 ∈ (0, 1]. Then D = ai , P σ1 ≥ · · · ≥ σr > 0 be the nonzero singular values deff = √ ai . (i)–(iii) follow from a1 = 1, ai ≤ 1. (iv) a2i ≤ ai gives D ≤ deff ; Cauchy–Schwarz gives deff ≤ r D.

Remark 63 (Stochastic estimation of D). ∥H∥2F is estimated by the Hutchinson method (Avron and Toledo, 2 d = 1 Pm ∥Hzk ∥2 with Rademacher vectors zk ∈ {−1, +1}dw . ∥H∥2 = σ 2 is estimated via 2011): ∥H∥ F

k=1

power iteration on H ⊤H; after T iterations the relative error is O((σ2 /σ1 )2T ) (Nocedal and Wright, 2006). Total cost: O((m + 2T ) · Backprop), memory O(P ). Algorithm 4 Stochastic estimation of D(v, w) via HVP

Require: Nodes v, w; number of probes m; iterations T ; operator HVPv,w Ensure: Estimate D̂(v, w) f 1: // Power iteration: σ̂12 ≈ ∥Hv,w ∥22 2: q ← randn(dw ); q ← q/∥q∥ 3: for t = 1, . . . , T do 4: p ← HVPv,w (q) 5: q ← HVP⊤ v,w (p) 6: q ← q/∥q∥ 7: end for 8: σ̂12 ← ∥HVPv,w (q)∥2 9: if σ̂12 < ε then return 0 10: end if 11: // Hutchinson: 12: S ← 0 13: for k = 1, . . . , m do 14: Sample zk ∈ {−1, +1}dw (Rademacher) 15: S ← S + ∥HVPv,w (zk )∥2 16: end for 17: return S/(m σ̂12 )

Appendix G.3. Path-based curvature analysis Appendix G.3.1. Path decomposition Definition 64 (DAG path and path Jacobian). For nodes v, w ∈ V , denote by P(v → c) the set of all directed paths from v to c ∈ Desc(v). For a path p = (u0 , u1 , . . . , uk ), the path Jacobian is Dp :=

k−1 Y

Dui+1 ←ui ∈ Rduk ×du0 .

(Appendix G.2)

i=0

Theorem 65 (Path decomposition of the Hessian). X X f Hv,w =

Infl(pv , pw ; c),

c∈Desc(v)∩Desc(w) pv ∈P(v→c) pw ∈P(w→c)

where Infl(pv , pw ; c) := Dp⊤v HfLc ,fc Dpw + (tensor terms). P Proof. From unrolling equation 2. The full Jacobian Dc←v = p∈P(v→c) Dp by linearity of differentiation. P P P P P f ⊤ L ⊤ L = Substituting: Hv,w c( pv Dpv ) Hfc ( pw Dpw ) + . . . = c pv ,pw Dpv Hfc Dpw + . . . Tensor terms decompose analogously. Definition 66 (Path resonance). R(pv , pw ) := ∥Infl(pv , pw ; c)∥F . P P Corollary 67 (Additivity of resonance over paths). R(v, w) ≤ c pv ,pw R(pv , pw ). Equality holds when all path contributions are sign-aligned. Appendix G.3.2. Theoretical justification of skip connections Definition 68 (Residual architecture). A skip connection between nodes u ≺ v adds an edge (u, v) with an identity (or linear) map: fv ← fv + Wskip fu . 31

Theorem 69 (Geometric stabilization via skip connections). Let Gbase be the base architecture, Gres the T GN one with added identity skip connections, and let L be convex. Assume ∥Hv,w ∥F ≤ η∥Hv,w ∥F with η < 1. ⊤ 2 ⊤ 2 ⊤ 2 Define B = ∆J ∇ L Jbase + Jbase∇ L ∆J + ∆J ∇ L ∆J. If B ⪰ 0, then for all v, w: (Appendix G.3)

Rres (v, w) ≥ (1 − η) Rbase (v, w). In the PSD regime (H T = 0): Rres ≥ Rbase .

Proof. With convex L, HGN = J ⊤∇2L J ⪰ 0. Adding identity skip edges expands the path set. HGN [Gres ] = HGN [Gbase ] + B. By B ⪰ 0: ∥HGN [Gres ]∥F ≥ ∥HGN [Gbase ]∥F . For the full R = ∥HGN + HT ∥F : by the reverse triangle inequality and the smallness condition on H T : Rres ≥ (1 − η)Rbase . Remark 70 (Global vs. block-wise monotonicity). The inequality is established for the full block matrix GN (global resonance). For individual blocks Hv,w (v ̸= w, not necessarily PSD), analogous monotonicity is not guaranteed: skip connections may redistribute curvature energy between blocks. However, on average over the graph the total Frobenius norm does not decrease. Theorem 71 (Resonance decay: vanilla vs. Pre-Activation ResNet). Consider an L-layer network with contracting activations (∥Dvi+1 ←vi ∥2 ≤ ρ < 1). (a) Vanilla (upper bound): Rvanilla (v0 , vL ) ≤ C · ρL → 0. (b) ResNet with identity skips every k layers (upper bound, independent of L): RResNet (v0 , vL ) ≤ C · ρk (does not worsen as L → ∞). (b′ ) ResNet (lower bound, PSD regime): If H T = 0 and ∇2L ⪰ 0: (Appendix G.4)

RResNet (v0 , vL ) ≥ (1 − ρk )L/k ∥∇2L∥F . k

The lower bound decays at rate ρk (compared with ln ρ1 for the vanilla network — exponentially slower when ρk ≪ 1). With residual branch scaling ∥DFi ∥2 ≤ αρ, α = c/L: (1 − (αρ)k )L/k → 1, i.e. the lower bound is independent of L. Proof. (a) Single path ∥Dp ∥2 ≤ ρL . (b) Each identity skip adds a path with Jacobian I. The shortest path from v0 to vL via skips consists of ⌈L/k⌉ hops; each hop traverses at most k base edges (≤ ρk ) plus a skip (∥D∥ = 1). By the triangle inequality the upper bound depends on k only. (b′ ) In Pre-Activation ResNet, the Jacobian of block j is I +Dres,j with ∥Dres,j ∥2 ≤ ρk . The full Jacobian QL/k factors as Dout←v0 = j=1 (I + Dres,j ). By multiplicativity of the minimum singular value: L/k

σmin (Dout←v0 ) ≥

L/k

σmin (I + Dres,j ) ≥

j=1

(1 − ∥Dres,j ∥2 ) ≥ (1 − ρk )L/k ,

j=1

where σmin (I +A) ≥ 1 − ∥A∥2 follows from ∥(I +A)x∥ ≥ ∥x∥ − ∥Ax∥ ≥ (1−∥A∥2 )∥x∥. ⊤ In the PSD regime: R(v0 , vL ) = ∥Dout←v ∇2L∥F ≥ (1 − ρk )L/k ∥∇2L∥F . 0 Remark 72 (Cross-path terms and the lower bound). The path decomposition (Theorem 65) writes the full P Jacobian as Dout←v0 = p∈P Dp , where P is the set of paths from v0 to out. Substituting into the GN block gives X ⊤ X HvGN = Dp ∇2L = Dp⊤ ∇2L, 0 ,out p

and individual path contributions Dp⊤ ∇2L may partially cancel: the inequality ∥ ′

⊤ ⊤ p Dp A∥F ≥ ∥Did A∥F The proof of (b ) circumvents this difficulty by applying σmin -multiplicativity to

does not hold in general. Q the factored product j (I + Dres,j ) rather than to the element-wise sum of paths, thereby handling crosspath interactions automatically through the minimum singular value of the full Jacobian without the need to control their signs. 32

Corollary 73 (Criterion of geometric stability). An architecture is geometrically stable (resonance does not decay with depth) iff for every pair of distant layers (v, w) there exists a path p with ∥Dp ∥2 ≥ ϵ > 0 independent of dist(v, w). Remark 74 (Architectural consequences of path analysis). (1) ResNet: ρ2 -stability via skipsP every 2 layers (verified in Exp. 1 and Exp. 6). (2) DenseNet: exponential number of paths maximizes R(pv , pw ). T (3) Transformer: attention creates “dynamic” skip connections (Exp. 5 verifies HQ,K ̸= 0). (4) U-Net: encoder–decoder skips maintain high C. Predictions (2) and (4) are structural consequences of Theorem 65 and remain to be validated empirically on the respective architectures. Proposition 75 (Optimal skip placement). To minimize resonance loss in an L-layer network with a budget of B additional edges, the optimal strategy is uniform skip placement with step k = ⌊L/B⌋, ensuring minv,w R(v, w) ≥ C · ρL/B . Proof sketch. By Theorem 77, resonance between layers v, w decays as (sρ)dist(v,w) . A skip edge of span k ensures that no pair of layers lies at graph distance exceeding k. Given B edges, the worst-case distance is minimized by uniform placement at intervals of k = ⌊L/B⌋, yielding minv,w R(v, w) ≥ C · ρ⌊L/B⌋ . Any non-uniform arrangement leaves a contiguous gap of length > k, strictly increasing the worst-case distance and reducing the minimum resonance. Appendix H. Additional theoretical results Appendix H.1. Structural sparsity of the Hessian f decay exponentially A key property of the inter-layer Hessian is its approximate sparsity: blocks Hv,w with the graph distance between nodes. Remark 76 (Computing blocks for a one-directional path). If a directed path v →∗ w exists in the DAG (i.e. f v is an ancestor of w), the block Hv,w satisfies f Hv,w =

⊤ f Du←v Hu,w .

u∈Ch(v)

This follows from the chain rule w.r.t. the first argument: every child u ∈ Ch(v) is closer to w in topological order, closing the recursion. Theorem 77 (Exponential curvature sensitivity (Hessian vanishing/exploding)). Let G = (V, E) be the DAG of a neural network with maximum out-degree s = maxv |Ch(v)|. Assume Jacobian norms are bounded: GN ∥Du←v ∥2 ≤ ρ for all edges (v, u) ∈ E, with sρ < 1. Then for the GN-resonance RGN (v, w) := ∥Hv,w ∥F : RGN (v, w) ≤ C GN · (sρ)dist(v,w) ,

(Appendix H.1)

GN where dist(v, w) is the shortest-path length in the undirected graph and C GN = ∥Hout,out ∥F = ∥∇2L∥F . For T piecewise-linear activations (ReLU, Leaky ReLU) the tensor component Hv,w ≡ 0 a.e. (Theorem 46), and f the bound extends to the full input resonance R(v, w) = ∥Hv,w ∥F .

Proof. Denote by ℓv the length of the longest directed path from v to out (topological depth to output). The proof proceeds by strong induction on ℓv + ℓw . Base case: ℓv + ℓw = 0 means v = w = out; RGN (out, out) = C GN = C GN · (sρ)0 . GN GN Inductive step: Suppose ℓv + ℓw ≥ 1. WLOG ℓv ≥ ℓw (otherwise use ∥Hv,w ∥F = ∥Hw,v ∥F ). Then ℓv ≥ 1, v ̸= out, and Ch(v) ̸= ∅. P GN ⊤ GN The GN recursion is closed (Theorem 7): Hv,w = u∈Ch(v) Du←v Hu,w for arbitrary (v, w), without restrictions on activation type. Bounding norms: GN GN GN ∥Hv,w ∥F ≤ |Ch(v)| · ρ · max ∥Hu,w ∥F ≤ sρ · max ∥Hu,w ∥F . u∈Ch(v)

u∈Ch(v)

For each u ∈ Ch(v) we have ℓu ≤ ℓv − 1 (the path v → u → · · · → out is one edge longer than from u), so GN ℓu + ℓw < ℓv + ℓw and by the inductive hypothesis ∥Hu,w ∥F ≤ C GN (sρ)dist(u,w) . By the triangle inequality in the undirected graph, dist(u, w) ≥ dist(v, w) − 1 (since u is a neighbor of v, i.e. dist(v, u) = 1). With sρ < 1: sρ · C GN (sρ)dist(u,w) ≤ sρ · C GN (sρ)dist(v,w)−1 = C GN (sρ)dist(v,w) .

Corollary 78 (Block-banded structure). For deep networks (L ≫ 1) with sρ < 1, the full Hessian ∇2θ L GN is approximately block-banded: GN-blocks Hv,w with dist(v, w) > k ∗ are negligible, where the threshold k ∗ k∗ is determined by (sρ) < ε for a given accuracy ε. For piecewise-linear activations (Proposition 10) this f extends to the full input blocks Hv,w . Definition 79 (Interaction radius). The interaction radius kε for a given accuracy ε > 0 is the minimum integer k such that f C · (sρ)k < ε · ∥Hv,v ∥F for all v ∈ V. f Blocks Hv,w with dist(v, w) > kε are set to zero, reducing the complexity from O(n2 ) to O(n · kε ) blocks.

Remark 80 (Practical significance). Theorem 77 explains the empirical observation of block-diagonal Hessian structure in deep networks and justifies “block-diagonal Hessian” approximations used in K-FAC and related methods. Remark 81 (Mathematical justification of the “edge of chaos”). Theorem 77 mathematically proves the “edge of chaos” phenomenon for curvature: vanilla networks without residual connections are structurally unstable and require strict balancing ρ = 1 (empirically achieved via He init (Hanin and Rolnick, 2019) and BatchNorm) to prevent exponential decay or explosion of the Hessian. While initialization can achieve ρ ≈ 1, during training of vanilla networks the spectral norm of layers inevitably deviates from 1. The theorem shows that vanilla networks are structurally unstable to such deviations (requiring infinitely precise balancing ρ = 1): • for ρ slightly below 1, the curvature signal decays exponentially: R(v, w) → 0; • for ρ slightly above 1 or s > 1 (branching), path contributions explode (gradient/curvature explosion). This agrees with mean-field theory (Schoenholz et al., 2017; Poole et al., 2016), where ρ = 1 corresponds to the critical phase-transition point. Theorem 71 rigorously shows that skip connections (ResNet) break this rigid requirement, guaranteeing R(v0 , vL ) ≥ C > 0 independently of L. Balduzzi et al. (2017) showed analogous O(2−L ) decay of gradient correlations. Appendix H.2. Negative curvature and saddle points The full Hessian, unlike the Gauss–Newton approximation, can detect saddle points through analysis of negative eigenvalues. Definition 82 (Negative curvature mass). For a symmetric matrix H with eigenvalues {λi }ni=1 , define the negative curvature mass: X m− (H) := |λi |. (Appendix H.2) λi <0

Theorem 83 (GN insensitivity to saddle points). Let L be a convex loss (e.g. MSE or cross-entropy). Then: 1. The Gauss–Newton component H GN is always positive semi-definite: m− (H GN ) = 0. 2. The full Hessian H f ull = H GN + H T may have m− (H f ull ) > 0, indicating a saddle point. 3. Negative curvature arises exclusively from the tensor component H T . 34

⊤ L L Proof. Item 1: H GN is a sum of matrices of the form Dc←v Hout Dc←w , where Hout ⪰ 0 for convex L (see proof of Theorem 7). P Item 2: the tensor component contains terms i Tu;v,w δu,i where δu,i can be negative (e.g. under misclassification) and the tensors T need not be PSD. Item 3: the tensor component H T contains cross-blocks Huf1 ,u2 (u1 ̸= u2 ) through tensor terms T combined with backpropagated residuals δu,i , which can be negative. The tensors T need not be PSD.

Remark 84 (Connection to the ReLU paradox). For piecewise-linear activations (ReLU, Leaky ReLU) Tv = 0 T f GN a.e. (the second derivative of a piecewise-linear function is zero), so Hv,w = 0 and Hv,w = Hv,w ⪰ 0 in − f activation space. Consequently m (H ) = 0 a.e.: the input Hessian contains no negative curvature. The parametric Hessian equation Appendix F.1 retains a residual term determined by the structure of node functions, but does not contain σ ′′ (see GN-Gap metric, Definition 11). Proposition 85 (Escape directions from saddle points). Let H f ull have spectral decomposition H f ull = P ⊤ i λi vi vi . Then the eigenvectors {vi : λi < 0} corresponding to negative eigenvalues define escape directions from the saddle point. Moving along ±vi (for λi < 0) is guaranteed to decrease the loss (for a sufficiently small step), since curvature is symmetric w.r.t. the sign of vi . This is a classical property of saddle points; the key contribution of the present work is proving that negative curvature is contained exclusively in the tensor component H T (Theorem 83). Proof. By the Taylor expansion: 2

L(θ + αvi ) = L(θ) + α∇L⊤ vi + α2 vi⊤ Hvi + O(α3 ). At a saddle point ∇L = 0, so

L(θ + αvi ) − L(θ) = α2 λi + O(α3 ). For λi < 0 and sufficiently small α > 0: L(θ + αvi ) < L(θ). Corollary 86 (Practical application). For efficient escape from saddle points, an optimizer can: (1) compute the full Hessian (not just the GN approximation); (2) find the smallest eigenvalue λmin and the corresponding eigenvector vmin ; (3) if λmin < −τ (for a threshold τ > 0), take a step along ±vmin . At a strict saddle (∇L = 0) the sign does not matter. In practice, in saddle neighborhoods (∇L ̸= 0 but ∥∇L∥ is small), choose the sign ensuring descent: −sign(∇L⊤ vmin ) · vmin . This decomposition provides a formal tool for analyzing negative curvature directions (cf. Dauphin et al. 2014). Appendix H.3. Functional-analytic properties of the Hessian Under Assumption 2: in the smooth case (A), ∇2θ L is continuous on RP ; in the non-smooth case (B), the Clarke Hessian coincides with the ordinary one a.e., and the AD-Hessian at singular points is an element of a conservative field (Bolte and Pauwels, 2021). Appendix H.4. Integration of specialized architectural components Theorem 87 (Integration of specialized layers). The following architectural components can be represented as DAG nodes and included in the framework: 1. Batch Normalization: represented as a node with two parameter types (scale and shift) and additional internal variables (batch statistics). 2. Residual connections (ResNet): modeled via parallel paths in the graph with subsequent merging. 3. Recurrent networks: mapped to a DAG by unrolling in time, where each time step is a separate subgraph with shared parameters. Note. Attention mechanisms (softmax + weighted sum) are formally representable as a DAG subgraph; an explicit derivation of the tensors T for Softmax, cross-blocks between Q, K, V, and architectural consequences is given in Appendix M (single-head case; the multi-head extension reduces to block diagonalization). 35

Proof sketch. For each layer type one defines the node functions gv and their first and second derivatives. For example, for Batch Normalization: x − µB + β, gv (x, γ, β) = γ p 2 σB + ϵ 2 where µB , σB are the batch mean and variance, γ, β are scale and shift parameters. The Jacobians Dv and second-derivative tensors Tv are computed by standard differentiation rules for each node type, after which the general formulas equation 2 and equation 5 apply.

Appendix H.5. Stochastic nodes Remark 88 (Reduction to a deterministic DAG). Stochastic nodes fv ∼ p(fv | fPa(v) , θv ) reduce to the deterministic formalism via the reparameterization trick: fv = h(fPa(v) , θv , ε), ε ∼ p(ε). After this, ε becomes an additional input node of the DAG without parameters, and formulas equation 2,equation 5 apply without modification to Eε [L]. Thus, the stochastic case requires no separate theory — it is a special case of the deterministic DAG formalism. Appendix I. Full Hessian computation algorithm Appendix I.1. Recursive HVP via graph recursion The graph recursion naturally yields a recursive algorithm for computing the Hessian–vector product (HVP), analogous to the classical Pearlmutter algorithm (Pearlmutter, 1994) but operating at the level of f input blocks Hv,w . For a direction r ∈ RP , the product ∇2θ L · r is computed in a single pass: X X X f Dw′ rθw′ + Tv;i δv,i rθv , (Appendix I.1) Hv,w [∇2θ L · r]θv = Dv⊤ w

w′

f where the sum over w runs over nodes with common descendants and Dv = ∂fv /∂θv . The product Hv,w · P f ( w Dw rθw ) does not require explicit storage of Hv,w ; a vector of dimension dv suffices. The total cost equals one forward+backward pass: O(Tfwd + Tbwd ).

Appendix I.2. Full algorithm Algorithm 5 Full Hessian computation for a neural network Require: DAG G = (V, E), functions {gv }, parameters {θv }, loss L Ensure: Full Hessian ∇2θ L 1: Forward pass: compute fv for all v ∈ V f 2: δout ← ∇fout L; Hout,out ← ∇2 L(fout ) 3: for v ∈ V \ {out} in reverse topological order do P f f f f ⊤ ⊤ 4: Hv,out ← u∈Ch(v) Du←v Hu,out ; Hout,v ← (Hv,out ) 5: end for f 6: Initialize Hv,w = 0 for v, w ̸= out 7: for v ∈ V in reverse topological order do 8: Compute δv by chain rule 9: for w ̸= out: Ch(v) ̸= ∅ or Ch(w) ̸= ∅ do f 10: Compute Hv,w via equation 2 11: end for 12: end for 13: for v ∈ V do 14: for w: ∃ paths v → u, w → u do 15: Compute Hθv ,θw via equation 5 16: end for 17: end for 18: // Weight sharing 19: for θ ∈ SharedParams do 20: Vθ ← NodesWithParam(θ) P 21: Hθ,θ ← v,w∈Vθ Hθv ,θw 22: Update block of H for θ 23: end for 24: Assemble blocks into ∇2θ L 25: return ∇2θ L

▷ Boundary blocks

▷ Internal blocks

Appendix J. Full proof of the computational complexity theorem P Proposition 89 (Computational complexity). Let |V | = n, P = v∈V pv , d = maxv dv , s = maxv |Pa(v) ∪ Ch(v)|. Then: 1. The time complexity of computing the full Hessian is O(nsd3 + nsd2 P + P 2 ) in the general case with dense tensors. 2. For networks with element-wise activation functions (ReLU, sigmoid), where tensors Tu;v and Tu;v,w are diagonal or sparse with O(d) cost, the total time reduces to O(nsd2 + nsdP + P 2 ). 3. The space complexity of storing the full Hessian is O(P 2 ). 4. For a fully connected DAG (s = O(n)) the time complexity is O(n2 d3 + n2 d2 P + P 2 ). f Proof. 1. Computing the input Hessian Hv,w : By formula equation 2, for each pair (v, w): compute Jacobians Du1 ←v and Du2 ←w for all (u1 , u2 ) ∈ Ch(v) × Ch(w): O(s2 d2 ) per pair; compute blocks Huf1 ,u2 (recursively) and multiply Du⊤1 ←v Huf1 ,u2 Du2 ←w : O(s2 d3 ); compute mixed-derivative tensor contractions for u ∈ Ch(v) ∩ Ch(w): O(sd3 ). Per pair: O(s2 d3 ). The number of pairs with non-empty Ch(v) ∪ Ch(w) is at most O(ns), giving O(ns3 d3 ) worst-case. For s = O(1) (chains, trees) this reduces to O(nd3 ), and the overall cost is O(nsd3 ). 2. Computing the parametric Hessian Hθv ,θw : Total for all pairs: O(nsdP ). If Tu;v are diagonal (element-wise activations), the cost reduces to O(nsdP + P 2 ). 3. Assembly: O(P 2 ) for placing blocks. Total: O(nsd3 + nsd2 P + P 2 ).

Corollary 90 (Summary of computational costs). Let n = |V |, P the total number of parameters, d = f ∈ Rdv ×dw is computed without maxv dv , s = maxv |Pa(v) ∪ Ch(v)|. (1) Metrics R, C, D: the block Hv,w building the full P × P Hessian; when even the block cannot be stored, the Hutchinson stochastic estimaf ∥2F = Ez [z ⊤A⊤Az], where each Az is one HVP, O(P ) time, tor (Avron and Toledo, 2011) is used: ∥Hv,w O(P ) memory. For the stable rank D, one additionally needs ∥H∥22 via power iteration (T iterations, 2T HVPs); total cost O((m + 2T ) · P ) (Algorithm 4). (2) Full Hessian: O(nsd3 + nsd2 P + P 2 ); with diagonal tensors (ReLU): O(nsd2 + nsdP + P 2 ). (3) HVP: O(Tfwd + Tbwd ); for fully connected layers: O(P ). Appendix K. Low-rank approximation and cost reduction f The low-rank structure of inter-layer blocks Hv,w is a mathematical property of DAG architectures following from Theorem 77 on exponential resonance decay and Corollary 56 on the rank bound. f Definition 91 (Low-rank factorization). For a block Hv,w ∈ Rdv ×dw , define the rank-r approximation: f f ⊤ Hv,w ≈ H̃v,w := Uv,w Sv,w Vv,w ,

(Appendix K.1)

where Uv,w ∈ Rdv ×r , Sv,w ∈ Rr×r diagonal, Vv,w ∈ Rdw ×r , and r ≪ min(dv , dw ). The optimal Frobenius-norm approximation is given by truncated SVD. f Proposition 92 (Recursive factor update). If Hu,u = Uu Su Vu⊤ is a rank-r factorization, then the GN f component of Hv,w inherits low-rank structure: GN Hv,w =

⊤ Du←v Uu Su Vu⊤ Du←w = Ũ S̃ Ṽ ⊤ ,

(Appendix K.2)

where Ũ = [Du⊤1 ←v Uu1 , . . .] has at most kr columns (k = |Ch(v)| · |Ch(w)|). A fixed rank r is maintained via repeated SVD truncation.

f Theorem 93 (Approximation error bound). Let σ1 ≥ σ2 ≥ · · · ≥ σdmin be the singular values of Hv,w (dmin = min(dv , dw )). Then for the rank-r approximation: v u dmin u X p f f ∥Hv,w − H̃v,w ∥F = t (Appendix K.3) σi2 ≤ dmin − r · σr+1 . i=r+1

In particular, if σr+1 ≤ εσ1 : f f ∥Hv,w − H̃v,w ∥F f ∥Hv,w ∥F

p ≤ ε dmin − r.

Proof. By the Eckart–Young–Mirsky theorem, truncated SVD gives the best rank-r approximation Pdmin in 2Frobef f σi . The nius norm. The error equals the Frobenius norm of the discarded part: ∥Hv,w − H̃v,w ∥2F = i=r+1 bound via σr+1 follows from σi ≤ σr+1 for i > r. Corollary 94 (Rank selection). For relative accuracy ε it suffices to choose rank r such that n o σk+1 ε r = min k : ≤√ . σ1 dmin − k Due to exponential resonance decay (Theorem 77), for pairs (v, w) with large dist(v, w) a very small rank r = O(1) suffices. Remark 95 (Complexity consequences). The low-rank structure guaranteed by Theorem 77 and Corollary 56 automatically reduces the cost of working with inter-layer blocks: storage O(dv dw ) → O((dv + dw )r); multiplication O(dv dw ) → O((dv + dw )r). Combined with the block-banded structure (Corollary 78), the full Hessian is stored in O(n · kε · d · r). Remark 96 (Cost reduction methods). Approaches for reducing the computational cost of the full Hessian: (1) Block approximation: computing only diagonal blocks Hθv ,θv reduces cost to O(nd3 +P d2 ). (2) Lowrank approximation: approximating off-diagonal blocks by rank-r products reduces cost to O(n2 d3 + n2 d2 r + P r). (3) Gauss–Newton approximation: using only the first term in equation 2 and equation 5 reduces cost and guarantees PSD. (4) Kronecker factorization: representing matrix blocks as Kronecker products of smaller matrices. Remark 97 (Stability and invariance). The relative recursion error under local perturbations of Jacobians (ϵD ) and tensors (ϵT ) grows linearly with depth: εv,w ≤ dist(v, w)(2ϵD + ϵT ) + O(ϵ2 ). At standard machine precision the protocol is stable for networks of depth L < 106 (float32). Geometric coupling C(v, w) is invariant under weight rescaling Wv 7→ αWv , Ww 7→ α−1 Ww . Detailed analysis: Appendix B. Remark 98 (Convergence of optimization methods). Using the full Hessian (or its regularized Clarke analog Ht + λI) in Newton’s method ensures quadratic/superlinear convergence in a neighborhood of regular minima. Details: Appendix C. Appendix L. Continuity of the Clarke Hessian Theorem 99 (Hausdorff continuity for piecewise-linear networks). Let F : Rd → Rm be a piecewise-linear 2 neural network (ReLU activations), ℓ ∈ C 2 (Rm , R), L = ℓ ◦ F . Then the set-valued map x 7→ ∂C L(x) is upper semi-continuous in the Hausdorff metric. 2 Proof sketch. ∂C L(x) = co{lim ∇2 L(xk ) : xk → x, xk ∈ D}. For piecewise-linear F the number of linear 2 regions is finite (Serra et al., 2018); inside each region ∇2 L is continuous. On a region boundary ∂C L(x) is the convex hull of finitely many limits (one per adjacent region), and upper semi-continuity follows by standard arguments (Rockafellar and Wets, 1998, Thm 5.19). Continuity of the AD-Hessian at smooth points given continuity of this set follows from (Rockafellar and Wets, 1998, Thm 1.17).

Appendix M. Example: Attention block (Softmax × Value) To confirm the applicability of the framework to Transformer architectures, we analyze a single-head Attention block as a DAG subgraph. Appendix M.1. DAG representation Single-head Attention: input X ∈ Rn×d , three linear nodes: Q = XWQ ,

K = XWK , d×dk

WQ , WK ∈ R

V = XWV ,

, WV ∈ Rd×dv .

√ P Softmax node: A = σ(QK ⊤ / dk ) ∈ Rn×n , where σ(Z)ij = exp(Zij )/ k exp(Zik ). Output node: O = AV ∈ Rn×dv . Remark 100 (DAG structure of Attention). In terms of our formalism, the Attention block is a DAG with 5 nodes: {vin , vQ , vK , vV , vout }, where vout combines Softmax and multiplication by V . Parameters: θQ = WQ , θK = WK , θV = WV . Parents of vout : Pa(vout ) = {vQ , vK , vV }. Appendix M.2. Softmax Jacobian and Hessian

√ The Softmax Jacobian w.r.t. logits z ∈ Rn (fixing a row of QK ⊤ / dk ): ∂σi = σi (δij − σj ) =: Sij , ∂zj

(Appendix M.1)

where S = diag(σ) − σσ ⊤ is a dense n×n matrix. The Hessian (tensor T in our notation) for the i-th output: [Tσ;z ]i,j,k =

∂ 2 σi = σi δij δik − δij σk − δik σj − δjk σj + 2σj σk . ∂zj ∂zk

Remark 101 (Density of Softmax cross-blocks). Since σi > 0 for all i, both S and Tσ;z are dense, meaning Softmax creates dense Hessian cross-blocks linking all sequence positions. This qualitatively distinguishes Attention from ReLU layers, where T is diagonal. Appendix M.3. Cross-blocks between Q, K, and V Applying formula equation 2 to node vout with Pa(vout ) = {vQ , vK , vV }: Proposition 102 (Attention cross-blocks). Cross-blocks of the input Hessian for the Attention block: f 1. HQ,V : curvature from first derivatives (Softmax Jacobian S, multiplied by V ). Since O = AV , differentiation w.r.t. Q passes through A, w.r.t. V directly. The block is dense (all entries nonzero) along the positional index. f 2. HQ,K : both affect A through QK ⊤ ; the cross-block contains the tensor component Tσ;z (Softmax Hessian), dense along the positional index. This is the main source of H T inside the Attention block. f f 3. HK,V : analogous to HQ,V (symmetric role of K in QK ⊤ ). f Sketch. By formula equation 2 the cross-block HQ,K contains terms of the mixed (item 2 of equation 2) and propagated (item 1) type. For the pair (Q, K): the common descendant vout has the mixed tensor Tout;Q,K that includes the Softmax Hessian [Tσ;z ]i,j,k from equation Appendix M.2. Since σi > 0 for all i, the tensor T and the resulting block are dense.

Appendix M.4. Consequences for architectural analysis Remark 103 (Structural consequences of Attention). 1. Loss of sparsity: Softmax destroys the Hessian sparsity inherent to ReLU networks. The blockbanded approximation (Corollary 78) is inapplicable inside an Attention block: all tokens exchange curvature through the dense tensor S. The effective dimension deff inflates along the sequence-length axis. 2. Q–K curvature: the main source of the tensor component H T in Transformers is the Softmax Hessian between Query and Key. This explains empirical observations of Transformer training instability without warmup or normalization: a large GN-Gap for the (Q,K) pair indicates significance of H T . 3. Attention as a fully connected graph: along the temporal/positional axis, Attention acts as a fully connected DAG subgraph, explaining global training coordination but making cross-blocks computationally dense (n2 nonzero entries). The full multi-head analysis (multi-head attention with parameter sharing across heads and output projection WO ) reduces to block diagonalization by heads followed by a linear projection and requires no new formalism. Appendix N. Additional experimental results Table Appendix N.1 reports the metrics Dfar and Cfar for BottleneckMLP with L = 8 (protocol identical to Exp. 2 in the main text, L = 6). Qualitatively the pattern persists: Dfar grows monotonically with du . init Quantitatively, Dfar at L = 8 is 8–20% lower for du ≥ 64, consistent with additional rank compression from the extra layers between the bottleneck and peripheral blocks. After training the difference is within seed final variance (Dfar differs by ≲ 10 % without systematic direction). Table Appendix N.1: Exp. 2: bottleneck ablation (L = 8, CIFAR-100). Same protocol as Table 3 (mean ±1σ over 5 seeds). du

init Dfar

final Dfar

init Cfar

final Cfar

Acc (%)

4 8 16 32 64 128 256

1.58 ± 0.43 2.43 ± 0.30 4.27 ± 0.57 6.36 ± 0.69 6.40 ± 0.60 7.90 ± 0.28 8.60 ± 0.51

2.44 ± 0.17 3.92 ± 0.39 5.01 ± 0.62 5.13 ± 0.45 5.41 ± 0.55 4.89 ± 0.53 5.73 ± 0.47

0.20 ± 0.08 0.34 ± 0.08 0.43 ± 0.03 0.51 ± 0.03 0.58 ± 0.02 0.65 ± 0.01 0.69 ± 0.01

0.33 ± 0.01 0.45 ± 0.04 0.56 ± 0.01 0.62 ± 0.01 0.68 ± 0.03 0.76 ± 0.02 0.79 ± 0.02

9.6‡ 14.1 17.7 19.1 20.5 21.4 22.2

512†

9.21 ± 0.65

4.89 ± 0.33

0.70 ± 0.02

0.81 ± 0.01

22.4

Appendix N.1. Exp. 1b: Verification of the theorem under spectral normalization Motivation. In the standard regime (He init, width = 64) ρmax > 1 and the sufficient condition of Theorem 26 is not satisfied (Section 5.1). For direct verification of the theorem in its own regime we repeat the Exp. 1 protocol with two modifications: (i) spectral normalization (Miyato et al., 2018) on all Linear layers of Plain MLP, guaranteeing ∥Wi ∥2 = 1; (ii) LayerNorm disabled to ensure piecewise linearity. With ReLU activation: ∥Dblock ∥2 = ∥diag(σ ′ ) · W ∥2 ≤ ∥W ∥2 = 1, so ρ ≤ 1 by construction; for a chain s = maxv |Ch(v)| = 1 and the condition sρ = ρ ≤ 1 holds (strictly at initialization: ρmax ≈ 0.91). Protocol. Plain MLP (width 64, depths L ∈ {8, 12, 16}, CIFAR-10, 50 epochs) with spectral normaliza¯ tion. Residual MLP as control (no SN, as in Exp. 1). All metrics R̄(d), C(d), ρ, λ1 computed by the same protocol. Results (Table Appendix N.2, Figure Appendix N.1). (a) ρmax ≤ 1 at initialization for all depths; after training ρmax ≈ 1.0 (slight excess ≤ 1.04 at L = 8 due to finite precision of spectral normalization). 40

(b) At initialization (ρ ≈ 0.91 < 1), R̄(d) decays exponentially with d for all depths — a direct verification of inequality equation 10. After training the decay rate is determined by λ1 : strong decay for L ≥ 12 (λ1 ≈ −1.0), weak for L = 8 (ρ → 1, boundary of the theorem condition). (c) The Lyapunov exponent λ1 < 0 for all configurations: λ1 ≈ −1.0 at init and −1.0–−0.2 after training — the network operates in the contracting regime. (d) Accuracy of Plain SN drops sharply with depth (36.3% at L = 8; 10.0% at L ≥ 12), confirming that strict curvature suppression (ρ ≤ 1) leads to vanishing curvature and loss of trainability. Residual MLP (control) stably reaches ≈ 51%. Table Appendix N.2: Exp. 1b: spectral normalization. Same metrics as Table 2. Plain SN guarantees ρ ≤ 1 and λ1 < 0; Residual MLP as control without SN (CIFAR-10, mean over 5 seeds). L

Arch.

ρinit max

ρfinal max

λinit 1

λfinal 1

Acc (%)

8 8 12 12 16 16

Plain SN Res. Plain SN Res. Plain SN Res.

0.91 1.25 0.91 1.20 0.91 1.17

1.00 1.33 0.89 1.25 0.90 1.20

−1.01 0.02 −1.01 0.01 −1.01 0.01

−0.24 0.04 −1.03 0.03 −1.03 0.02

36.3 51.3 10.0 51.3 10.0 51.0

L = 8 init L = 12 init L = 16 init L = 8 final L = 12 final

R̄(d)

10−2

10−5

10−8 0

Distance d Figure Appendix N.1: Exp. 1b: resonance R̄(d) under spectral normalization (Plain SN). Solid curves — initialization (ρ ≈ 0.91 < 1): exponential decay over ∼5.5 orders of magnitude at L = 16, directly verifying inequality equation 10. Dashed — after training: at L = 8 (ρ → 1.0) R plateaus at ≈0.94 (boundary of the theorem condition); at L = 12 decay persists (λ1 < 0). Error bands: ±1σ over 5 seeds.

Appendix N.2. Exp. 1: depths L = 16 and L = 32 Table Appendix N.3 extends the Exp. 1 results to depths L ∈ {16, 32}. Exponential decay of R̄(d) persists up to L = 32: the ratio R̄(0)/R̄(L−1) reaches 5.9× (L = 32, plain, init), and the normalized ¯ coupling C(L−1) drops to 0.145. Residual MLP demonstrates stability of R̄(d) for all depths and training ¯ phases (C(L−1) > 0.91). For L = 32 (plain, final), R̄(d) exhibits a U-shaped pattern: decay from 1.75 to a minimum of ≈0.33 at d = 14, then rise to 0.52 at d = 31. The cause is non-uniform self-curvature: layers near the output train more and have larger R(v, v); at large d few pairs remain, biased toward high-resonance output layers. The ¯ normalized coupling C(d) corrects this artifact and decays monotonically (R2 > 0.91).

Table Appendix N.3: Exp. 1: metrics for depths L ∈ {16, 32} (CIFAR-10, width = 128, mean ±1σ over 5 seeds; stochastic ¯ estimation: 100 Rademacher probes, subsample 64). R̄(0)/R̄(L−1) — resonance of nearest/most distant pairs; C(L−1) — coupling at maximum distance. ¯ L Arch. Phase ρmax R̄(0) R̄(L−1) C(L−1) 16 16 16 16 32 32 32 32

Plain Plain Res. Res. Plain Plain Res. Res.

init final init final init final init final

2.59 ± 0.19 2.56 ± 0.03 1.17 ± 0.00 1.21 ± 0.01 2.71 ± 0.13 2.62 ± 0.03 1.12 ± 0.00 1.13 ± 0.00

0.072 ± 0.018 4.852 ± 1.219 0.121 ± 0.005 2.156 ± 0.111 0.067 ± 0.018 1.750 ± 0.530 0.120 ± 0.004 2.176 ± 0.082

0.021 ± 0.001 0.80 ± 0.13 0.12 ± 0.00 1.44 ± 0.08 0.011 ± 0.003 0.52 ± 0.16 0.12 ± 0.00 1.48 ± 0.05

0.263 0.624 0.983 0.887 0.145 0.403 0.986 0.915

Appendix N.3. Exp. 4: GN-Gap by branch depth k ∈ {1, 2, 3} Table Appendix N.4 extends the Exp. 4 results from the main text (k = 2) with data for k = 1 and k = 3. The qualitative picture is independent of k: cat+SiLU (nonlinear merging, σ ′′ ̸= 0) exhibits GapGN > 1 at initialization for all k; remaining configurations stay at machine precision (∼ 10−8 ). The variation of Gap across k for cat+SiLU is less than 3% at init and ∼15% at final, confirming the robustness of the conclusion about the dominant role of merging type and activation, not branch depth. Table Appendix N.4: Exp. 4: GN-Gap at the merging node of Diamond MLP for k ∈ {1, 2, 3} (width = 32, CIFAR-10, mean over 5 seeds). k

Configuration

Gapinit

Gapfinal

∥T ∥init F

∥T ∥final F

1 1 1 1

sum+ReLU sum+SiLU cat+ReLU cat+SiLU

6.5·10−8 6.5·10−8 8.1·10−8 1.35

8.7·10−8 9.2·10−8 1.2·10−7 0.078

3.6·10−9 3.5·10−9 6.0·10−10 6.6·10−3

2.5·10−8 2.6·10−8 5.6·10−8 3.3·10−2

2 2 2 2

sum+ReLU sum+SiLU cat+ReLU cat+SiLU

6.7·10−8 7.0·10−8 8.4·10−8 1.33

9.3·10−8 9.6·10−8 1.2·10−7 0.082

4.1·10−9 4.1·10−9 6.6·10−10 6.8·10−3

2.0·10−8 2.0·10−8 3.3·10−8 2.5·10−2

3 3 3 3

sum+ReLU sum+SiLU cat+ReLU cat+SiLU

6.5·10−8 7.1·10−8 8.3·10−8 1.32

8.9·10−8 9.4·10−8 1.4·10−7 0.090

3.8·10−9 4.1·10−9 7.6·10−10 7.0·10−3

1.7·10−8 1.9·10−8 3.0·10−8 2.1·10−2

Appendix N.4. Exp. 3: decomposition of ∥H T ∥F and ∥H GN ∥F by distance Table Appendix N.5 reports the Frobenius norms of the tensor and Gauss–Newton components of the inter-layer Hessian by distance d for three smooth activations (init, L = 6, width = 64). Both components decay exponentially with d; at maximum distance d = L−1 = 5, ∥H T ∥F → 0 (structural zero: at the top level of recursion equation 2 tensor terms are absent; see Proposition 14).

Table Appendix N.5: Exp. 3: ∥H T ∥F and ∥H GN ∥F by distance d (init, L = 6, width = 64, CIFAR-10, mean over 5 seeds). At T d = 5 the tensor component vanishes (Hd=L−1 ≡ 0). d

Activation

∥H T ∥F

∥H GN ∥F

GapGN

0 0 0

Softplus SiLU GELU

1.50·10−2 2.65·10−2 4.22·10−2

1.22·10−1 1.17·10−1 1.17·10−1

0.123 0.227 0.361

2 2 2

Softplus SiLU GELU

1.24·10−3 2.15·10−3 3.45·10−3

1.05·10−2 9.59·10−3 9.70·10−3

0.119 0.224 0.357

4 4 4

Softplus SiLU GELU

7.68·10−5 1.32·10−4 2.15·10−4

8.63·10−4 7.88·10−4 8.04·10−4

0.089 0.168 0.269

5 5 5

Softplus SiLU GELU

< 10−10 < 10−10 < 10−10

2.37·10−4 2.08·10−4 2.15·10−4

< 10−7 < 10−7 < 10−7

Appendix N.5. Exp. 4: cross-branch resonance R̄AB sum+ReLU init

cat+ReLU init

sum+SiLU init

cat+SiLU init

sum+ReLU final

cat+ReLU final

sum+SiLU final

cat+SiLU final

R̄AB (dgraph )

100

10−1

10−2 2

Graph distance dgraph Figure Appendix N.2: Exp. 4: cross-branch resonance R̄AB (dgraph ) in Diamond MLP (k = 2). Solid curves — initialization (decay with distance, H4.2a); dashed — after training (growth, H4.2b). All 4 configurations show R growth with dgraph at final (3.5–5.7× from d = 2 to d = 4). Error bands: ±1σ over 5 seeds.

Appendix N.6. Ablation: exact vs. stochastic estimation To assess the sensitivity of the metrics to the Hessian computation method, we compare exact decomposition (Exp. 3) with Hutchinson stochastic estimation using 30 Rademacher probes on the same configurations. Table Appendix N.6 reports the relative discrepancy of GN-Gap for smooth activations. The discrepancy does not exceed 0.82% at init and 0.33% at final — the stochastic estimator systematically overestimates Gap by less than 1%, confirming the adequacy of the Hutchinson estimator with 30 probes for diagnosing the tensor component. In Exp. 1 at L ≥ 12 the stable rank D < 1 for distant pairs, indicating degeneration: the ratio ∥H∥2F /∥H∥22 falls below unity due to estimation noise in ∥H∥F . For deep or wide networks we recommend scaling the number of probes as m ∝ min(dv , dw ) or using exact block computations when dimensionality permits. The control run of Exp. 2 with dbase = du = 512 (uniform network without bottleneck; Table 3, † ) shows init final Dfar = 11.2 and Dfar = 5.2 — in the absence of a narrow layer the stable rank is determined by network width and number of classes, consistent with the theoretical limit D ≤ min(du , K−1). 43

Table Appendix N.6: Ablation: GN-Gap discrepancy between exact and stochastic (30 probes) estimation (Exp. 3, L = 6, width = 64, mean over 5 seeds). Activation

Gapinit exact

Gapinit stoch

∆init (%)

∆final (%)

Softplus SiLU GELU

0.122 0.225 0.358

0.123 0.226 0.359

0.82 0.46 0.43

0.25 0.33 0.31

Appendix N.7. Exp. 5: Toy-Attention protocol and additional details S×d Architecture. (S = 8, d = 16), single-head self-attention O = √ ToyAttentionModel: input X ∈ R ⊤ Softmax(QK / d) V , mean-pool, Linear(d, 1). Parameters: 3d2 + d + 1 = 785. ToyReluMLP: input X ∈ RS×d , 3×[Linear(d, d)+ReLU] applied identically at each position (shared weights), mean-pool, Linear(d, 1). Parameters: 3(d2 +d) + d + 1 = 833. Both models process input per-position and aggregate via mean-pool; the only structural difference is the nonlinearity (Softmax attention vs. position-wise ReLU). Synthetic data. Teacher network with fixed random weights (seed = 0, independent of training seeds): S

1X tanh x⊤ s Wt wr + ε, S s=1

Wt ∈ Rd×d , wr ∈ Rd , ε ∼ N (0, 0.01).

Wt,ij ∼ N (0, 1/d), wr,j ∼ N (0, 1/d). Training set: 2048 samples; validation: 512; targets standardized per training split. GN-Gap computation. Both models expose a forward_with_intermediates method returning the graph of intermediate activations. For Attention, the pair (v, w) = (Q, K); for ReLU-MLP, (block0 , block1 ). The cross-Hessian ∂ 2 L/∂fv ∂fw is computed exactly via second-order autograd; the Gauss–Newton compoPB (b)⊤ (b) GN = (2/B) b=1 Jv Jw (MSE loss with mean reduction, scalar output). nent is Hv,w Appendix N.8. Exp. 6: ResNet-18 protocol and additional results Architecture. SegmentedResNet18: standard torchvision ResNet-18 (∼11M parameters) with 5 layer names (stem, layer1–layer4) and 6 segments. seg5 = avgpool + flatten + fc is purely linear (no activation). SegmentedPlainResNet18: identical parameterization with identity shortcuts replaced by sequential conv– BN–σ paths; stride-based downsampling preserved. Training details. SGD (η = 0.1, momentum 0.9, weight decay 5·10−4 ), cosine schedule, 100 epochs, batch 128. CIFAR-10 augmentation: random crop 32×32 with padding 4, random horizontal flip, normalization to channel mean/std. Five seeds {42, . . . , 46}. Stochastic estimation. Hutchinson estimator: m = 30 Rademacher probes, subsample of 32 examples. ˆ 2 < 10−24 ⇒ D = NaN (segment contributes Power iteration for ∥H∥2 : T = 50 steps. NaN threshold: ∥H∥ 2

negligible curvature). GN-Gap estimator: m = 30 probes, common probe vector, analytical J ⊤ L′′ J computation. Mid-checkpoint results. Table Appendix N.7 provides the full set of metrics at the mid checkpoint (epoch 50), complementing the init and final data in the main text.

Table Appendix N.7: Exp. 6 (mid, epoch 50): metrics R, C, D by distance d (mean over 5 seeds). Test accuracy: ReLU ResNet 80.8 %, ReLU Plain 77.0 %, SiLU ResNet 79.9 %, SiLU Plain 78.4 %. Condition

R(d)

C(d)

D(d)

ReLU ResNet

0 1 2 3 4

1.05 0.87 0.78 0.72 0.39

1.00 0.93 0.85 0.79 0.73

2.87 2.73 2.65 2.85 2.91

ReLU Plain

0 1 2 3 4

3.81 2.16 1.80 1.73 0.41

1.00 0.90 0.81 0.71 0.66

2.52 2.40 2.29 2.47 3.06

SiLU ResNet

0 1 2 3 4

2.60 1.68 1.27 1.64 0.49

1.00 0.90 0.83 0.74 0.61

3.47 3.33 3.30 3.30 2.95

SiLU Plain

0 1 2 3 4

3.89 2.30 1.61 1.87 0.45

1.00 0.87 0.79 0.72 0.60

2.79 2.64 2.29 2.46 2.51

GN-Gap across checkpoints. Table Appendix N.8 reports GN-Gap(d) for all three checkpoints. For ReLU conditions, Gap < 10−5 at all checkpoints, consistent with Proposition 10. For SiLU conditions, the training dynamics of Gap are architecture-dependent: in Plain networks, Gap grows from near-zero at init to 0.1–0.2 at small distances (d ≤ 1) as nonlinear curvature emerges, but slightly declines at larger d; in ResNet, skip connections yield a substantial init Gap (0.4–0.6 at d ≥ 1) that decreases during training as H GN curvature dominates. At d = 4 (linear head), Gap ∼10−6 for all conditions—the analogy with ReLU is explained by linearity of avgpool + fc (Remark 39). Table Appendix N.8: Exp. 6: GN-Gap(d) across checkpoints (mean over 5 seeds). (H T ≡ 0).

† At d = 4 the segment is purely linear

Condition

d=0

d=1

d=2

d=3

d = 4†

ReLU ResNet ReLU Plain SiLU ResNet SiLU Plain

<10−5

<10−5 0.426 0.003

<10−5 0.526 0.028

<10−5 0.619 0.080

<10−5 0.642 0.132

<10−5 <10−5 <10−5 <10−5

mid

ReLU ResNet ReLU Plain SiLU ResNet SiLU Plain

<10−5 <10−5 0.274 0.122

<10−5 <10−5 0.144 0.062

<10−5 <10−5 0.137 0.041

<10−5 <10−5 0.136 0.056

<10−5 <10−5 <10−5 <10−5

final

ReLU ResNet ReLU Plain SiLU ResNet SiLU Plain

<10−5 <10−5 0.426 0.213

<10−5 <10−5 0.241 0.118

<10−5 <10−5 0.220 0.058

<10−5 <10−5 0.150 0.054

<10−5 <10−5 <10−5 <10−5

Checkpoint init

Record · ID 10345 · SHA-256 22bfcdc2eb30e1b8
Conceptio Open Knowledge Archive — every document is proof-bundled with source, license, and retrieval metadata.