QUANTITATIVE GAUSSIAN-PROCESS LIMITS OF TENSOR PROGRAMS ANDREA AGAZZI, ELOY MOSIG GARCÍA, AND DARIO TREVISAN Abstract. We study the infinite-width Gaussian-process limit of random neural networks through the lens of tensor programs, and we provide a quantitative convergence theory in Wasserstein distance. Our main result gives explicit finite-width error bounds, of order inverse square-root of the widths between finite-network executions and their Gaussian-process limits. The framework is architecture-agnostic and covers feed-forward models together with weight-sharing schemes relevant for recurrent and transformer-type architectures.
arXiv:2607.06290v1 [cs.LG] 7 Jul 2026
1. Introduction A depth-M Multilayer Perceptron (MLP) with input x ∈ Rd and layer widths (n0 , . . . , nM ) is defined as the map h(M ) : Rn0 → RnM , where h(0) (x) := x and h(ℓ+1) (x) := ϕℓ W (ℓ) h(ℓ) (x) , (1.1) for ℓ = 0, . . . , M − 1, with activation functions ϕℓ : R → R acting componentwise on their input, and parameters (or weights) W (ℓ) ∈ Rnℓ+1 ×nℓ for ℓ ∈ {0, . . . M − 1}. These parameters are typically initialized randomly, drawn independently from a common (layer-dependent) distribution. A standard practice is to choose such distribution as Gaussian centered with O(n−1 ℓ ) variance, for instance (1.2)
(ℓ)
Wij ∼iid N (0, n−1 ℓ ).
Under this choice of scaling, when the widths n1 , . . . , nM −1 go to infinity, the pre-activations g (ℓ) := W (ℓ) h(ℓ) converge to a centered Gaussian process whose kernel is given by an explicit recursion in ϕ and in the input Gram matrix. This is the neural-network Gaussian-process (NNGP) limit, first identified by [Nea96] for M = 2 and later studied for generic integer M ≥ 2 by [Lee+20]. Quantitative versions of these Central Limit Theorem (CLT) results were obtained in [BT22; Fav+25; Tre23] and more recently in [Cel26; GKR26]. This work extends such quantitative convergence results to a significantly richer class of functions called tensor programs defined below. 1.1. Tensor programs. As discussed in [Yan21], the above function class, as well as many other neural network architectures, can be interpreted as a specific instance of abstract algorithms, called Netsor programs. These are a subclass of the more general tensor programs [Yan20b; YH20], consisting of a sequence of lines where variables are declared and then combined via elementary operations. In particular, a Netsor program T is defined by specifying its inputs and its operations as follows: Inputs. The inputs1 of a Netsor program are associated to two different variable types: • (Input) H-vars are a collection of fixed real vectors xk ∈ Rdk
X := {xk }k
generalizing the input data on which the network acts. The dimensions of this collection are stored in set nin := {dk }k . These variables belong to the class of Date: July 8, 2026. 1a set of variables that are declared at the beginning of the program, and whose definition does not rely on other variables. 1
2
ANDREA AGAZZI, ELOY MOSIG GARCÍA, AND DARIO TREVISAN
H-vars, i.e., real, vector-valued quantities interpretable as hidden states of the network. • A-vars denote an abstract collection of matrix-valued parameters W := {W(j) }j , representing the weights (and possibly biases) of the network. Conforming to the notation in [Yan21], we recall the type and dimension of a newly declared variable through the : notation2, so that any A-var W of dimension n × m will be declared as W : A(n, m), while H-vars (and G-vars defined below) of dimension n will be respectively introduced as h : H(n) and g : G(n). Operations. The variables defined above can be combined through the following operations • MatMul: The Matrix-vector multiplication operation combines a W : A(n, m) and a h : H(m). The result of this operation is a new n-dimensional vector-valued H-var behaving as a Gaussian when the dimension being contracted diverges. Denoting the subset of H-vars satisfying this property as G-vars (abstracting the concept of pre-activations), we write an instance of this operation as MatMul
g := Wh : G(n).
• NonLin: Componentwise composition of h1 , . . . , hm : H(n) with a (possibly nonlinear) function ϕ : Rm → R, yielding a new H-var of dimension n. This operation is declared as NonLin h := ϕ(h1 , . . . , hm ) : H(n) The execution of each of the above operations yields a line of the Netsor program. We denote by L ∈ N the total number of lines in a program. Hidden widths and structural constants. The hidden widths of the program are defined as the dimensions over which the program performs internal matrix contractions. More precisely, a dimension n is called a hidden width if there exists a G-var of shape n which is produced by a MatMul line and whose downstream descendants in the graph are later used as the parent of another MatMul line. Equivalently, n is internal to a chain of matrix multiplications. We denote by nhidd the collection of hidden widths, which will be considered as asymptotically large in this work. By contrast, recall that input dimensions nin are kept fixed. Dimensions that are neither input nor hidden may be referred to as output dimensions nout . Furthermore, quantities that remain fixed throughout the asymptotic limit, such as the input space X , the input dimensions nin , the activations ϕ appearing in NonLin lines, and the order of the Wasserstein distance, are collectively referred to as structural constants. We then use c to denote a generic bounding constant that depends exclusively on these structural parameters and is independent or uniform in the hidden widths nhidd . To keep notation simple, we conventionally allow the value of c to change from line to line. Example 1.1 (Shallow neural network). The following 3-line Netsor program yields (1.1) with M = 2. In this case, nin = {d}, nhidd = {n}, nout = {dout }. Program graphs. To better visualize the above concepts, it is useful to associate to a Netsor program a directed graph, encoding the flow of variables through the program lines. The nodes of the graph are the H-vars (hence also the G-vars) labelled by their type and shape. A MatMul line g = Wh contributes a solid directed edge h −→ g labelled by the corresponding A-var W. Thus, repeated use of the same A-var is represented by repeated occurrences of the same edge label. A NonLin line h = ϕ(h1 , . . . , hk ) contributes directed edges hi −→ h which we draw as dashed lines with ϕ as label. 2This notation is borrowed from the variable declaration syntax in computer science and should be
thought of as an alternative to the symbol ∈.
QUANTITATIVE GAUSSIAN-PROCESS LIMITS OF TENSOR PROGRAMS
3
Netsor program 1 Shallow neural network evaluated on a single input Input: x : H(d) Input: W(0) : A(n, d) Input: W(1) : A(dout , n) 1 g(1) = W(0) x : G(n) 2 h(1) = ϕ(g(1) ) : H(n) 3 g(2) = W(1) h(1) : G(dout )
x : H(d)
W(0)
g(1) : G(n)
▷ MatMul ▷ NonLin ▷ MatMul ϕ
input
h(1) : H(n) hidden width n
W(1)
g(2) : G(dout ) output
Fig. 1. Program graph of a shallow fully connected network evaluated on a single input. The input dimension is d ∈ nin , the hidden width is n ∈ nhidd , and the output dimension is dout ∈ nout . Program execution. For a given choice of inputs X , the execution of a program consists in the realization of its H-vars h(j) (including its G-vars g(j) ) as vector-valued random variables, with a given distribution on a common probability space. At the level of generality considered so far, the execution of a Netsor program and the statistical properties of its H-vars can be interpreted in different ways. In this work, we focus on the following two execution rules. In both executions, NonLin is intended as the component-wise application of ϕ, hence they differ only on the concrete implementations of the Netsor operation. • Finite-width: this execution corresponds to the classical evaluation of the feedforward pass of a neural network, with Gaussian initialization of its weights. For every input A-var W : A(n, m), we realize Wij ∼ N (0, 1/m), all independently of each other. Then, given a H-var h ∈ Rm and an A-var W ∈ Rn×m , the MatMul line g = Wh is interpreted as a standard matrix-vector multiplication, defining a random vector g = W h with values in Rn . • Infinite-width: this execution replaces the finite-width matrix multiplications by their Gaussian central-limit approximation, and interprets G-vars ḡ (1) , . . . , ḡ (G) appearing in the program as jointly Gaussian random variables. Precisely, the output of a MatMul line g := Wh with W : A(n, m) and h : H(m) is defined as Gaussian random variable ḡ taking values in Rn , such that, for any previously introduced G-var g′ with corresponding Gaussian ḡ ′ (including the case g′ = g), the cross-covariance Σḡḡ′ is defined as ( 0 if g′ = W′ h′ for an A-var W′ ̸= W (1.3) Σḡḡ′ := h̄∗ h̄′ Idn if g′ = Wh′ , with same A-var W and h′ : H(m)
(1.4)
where for two executions h̄, h̄′ ∈ Rm of H-vars we introduce the infinite-width product notation: 1 h ⊤ ′i E h̄ h̄ , h̄∗ h̄′ := m where with ⊤ we denote the matrix transpose.
1.2. Main result. The two executions presented above are connected by our main result, presenting quantitative distributional convergence estimates in the wide limit min n → ∞ ,
n∈nhidd
4
ANDREA AGAZZI, ELOY MOSIG GARCÍA, AND DARIO TREVISAN
for the finite-width execution h(ℓ) of a generic Netsor program to its infinite-width counterpart h̄(ℓ) . The convergence is established in Wasserstein distance, defined for any p ≥ 1 for a pair of probability distributions µ, ν with finite p-th moments as Z p Wp (µ, ν) := inf ∥x − y∥p dγ(x, y) γ∈Γ(µ,ν)
where ∥ · ∥ denotes the Euclidean norm and Γ(µ, ν) is the set of all couplings of µ and ν. By slight abuse of notation, for a pair of random vectors X, Y with distributions PX , PY we write Wp (X, Y ) := Wp (PX , PY ). Theorem 1.2 (Quantitative CLT). Consider a Netsor program with L lines, denote the H-vars defined by the program with h(ℓ) : H(nℓ ), and write h(ℓ) ℓ=1,...,L , h̄(ℓ) ℓ=1,...,L , for their finite-width and infinite-width executions. Assume that all functions appearing in NonLin lines are Lipschitz continuous. Then, for every p ≥ 1, there exists a constant c < ∞, depending only on p and on the structural constants of the program, such that ! ! (ℓ) (ℓ) X h̄ 1 h ≤c √ . , √ (1.5) Wp √ nℓ nℓ m m∈n ℓ=1,...,L
ℓ=1,...,L
hidd
Our second main result can be regarded as a law of large numbers (LLN) for empirical kernels built from H-vars, and follows from Theorem 1.2, but it also provides a key ingredient in the induction argument for its proof. Corollary 1.3 (Kernel LLN). Assume the hypotheses of Theorem 1.2. Then for any pair of H-vars h, h′ : H(n), X 1 1 ⊤ ′ √ . ≤c h h − h̄∗ h̄′ n m Lp m∈n ∪n hidd
out
where c < ∞ depends only on p and on the structural constants of the program. We remark that the set nout appears in the right-hand side only when one takes empirical kernels over output dimensions, i.e., n ∈ nout . Scope and comparison with existing limits. Let us emphasize that Theorem 1.2 is not merely a reformulation of existing quantitative central limit theorems for fully connected feed-forward networks. The tensor program formulation allows one to treat architectures in which the same weight matrix is used several times in the computation, possibly across different layers or different input branches. This includes, for example, weight-sharing patterns arising in recurrent networks and in the joint evaluation of a network on several inputs (see the examples proposed in Subsection 2.1). In such cases, the relevant finite-width dependencies are not captured by a layer-by-layer independent CLT, and the conditional Gaussian structure of tensor programs becomes essential. The result should also be compared with Yang’s master theorem [Yan21, Theorem 5.4] for tensor programs. That theorem gives a qualitative law-of-large-numbers description of the infinite-width limit for a very broad class of programs and nonlinearities. Our theorem is more restrictive in its regularity assumptions, since we require Lipschitz continuity, but it gives an explicit quantitative Wasserstein estimate between the finite-width execution and the infinite-width Gaussian-process execution. In this sense, the present result can be viewed as a quantitative counterpart, for a restricted, but architecture-rich, class of tensor programs, to the qualitative tensor-program limit theory. Finally, the basic Netsor language does not cover all operations needed for modern architectures, such as transformers. Subsection 2.2 introduces a restricted extension, denoted NetsorK , in which scalar kernel variables and scalar-parametric nonlinearities are
QUANTITATIVE GAUSSIAN-PROCESS LIMITS OF TENSOR PROGRAMS
5
allowed. This extension is designed to cover attention-type layers while keeping the quantitative proof based on kernel LLN estimates and finite-dimensional stability bounds. Structure of the proof. The proof of Theorem 1.2 is a line-by-line induction on the program. One mild subtlety is that the kernel LLN in Corollary 1.3 is both a consequence of the quantitative convergence theorem and an ingredient in its proof. This does not lead to a circular argument: at the induction step from lines 1, . . . , r − 1 to line r, we use the theorem only for the truncated program up to line r − 1. This gives the quantitative coupling of all previously constructed H-vars. Combining this coupling with a standard quantitative law of large numbers for the independent coordinates of the infinite-width execution yields the kernel LLN for all kernels that can appear before line r. Once this kernel control is available, the induction step is as follows. If line r is a NonLin instruction, the estimate is propagated directly by the Lipschitz stability of the activation. If line r is a MatMul instruction, both the finite-width and infinite-width executions are conditionally Gaussian given the previous H-vars. We couple these two conditional Gaussian laws by using the same standard Gaussian noise. The resulting error is reduced to estimating the difference between the corresponding conditional means and square-roots of conditional covariances. The structure lemmas 3.5 and 3.6 show that these conditional means and covariances are finite-dimensional functions of the Gram matrices of the parent H-vars associated with the same A-var. After the non-degeneracy reduction, these functions are differentiable at the infinite-width Gram point. The kernel LLN for the truncated program therefore controls the difference between the finite-width and infinite-width conditional parameters. This closes the induction and proves Theorem 1.2. Applying the same argument to the full program then gives Corollary 1.3. Related works. The convergence of neural networks to Gaussian Processes in the wide limit was first established for shallow, fully connected neural networks by [Nea96] and subsequently extended to deep fully connected networks [GM+18]. These qualitative results, coupled with the introduction of the Neural Tangent Kernel (NTK) [JGH20; Aro+19], provided a rigorous foundation for analyzing the training and generalization of infinite-width networks. To systematically encompass the vast landscape of modern architectures, the Tensor Programs formalism was introduced [Yan20a; Yan21; Yan20b; YH20]. This framework provides a universal algebraic language to prove that the forward and backward passes of essentially any standard architecture converge to Gaussian Process and NTK limits. Importantly, the Tensor Programs formalism also clarifies the distinction between the NTK regime and the feature learning regime [YH21], where representations actively evolve during training. In this work, we focus on the former setting. Parallel to the infinite-width asymptotic, the infinite-depth and joint wide-and-deep limits have been extensively explored to understand hierarchical feature propagation. For instance, [PF20] demonstrated that residual networks converge to diffusion processes in the infinite-depth limit, while [Hay24] formalized the commutativity of sequential width and depth limits. More recently, [Aga+26; KGR26] established quantitative scaling limits for residual attention architectures in the wide, deep, and particle mean-field limits. While foundational results established qualitative asymptotic limits, determining the precise rate of convergence at finite width has become a central focus. From a rigorous probabilistic standpoint, quantitative CLTs characterizing the convergence of finite-width network architectures to their Gaussian wide limits were established for fully connected networks in [BT22; Tre23; BFF25; Fav+25] by using optimal transport theory in the first two cases, second-order Poincaré inequalities in the third, and Stein-Malliavin method in the last. Notably, also functional quantitative CLTs were proven in [EMS21; Cam+23] and extended to deep networks and strengthened in [Fav+25]. Moving beyond initialization,
6
ANDREA AGAZZI, ELOY MOSIG GARCÍA, AND DARIO TREVISAN
[MAT26] extended the optimal transport quantitative bounds from [BT22] to positive training time in the shallow MLP case. Very recently, these quantitative guarantees have been significantly broadened. [Cel26] provided convergence rates with an explicit dependence on general weight distributions and hyperparameters, while [GKR26] bypassed traditional Gaussian assumptions to prove universality in deep networks via the Lindeberg exchange principle. Concurrently, efforts have been made to quantitatively analyze more intricate architectures. Prior to these quantitative efforts, [Hro+20] rigorously derived the exact NNGP and NTK kernels for deep attention networks in both d−1 and d−1/2 scalings. Building upon this, [SKI26] recently leveraged the Tensor Programs framework to analyze the infinite-width limit of attention layers. Notably, they extended the Master Theorem in [Yan21] to attention architectures with the d−1/2 scaling in the dot product, as opposed to the d−1 scaling used by [Yan21]. Finally, we mention a different but related direction, concerning Gaussian-process limits for quantum neural networks. In that setting, the random function generated by a parametrized quantum circuit is typically an expectation value of observables over many qubits, and Gaussian limits arise in a large-width or many-qubit regime under suitable locality or weak-correlation assumptions. Recent work [GDP25; MH+25] has shown that both untrained and trained quantum neural networks may converge to Gaussian processes, and that this convergence can be quantified in Wasserstein distance. Although the mechanisms and architectures differ substantially from the tensor-program setting considered here, these results share the same general objective: to turn qualitative Gaussian-process limits for random neural models into explicit finite-size approximation bounds. Structure of this work. Section 2 illustrates the scope of the main theorem through several examples, including shallow and deep fully connected networks, recurrent architectures, and its extension to NetsorK programs, which adds scalar kernel variables and scalarparametric nonlinearities. Section 3 develops the structural ingredients of the argument: conditional-law formulas for finite- and infinite-width executions, and the reduction to non-degenerate representative programs. Section 4 proves the main results for both base Netsor programs and extended NetsorK programs. Finally, Section 5 presents numerical experiments illustrating finite-width convergence to the Gaussian-process limit across several architectures. 2. Examples and Extensions Before we move to the proof of our main results, we illustrate how they apply to the forward pass of a large family of architectures. 2.1. Examples. We start with simple examples, moving towards more complex ones. Shallow fully connected network. The previously defined Netsor program 1 showcases a single-hidden layer neural network. The corresponding graph is shown in Figure 2. In this case, Theorem 1.2 yields: ! !! h̄(1) ḡ (2) c h(1) g (2) √ ,√ , √ ,√ ≤√ , (2.1) Wp n n n dout dout √ where h̄(1) / n is the scaled componentwise image via ϕ of the random vector ḡ (1)√with entries i.i.d. centered Gaussian variables with variance equal to x⊤ x and G = ḡ (2) / dout is the Neural Network Gaussian process (NNGP) for a shallow network [JGH20; Lee+20; BT22]. ḡ (1) and ḡ (2) are to be interpreted as centered Gaussian variables with covariance as in Equation (1.3).
QUANTITATIVE GAUSSIAN-PROCESS LIMITS OF TENSOR PROGRAMS
x : H(d)
W(0)
g(1) : G(n)
ϕ
input
h(1) : H(n) hidden width n
W(1)
7
g(2) : G(dout ) output
Fig. 2. Program graph of a shallow fully connected network. The dimension n is a hidden width because it is produced by the first MatMul line and later contracted by the readout MatMul line. The input dimension d and output dimension dout are structural. Deep fully connected network. Figure 3 and the corresponding Netsor program 2 extend the shallow architecture to any M -layer deep feed-forward network. The H-vars propagate through a sequence of matrix multiplications parameterized by independent A-vars W(ℓ) and componentwise nonlinearities ϕℓ . In this case Theorem 1.2 provides joint convergence of the G = (G (ℓ) )ℓ=1,...,M ∼ pre-activations to a centered Gaussian process: (M ) n×n N 0, K Iddout , where Idn denotes the identity matrix in R , recursively defining the deep NNGP kernel: c(0) ⊤ c(ℓ) h (ℓ) (ℓ) ⊤ (ℓ) (ℓ) i K (1) = x x, K (ℓ+1) = E ϕ (G ) ϕ (G ) ∈ R, d nℓ (ℓ) where Wij ∼ N 0, c(ℓ) /nℓ , for each i, j ≤ nℓ and ℓ ≤ M − 1. Note that the NNGP kernel may also be written as a function of the H-vars of the program via the infinitewidth product ∗ defined in Equation (1.4): K (ℓ+1) = c(ℓ) h̄(ℓ)∗ h̄(ℓ) , for each ℓ ≤ M. √ Note that, in our notation, G (ℓ) = ḡ (ℓ) / nℓ for each ℓ, so the definition of the NNGP at each layer is consistent with the shallow network example. Restricting the bound in Theorem 1.2 to the G-vars of the program recovers the bound from the main result in [BT22]: ! M −1 (ℓ) X 1 g √ . , G ≤ c (2.2) Wp √ nℓ ni ℓ=1,...,M
i=1
Note that for the fully connected architecture this rate has been improved in [Tre23] to −1 O (mini ni ) . Netsor program 2 Deep feed-forward network Input: x : H(d) Input: W(0) : A(n1 , d) Input: W(ℓ) : A(nℓ+1 , nℓ ) for ℓ = 1, . . . , M − 2 Input: W(M −1) : A(dout , nM −1 ) 1 h(0) = x : H(d) for ℓ = 1, . . . , M do 2 g(ℓ) = W(ℓ−1) h(ℓ−1) : G(nℓ ) 3 h(ℓ) = ϕ(g(ℓ) ) : H(nℓ ) end for
▷ MatMul ▷ NonLin
Time-unrolled RNN with residue. The graph in Figure 4 and the corresponding Netsor program 3 illustrate a recurrent neural network unrolled over two time steps. The structurally defining feature of this architecture is the reuse of the A-vars Wx and Wh across the temporal sequence, which introduces dependencies between the generated H-vars, as
8
ANDREA AGAZZI, ELOY MOSIG GARCÍA, AND DARIO TREVISAN input
x : H(d)
hidden n1
W
(0)
g(1) : G(n1 )
ϕ1
h(1) : H(n1 )
W(1) intermediate layers
hidden n2
g(2) : G(n2 )
ϕ2
h(2) : H(n2 )
W
(2)
...
W(M −2)
g(M −1) : G(nM −1 )
ϕM −1
h(M −1) : H(nM −1 )
W(M −1)
g(M ) : G(dout )
ϕM
y : H(dout ) output
Fig. 3. Program graph of a deep feed-forward network with activated output. The internal widths are n1 , . . . , nM −1 . The terminal output dimension dout is structural, even if the final program line is obtained by a componentwise nonlinearity. opposed to the two previous examples. In particular, the results for fully connected architectures are no longer applicable in this case. Under the Netsor formalism, Theorem 1.2 guarantees that the unrolled hidden states jointly converge to the recurrent NNGP [Yan21] with the expected rate. Specifically, we have y c h1 h2 √ ,√ ,√ , (Z1 , Z2 , G) ≤ √ , (2.3) Wp n n dout n √ where G = ȳ/ dout denotes the recurrent NNGP associated with this architecture and Z1 = ϕ(U1 ), Z2 = ψ(V, U2 ) with h i 1 1 ⊤ 1 ⊤ U1 ∼ N 0, x⊤ ϕ(U ) ϕ(U ) . x x x E , U ∼ N 0, , and V ∼ N 0, 1 1 1 2 2 Z d 1 d 2 n √ Then the NNGP G = ȳ/ dout is, explicitly, a Rdout -valued random variable with i.i.d. components centered Gaussian with variance: h i ΣG = Iddout E ψ(V, U2 )⊤ ψ(V, U2 ) . Note that this agrees with the choice of covariance in the infinite-width execution (1.3). Netsor program 3 Time-unrolled recurrent architecture Input: x1 , x2 : H(d) Input: Wx : A(n, d) Input: Wh : A(n, n) Input: Wo : A(dout , n) 1 g1 = Wx x1 : G(n) 2 h1 = ϕ(g1 ) : H(n) 3 g3 = Wh h1 : G(n) 4 g2 = Wx x2 : G(n) 5 h2 = ψ(g2 , g3 ) : H(n) 6 y = Wo h2 : G(dout )
▷ MatMul ▷ NonLin ▷ MatMul ▷ MatMul ▷ NonLin ▷ MatMul
QUANTITATIVE GAUSSIAN-PROCESS LIMITS OF TENSOR PROGRAMS
input time 1
x1 : H(d)
Wx
g1 : G(n)
ϕ
h1 : H(n)
Wh
9
g3 : G(n)
ψ input time 2
x2 : H(d)
Wx
g2 : G(n)
ψ
h2 : H(n) Wo
y : G(dout ) output
Fig. 4. Time-unrolled recurrent architecture. The same input weight Wx is reused across time, and Wh maps the previous hidden state h1 into the next preactivation g3 . The postactivation is defined as h2 = ψ(g2 , g3 ). The recurrent dimension n is an internal hidden width. Two-input shallow network. To analyze the covariance structure of the outputs given multiple data points, Figure 5 depicts a shallow network evaluated on two distinct inputs x1 and x2 . The corresponding Netsor program 4 is obtained by duplicating the lines of the shallow network program 1. Theorem 1.2 yields the expected estimate: 1 y2 c y √ ,√ ,G ≤ √ , (2.4) Wp n dout dout where G = G 1 , G 2 with G j the NNGP for the shallow architecture applied to the inputs xj , for j = 1, 2. Netsor program 4 Two-input shallow network Input: x1 , x2 : H(d) Input: W(0) : A(n, d) Input: W(1) : A(dout , n) 1 g1 = W(0) x1 : G(n) 2 g2 = W(0) x2 : G(n) 3 h1 = ϕ(g1 ) : H(n) 4 h2 = ϕ(g2 ) : H(n) 5 y1 = W(1) h1 : G(dout ) 6 y2 = W(1) h2 : G(dout )
▷ MatMul ▷ MatMul ▷ NonLin ▷ NonLin ▷ MatMul ▷ MatMul
2.2. Extensions. The Netsor+ language introduced in [Yan21] allows rather general scalar-valued operations (Moment) yielding additional parameters that can be used in functions of NonLin lines. This extension covers in particular attention layers, where scalar variables are used mainly to store empirical kernels, and these kernels are then fed into scalar nonlinearities such as the softmax. In this section, we describe a similar extension, denoted here by NetsorK , that can be obtained with minimal modification of our arguments, but is still sufficient to cover attention and normalization of layers. Our main result Theorem 1.2 extends naturally to this extended language, effectively covering attention-based architectures.
10
ANDREA AGAZZI, ELOY MOSIG GARCÍA, AND DARIO TREVISAN
input 1
x1 : H(d)
W(0)
g1 : G(n)
ϕ
h1 : H(n)
W(1)
y1 : G(dout ) output 1
input 2
x2 : H(d)
W(0)
g2 : G(n)
ϕ
h2 : H(n)
W(1)
y2 : G(dout ) output 2
Fig. 5. A shallow network evaluated on two inputs. The same A-vars are reused for the two inputs, represented by repeated edge labels W(0) and W(1) . The width n is hidden and is common to both branches. NetsorK programs. A NetsorK program is a Netsor program, augmented with Cvariables: in the finite-width execution they are real-valued random variables, while in the infinite-width execution they are deterministic constants. In addition to the usual MatMul and NonLin instructions, we allow the following operations: • Kernel variables. If h, h′ : H(n), then (Kernel)
s := h⋆ h′ : C.
• Nonlinearities with scalar parameters. If ϕ : Rk × Ru → R, if h1 , . . . , hk : H(n) and s1 , . . . , su : C, then (NonLin+ )
h := ϕ(h1 , . . . , hk ; s1 , . . . , su ) : H(n),
where ϕ is to be interpreted as acting componentwise on the H-vars. • Scalar maps. If ψ : Ru → R and s1 , . . . , su : C, then (Scal)
s := ψ(s1 , . . . , su ) : C.
Extended program execution. The execution rules for NetsorK naturally extend those of Netsor. Keeping MatMul and NonLin unchanged, the execution of a Kernel line s = h⋆ h′ is defined in finite-width interpretation as s = n1 h⊤ h′ , and in infinite-width as s̄ = h̄∗ h̄′ . The Scal and NonLin+ instructions are executed as componentwise compositions. Moreover, we assume that input C-variables, if present, are to be executed as deterministic constants. Hence, every C-var generated by a NetsorK program is a measurable function of the previously generated H-vars. In particular, scalar variables do not introduce additional observations of the weights beyond those already encoded by the filtration. Consequently, the conditional Gaussian formulas for MatMul lines remain identical to those of the base Netsor language. We impose the following regularity assumption on NetsorK programs: functions ϕ appearing in NonLin+ lines satisfy local Lipschitz estimate with linear growth. Namely, there exists a constant C < ∞ such that, for all x, y ∈ Rk and s, t ∈ Ru , (2.5)
|ϕ(x; s) − ϕ(y; t)| ≤ C(1 + ∥s∥ + ∥t∥) ∥x − y∥ + C(1 + ∥x∥ + ∥y∥) ∥s − t∥ .
Remark 2.1 (Relation with Yang’s Moment instruction). The distinction between NetsorK and the Netsor+ language of [Yan20b] is mainly analytic rather than algebraic. In Yang’s formulation, scalar variables are generated by a general Moment instruction of the form n
(Moment)
s :=
1X ψ gi1 , . . . , gik ; s1 , . . . , su : C. n i=1
QUANTITATIVE GAUSSIAN-PROCESS LIMITS OF TENSOR PROGRAMS
11
Whenever the componentwise quantity can itself be represented as the i-th coordinate of an H(n)-var, the corresponding moment can be written as a kernel against the constantone vector, 1n := (1, . . . , 1) ∈ Rn . Thus, at a purely formal level, Kernel together with scalar operations can reproduce any Moment instructions. The reason for isolating the restricted language NetsorK is that the quantitative argument requires stability estimates compatible with the norms controlled by the main theorem. Kernel variables are controlled directly by Corollary 1.3, while subsequent nonlinear maps are required to satisfy the Lipschitz condition (2.5). By contrast, a fully general Moment instruction with a locally Lipschitz integrand of polynomial growth would require additional empirical high-moment estimates. Hence, the restriction to Kernel variables serves to state a quantitative theorem under transparent and checkable analytic assumptions. Extended program graphs. For NetsorK programs, the directed program graph must also record scalar variables. In addition to H- and G-nodes, we include nodes of type C. A Kernel line s = h⋆ h′ contributes directed edges h −→ s, h′ −→ s, which we draw as dotted arrows with the label ⋆. A Scal line s = ψ(s1 , . . . , su ) contributes directed edges si → s, labelled by ψ. Finally, a scalar-parametric function contributes edges both from the H-variables and from the scalar variables to the output h. The definition of a hidden width expands accordingly: a dimension n is a hidden width if there exists an H-var of shape n produced by a MatMul line whose downstream descendants are later contracted as the parent of a further MatMul line, or within a Kernel line s = h⋆ h′ . We continue to denote by nhidd the collection of such widths. This convention is important for attention layers. The query and key vectors, denoted in the program by q and k respectively, have dimension u ∈ nhidd , and the pre-softmax scores wij = qi⋆ kj ,
for each i, j ≤ u,
are obtained by contracting over that dimension. Hence u is a hidden width even if the corresponding variables are not later used as parents of another MatMul line. Dimensions not belonging to nhidd are either considered input nin , hence structural, or output nout . Example 2.2 (Single-head attention). The following program describes a single layer of single-head attention with deterministic kernel scaling. The nonlinear functions involved are x−µ LayerNormε (x; µ, σ) := √ , σ2 + ε
exp(zj ) . m=1 exp(zm )
SoftMax(z1 , . . . , zu )j := Pu 1
u
i1
iu
(2.5) is tailored to the attention output operation (v , . . . , v ; a , . . . , a ) 7−→ PuCondition ij v j , which is not globally Lipschitz jointly in (v, a), but satisfies the locally Lipa j=1 schitz bound with C = 1. For layer normalization LayerNormε , the nonlinearity also satisfies (2.5) with C depending on ε. For readability, Figure 6 displays only the dependencies associated with a fixed query token i and key/value token j; the full attention layer is obtained by repeating the same operation for all i, j ≤ u. Note that query and key vectors have dimension n, and the pre-softmax score is obtained by contracting over that dimension. Hence, n is a hidden width even if the variables are not later used as parents of another MatMul instruction. Successive iterations of the attention block in Algorithm 5 make LayerNorm act on noninput H-vars. In this case the statistics µ and σ can be computed via the new instructions
12
ANDREA AGAZZI, ELOY MOSIG GARCÍA, AND DARIO TREVISAN query
qi : G(n) input tokenW
xi : H(d)
⋆
q
Wk
wij : C ⋆
kj : G(n)
score
SoftMax ij
Wv
a :C P
attention weight
ij j ja v
j
v : G(n)
yi : H(n) attention output
Fig. 6. Schematic program graph for a single-head attention layer in NetsorK . Solid arrows are MatMul operations labelled by A-vars, dotted arrows indicate kernel/scalar operations, and dashed arrows indicate scalarparametric nonlinearities. The width n is produced by the projections and contracted in the kernel scores wij = qi⋆ kj . Netsor program 5 A single attention layer in NetsorK Input: {xi }ui=1 : H(d) Input: {µi , σ i }ui=1 : C Input: Wq , Wk , Wv : A(n, d)
▷ input LayerNorm statistics, with σ i > 0 ▷ 1. Layer normalization
for i = 1, . . . , u do 1 x̃i := LayerNormε (xi ; µi , σ i ) : H(d) end for
▷ NonLin+ ▷ 2. Projections
for i = 1, . . . , u do 2 qi := Wq x̃i : G(n) 3 ki := Wk x̃i : G(n) 4 vi := Wv x̃i : G(n) end for
▷ MatMul ▷ MatMul ▷ MatMul ▷ 3. Pre-softmax scores
for i, j = 1, . . . , u do 5 wij := qi⋆ kj : C end for
▷ Kernel ▷ 4. Softmax
for i = 1, . . . , u do for j = 1, . . . , u do 6 aij := SoftMax(wi1 , . . . , wiu )j : C end for end for
▷ Scal
▷ 5. Attention output for i = 1, .P . . , u do u ij j 7 yi := j=1 a v : H(n) end for
▷ NonLin+
with: (Kernel)
µ = h⋆ 1n
: C,
(NonLin+ )
hc = h − µ
: H(n),
(Kernel)
σ
2
= h⋆c hc
: C,
(Scal)
σ=
σ2
: C.
√
QUANTITATIVE GAUSSIAN-PROCESS LIMITS OF TENSOR PROGRAMS
13
Remark 2.3. With the 1/n normalization in the Kernel rule, the attention scores have deterministic infinite-width limits. This is the√law-of-large-numbers scaling. A different scaling, for instance the fluctuation scaling 1/ n for centered kernels as the one used in [SKI26], would lead to Gaussian limits for the scalar variables and requires a quantitative CLT rather than the kernel LLN used here. Extended convergence result. With the above formalization, we state the extension of Theorem 1.2 to NetsorK programs, under an explicit non-degeneracy assumption for G-vars in the infinite-width limit. Theorem 2.4. Consider a NetsorK program such that nonlinearities in NonLin, Scal are Lipschitz and those in NonLin+ satisfy the local Lipschitz bound (2.5). Let h(ℓ) ℓ=1,...,H , s(j) j=1,...,S , h̄(ℓ) ℓ=1,...,H , s̄(j) j=1,...,S , be the finite-width and infinite-width executions of its H-vars and C-vars, where H and S denote, respectively, the number of H-vars and C-vars defined in the program, and assume that its G-vars are jointly non-degenerate Gaussians in the infinite-width execution. Then, for every p ≥ 1, there exists a constant c < ∞, depending only on p and on the structural constants of the program, such that ! ! (ℓ) (ℓ) h̄ h , √ , s(j) , s̄(j) j=1,...,S Wp √ nℓ nℓ j=1,...,S ℓ=1,...,H ℓ=1,...,H (2.6) X 1 √ . ≤c m m∈n hidd
Remark 2.5. In the basic Netsor setting we used an equivalence reduction to remove degenerate MatMul lines. For NetsorK this reduction is more delicate: a scalar kernel variable is deterministic in infinite-width but random at finite width, so a linear relation that holds identically in the infinite-width execution need not hold identically in the finitewidth execution. Specifically, the pointwise lifting of infinite-width linear identities to finite-width linear identities utilized in the proof of Proposition 3.11 fails when these identities depend on the fluctuations of empirical scalar kernels. For this reason we state the extension under an explicit non-degeneracy assumption for G-vars in the infinite-width execution. The proof of this theorem requires an extension of our inductive line-by-line coupling arguments to encompass scalar variable executions and is given in Subsection 4.3. 3. Technical lemmas In this section we collect some technical results that are used in the proof of Theorem 1.2. We begin by setting some notation: ∥ · ∥, ∥ · ∥op and ∥ · ∥Lp denote the Euclidean norm, operator norm and Lp norm, respectively. Unless specified, the norm inner to the Lp norm is the Frobenius norm. As anticipated in the introduction, c denotes a positive constant that may change from line to line. Let S, T in N. For any block matrix A B M= ∈ R(S+T )×(S+T ) , C D we define its generalized Schur complement (3.1)
M/D := A − BD+ C,
where the superscript + denotes the Moore-Penrose pseudoinverse (3.2)
M + = lim(M ⊤ M + δIdT )−1 M ⊤ , δ↓0
14
ANDREA AGAZZI, ELOY MOSIG GARCÍA, AND DARIO TREVISAN
for any M ∈ RS×T . If the block matrix M as above is symmetric positive semidefinite, then so is M/D. This extension matches the classical Schur complement when D is invertible, and is the relevant form for Gaussian conditioning with singular covariance matrices, see Equation (3.7) below. We also recall Greville’s formula for the pseudoinverse of a rank-one update of a matrix [UK97; Gre60]. Let B ∈ RS×T be partitioned as B = [A, a], where A consists of the first T − 1 columns of B and a is the last column. Then the pseudoinverse of B is given by + A − A+ ab + (3.3) B = , b where (3.4)
( (a − AA+ a)+ b= (1 + a⊤ (A⊤ A)+ a)−1 a⊤ (A⊤ A)+
if (IdS − AA+ )a ̸= 0, otherwise.
This formula adapts the well known Woodbury-Sherman-Morrison [SM50] formula to rectangular matrices and their pseudoinverses. Furthermore, the generalized Schur complement, the pseudoinverse, and Greville’s formula can be understood from an operator-theoretic point of view, see e.g. [And05]. Whenever B represents a linear operator between two finite-dimensional Hilbert spaces H and H′ , its pseudoinverse coincides with the matrix pseudoinverse as long as orthonormal bases of H and H′ are chosen to represent the operator, where A⊤ in (3.2) becomes the adjoint operator to A, defined as ⟨Ax, y⟩H′ = ⟨x, A⊤ y⟩H ,
for all x ∈ H and y ∈ H′ ,
whenever A is a linear operator from H to H′ . Note that the pseudoinverse explicitly depends on the inner products of H and H′ . Under these considerations, Greville’s formula (3.3) can be applied to finite-dimensional Hilbert spaces. 3.1. Random variables and moment bounds. Given p ≥ 1 and a probability space (Ω, A, P), we write X ∈ Lp (Ω, A, P; RS ) to denote a random variable X : Ω → RS such that ∥X∥Lp := E [∥X∥p ]1/p < ∞. Since we are mostly interested in convergence in law, we often do not specify the underlying probability space for random variables, and also write for brevity X ∈ RS instead of law X : Ω → RS , X = Y to denote that two random variables X and Y have the same law. The following elementary concentration result will be paramount for proving Theorem 1.2. Lemma 3.1. Let p ≥ 1, let X be a random variable with values in RS and let f : RS → RT be differentiable at x̄ ∈ RS . Then, for some constant c = c(f, x̄) < ∞ it holds (3.5)
∥f (X) − f (x̄)∥Lp ≤ c (∥X − x̄∥Lp + (∥f (X)∥L2p + ∥f (x̄)∥) ∥X − x̄∥L2p ) .
Proof. By differentiability of f at x̄, there exists ε > 0 and M < ∞ such that, if ∥x − x̄∥ ≤ ε, then ∥f (x) − f (x̄)∥ ≤ M ∥x − x̄∥ . Consider the event A =: {∥X − x̄∥ ≤ ε}, so that by Markov inequality: ∥1Ac ∥L2p ≤ ε−1 ∥X − x̄∥L2p . Then, ∥f (X) − f (x̄)∥Lp ≤ ∥(f (X) − f (x̄)) 1A ∥Lp + ∥(f (X) − f (x̄)) 1Ac ∥Lp ≤ M ∥X − x̄∥Lp + ∥(∥f (X)∥ + ∥f (x̄)∥) 1Ac ∥Lp
≤ M ∥X − x̄∥Lp + (∥f (X)∥L2p + ∥f (x̄)∥) ∥1Ac ∥L2p ≤ c (∥X − x̄∥Lp + (∥f (X)∥L2p + ∥f (x̄)∥) ∥X − x̄∥L2p )
QUANTITATIVE GAUSSIAN-PROCESS LIMITS OF TENSOR PROGRAMS
with c := max M, ε−1 .
15
□
When considering matrix valued random variables W ∈ RS×T , we will say that the covariance ΣW is separable if there exist L ∈ RS×S , R ∈ RT ×T such that ΣW = L ⊗ R. We will also rely on Rosenthal’s inequality to bound the moments of sums of independent, centered random variables [Ros70]. Let p ≥ 2, and let X1 , . . . , Xn ∈ Lp (Ω; RS ) be such random variables. Then, there exists a constant cp depending only on p such that !p/2 p n n n X X X (3.6) ≤ cp max ∥Xi ∥pLp , Xi ∥Xi ∥2L2 . p i=1
L
i=1
i=1
We now record a uniform moment bound for the infinite-width execution of Netsor programs. Lemma 3.2 (Uniform infinite-width moment bounds). Consider a Netsor program such that all activations ϕ appearing in NonLin lines have polynomial growth. Then, for every p ≥ 1, there exists a constant c < ∞, depending only on p and on the structural constants of the program, such that for every H-var h : H(n) generated by the program, every coordinate i ≤ n, and every infinite-width execution, E |h̄i |p ≤ c. In particular, c is independent of the hidden widths. Proof. We prove the claim by induction over the program lines. Since the program is finite, it is enough to show that each step propagates bounds depending only on the previous bounds and on the structural data of the program, but not on the widths. Input H-vars are deterministic. Hence, if x : H(d) is an input, then |x̄i |p ≤ max max |xi |p , x∈X 1≤i≤d
which is finite and independent of the hidden widths. Assume now that the bound has been proved for all H-vars available before a given line. If the line is of type MatMul, say g = Wh with h : H(m) and g : G(n), then by the infinite-width execution rule each coordinate ḡi is centered Gaussian with variance m 1 X E[|h̄r |2 ]. Var(ḡi ) = h̄∗ h̄ = m r=1
By the induction hypothesis with exponent p = 2, Var(ḡi ) ≤ c, uniformly in m, n and in i. Therefore, by finiteness of Gaussian moments, there exists a positive constant cp with: p/2 ≤c E[|ḡi |p ] ≤ cp Var(ḡi ) This gives the desired uniform bound for the new G-var. If the line is of type NonLin, say h = ϕ(h1 , . . . , hk ), then, componentwise, h̄i = ϕ(h̄1,i , . . . , h̄k,i ). By polynomial growth of ϕ, there exist constants c, q < ∞ such that ! k X |ϕ(x1 , . . . , xk )| ≤ c 1 + |xr |q . r=1
16
ANDREA AGAZZI, ELOY MOSIG GARCÍA, AND DARIO TREVISAN
Hence, for every p ≥ 1, p
|h̄i | ≤ cp
1+
k X
! pq
|h̄r,i |
.
r=1
Taking expectations and applying the induction hypothesis with exponent pq gives E[|h̄i |p ] ≤ c, where the constant depends only on the previous moment bounds and on the polynomialgrowth constants of ϕ, but not on the widths. □ 3.2. Conditional laws. The proof of Theorem 1.2 relies on tracking the laws of the variables introduced by the program, conditionally on the σ-algebra generated by the previously defined H-vars. The key point is that in both finite-width and infinite-width executions, every newly created G-var at a MatMul line is conditionally Gaussian, and its conditional mean/covariance can be computed recursively from earlier lines. In this subsection we introduce such recursive formulas in Lemmas 3.5 and 3.6, but first we recall some definitions and results from theory of Gaussian variables. Gaussian variables. Let S, T be natural numbers. If (X, Y ) is jointly Gaussian with values in RS × RT , with mean µX,Y = (µX , µY ) and covariance ΣX ΣXY Σ(X,Y ) = ∈ R(S+T )×(S+T ) , ΣY X ΣY then the marginal Y is Gaussian and, conditionally on Y = y, the law of X is Gaussian with (3.7) µX|Y =y = µX + ΣXY Σ+ Y (y − µY ),
ΣX|Y =y = Σ(X,Y ) /ΣY = ΣX − ΣXY Σ+ Y ΣY X ,
where Σ+ Y denotes the Moore-Penrose pseudoinverse and Σ(X,Y ) /ΣY is the generalized Schur complement (3.1). In particular, y 7→ µX|Y =y is affine and y 7→ ΣX|Y =y is constant. These properties also characterize joint Gaussianity [Rao73]: Lemma 3.3. Let (X, Y ) be a random variable in RS × RT such that Y is Gaussian and, for a.e. y ∈ RT , the conditional law of X given Y = y is Gaussian with y 7→ µX|Y =y affine and y 7→ ΣX|Y =y constant. Then (X, Y ) is jointly Gaussian. The formulas (3.7) extend verbatim to matrix-valued Gaussian variables, once covariances are interpreted as 4-index tensors. Concretely, fix S1 , S2 , T1 , T2 ∈ N and let (X, Y ) be jointly Gaussian with values in RS1 ×S2 ×RT1 ×T2 . Viewing X and Y as vectors of lengths S1 S2 and T1 T2 respectively, all standard vector-Gaussian theory applies, Regarding X and Y as matrices the cross-covariance ΣXY = E (X − µX ) ⊗ (Y − µY ) ∈ RS1 ×S2 ×T1 ×T2 is now a 4-index tensor, where ⊗ denotes the Kronecker product, that is, the outer product of the two matrices regarded as vectors (i.e. (ΣXY )i1 i2 j1 j2 = E[(X − µX )i1 i2 (Y − µY )j1 j2 ]). The products appearing in (3.7) are then tensor contractions over the Y -indices: given additionally Z with values in RU1 ×U2 , the product ΣXY ΣY Z contracts the (T1 , T2 ) indices shared by ΣXY ∈ RS1 ×S2 ×T1 ×T2 and ΣY Z ∈ RT1 ×T2 ×U1 ×U2 , yielding an element of RS1 ×S2 ×U1 ×U2 . Equivalently, identifying each matrix space with the corresponding vector space via lexicographic ordering of indices, the tensor contraction coincides with the usual matrix product on the vectorized representatives, so (3.7) holds. In particular, given a Gaussian random matrix W ∈ Rn×m and a deterministic vector h ∈ Rm , the product W h is Gaussian with µW h = µW h,
ΣW h = ΣW (hh⊤ ) ∈ Rn×n ,
QUANTITATIVE GAUSSIAN-PROCESS LIMITS OF TENSOR PROGRAMS
17
where the right-hand side uses the tensor contraction just described (contracting the mdimensional index of ΣW ∈ Rn×m×n×m against hh⊤ ∈ Rm×m , yielding Rn×n ). The joint variable (W, W h) ∈ R(n+1)×m is Gaussian, and (3.7) gives µW |W h=g = µW + ΣW,W h Σ+ W h (g − µW h),
ΣW |W h=g = ΣW − ΣW,W h Σ+ W h ΣW h,W .
In the separable covariance case ΣW = L ⊗ R (where L ∈ Rn×n , R ∈ Rm×m ), the constituent tensors factorize explicitly: ΣW,W h = L ⊗ (Rh),
ΣW h = L (h⊤ Rh),
ΣW h,W = L ⊗ (h⊤ R).
Substituting into the conditioning formulas yields (3.8) (3.9)
µW |W h=g = µW + (g − µW h) (h⊤ Rh)+ (Rh)⊤ , ΣW |W h=g = L ⊗ R − Rh (h⊤ Rh)+ h⊤ R .
In particular, the separable structure of ΣW is preserved under conditioning on W h = g. Recall that since ΣW is a covariance tensor R is positive semidefinite. Let us notice that, since h⊤ Rh ≥ 0 is scalar, the pseudoinverse (h⊤ Rh)+ is actually the inverse if h⊤ Rh > 0 and 0 otherwise. This is quite expected since in the latter case W h is constant and conditioning on W h = g leaves the law of W unchanged. Let us remark that the formulas above do not require invertibility of L. Indeed, if L is not invertible, terms of the form LL+ L appear in the Schur complement, but it holds LL+ L = L. Finite-width execution. To describe the conditional law of G-vars, it is useful to describe first the conditional law of A-vars. Indeed, each A-var W : A(n, m) is initialized as a 1 Gaussian random matrix W ∈ Rn×m with separable covariance ΣW = m Idn ⊗ Idm . Given a Netsor program, denote by Hℓ the set of H-vars defined up to line ℓ. For the σ-algebra Fℓ := σ(Hℓ ), conditioning amounts to linear constraints of the form g = W h, where g, h ∈ Hℓ . Hence, the conditional law of W remains Gaussian, with updates given by applications of (3.8) and (3.9) at each MatMul operation. In particular, since W : A(n, m) has separable initial covariance, then we also have ΣW |Fℓ = LW |Fℓ ⊗ RW |Fℓ ,
ℓ ≥ 0,
with suitable LW |Fℓ and RW |Fℓ , that we describe below in (3.12). We start with explicitly writing the update rules for MatMul lines. If line ℓ is a MatMul operation g := Wh with finite-width execution g = W h, then h is Fℓ−1 measurable µW |Fℓ , LW |Fℓ and RW |Fℓ are special cases of (3.8) and (3.9). These read as follows + (3.10) µW|Fℓ = µW|Fℓ−1 + g − µW|Fℓ−1 h h⊤ RW|Fℓ−1 h h⊤ RW|Fℓ−1 . The left part of the covariance tensor remains unchanged: 1 (3.11) LW|Fℓ = LW|Fℓ−1 = Idn , m while the right part satisfies: + (3.12) RW|Fℓ = RW|Fℓ−1 − RW|Fℓ−1 h h⊤ RW|Fℓ−1 h h⊤ RW|Fℓ−1 . Moreover different A-vars remain independent. We obtain as a consequence the conditional law of G-vars. Indeed, these are defined only via MatMul operations: therefore, if line ℓ is as above a MatMul operation g := Wh with finite-width execution g = W h, then h is Fℓ−1 -measurable and, conditionally on Fℓ−1 , g is Gaussian with (3.13)
µg|Fℓ−1 = µW |Fℓ−1 h,
18
ANDREA AGAZZI, ELOY MOSIG GARCÍA, AND DARIO TREVISAN
and (3.14)
Σg|Fℓ−1 = LW |Fℓ−1
h⊤ RW |Fℓ−1 h .
We notice that NonLin lines do not enlarge the σ-algebras, hence Fℓ = Fℓ−1 for such lines, hence conditional laws remain unchanged in such case. Thus, conditional laws are algebraically characterized by following the update rules (3.10) and (3.12) along the program lines. However, some structural properties are more evident by providing also a geometric description. To this aim, for any program line ℓ, we introduce a Gram matrix associated to each W : A(n, m), defined in the following way. We collect all the MatMul lines before (possibly including) line ℓ, in which W appears: g(1) = Wh(1) , . . . , g(k) = Wh(k) , and we collect such G-vars and H-vars in the (random) matrices, respectively in Rn×k and Rm×k , (ℓ) (ℓ) , , HW := h(i) (3.15) GW := g (i) 1≤i≤k
1≤i≤k
so that the Gram matrix (associated to the h-vars for W up to line ℓ) in Rk×k , is given by (ℓ),⊤ (ℓ) (i),⊤ (j) , (3.16) HW HW = h h 1≤i,j≤k and the orthogonal projection on the (random) subspace in Rm , spanned by h(i) 1≤i≤k , is given by (ℓ) (ℓ),⊤ (ℓ) + (ℓ),⊤ (ℓ) (ℓ) + (3.17) HW HW HW HW = HW HW ∈ Rm×m . With this notation, we will prove the structure lemma 3.5 for conditional laws. Example 3.4. Let W : A(2, 3), and let its finite-width execution be a Gaussian matrix W ∈ R2×3 with independent entries of variance 1/3. Suppose that, before line ℓ, the same A-var has already appeared in two MatMul lines g1 = Wh1 ,
g2 = Wh2 ,
h1 , h2 ∈ R3 ,
and assume, for simplicity, that the finite executions h1 , h2 are linearly independent. After conditioning on Fℓ−1 , the vectors h1 , h2 , g1 , g2 are fixed. Using the matrices (ℓ−1) (ℓ−1) (ℓ−1) + (ℓ−1) (ℓ−1) −1 (ℓ−1) HW = (h1 , h2 ), GW = (g1 , g2 ), HW = (HW )⊤ HW (HW )⊤ . we can decompose the finite execution of a new MatMul line g = Wh, into the Fℓ−1 measurable part (the one already “tested” by W ) and its orthogonal residual: (ℓ−1) (ℓ−1) + (ℓ−1) (ℓ−1) + h = HW HW h + Id3 − HW HW h (ℓ−1)
= HW α + r = α1 h1 + α2 h2 + r, (ℓ−1) +
α = HW
h = (α1 , α2 )⊤ ,
(ℓ−1)
RW |Fℓ−1 = Id3 − HW
where (ℓ−1) +
HW
,
r = RW |Fℓ−1 h.
Note that here RW |Fℓ−1 is the projection on the orthogonal complement of span{h1 , h2 }. Applying W to the above formula allows a similar decomposition of g g = W h = α1 W h1 + α2 W h2 + W r = α1 g1 + α2 g2 + W r. where the first two terms are Fℓ−1 -measurable (i.e., given under the conditioning), while W r is fresh Gaussian noise. Based on this decomposition, the conditional distribution of g is given by I2 ⊤ (ℓ−1) (ℓ−1) + (3.18) g | Fℓ−1 ∼ N GW HW h, h RW |Fℓ−1 h . 3
QUANTITATIVE GAUSSIAN-PROCESS LIMITS OF TENSOR PROGRAMS
19
This computation foreshadows the content of the structure lemmas (Lemma 3.5 and Lemma 3.6) below, characterizing in full generality the information that previous uses of the same A-var provide on how W acts on span{h1 , h2 }. As outlined above, these lemmas show that a new product W h is predictable on the projection of h onto this span, but remains Gaussian through the residual component orthogonal to it, resulting in the formula (3.18). If the residual vanishes, then h ∈ span{h1 , h2 }, so W h = α1 g1 + α2 g2 . In this case, the execution of the program can be equivalently rewritten in nondegenerate form by replacing the MatMul line with a NonLin line, as formalized in Lemma 3.9 below. Lemma 3.5 (Structure lemma for finite-width execution). Consider an A-var W : A(n, m) in a Netsor program and let Fℓ = σ(Hℓ ) be the σ-algebra generated by the H-vars defined up to line ℓ by the finite-width execution of the program. With the notation (i) introduced above, RW |Fℓ is the projector on the orthogonal subspace to the span of h 1≤i≤k , (ℓ)
(ℓ) +
RW |Fℓ = Idm − HW HW , and (ℓ)
(ℓ) +
µW |Fℓ = GW HW . If moreover line ℓ is MatMul g = Wh, then, conditionally upon Fℓ−1 , g is Gaussian with mean (ℓ−1) (ℓ−1) + µg|Fℓ−1 = GW HW h, and covariance Idn ⊤ (ℓ−1) (ℓ−1) + h , (3.19) Σg|Fℓ−1 = h h − h⊤ HW HW m where h ∈ Fℓ−1 is the finite-width execution of h. In particular, Lemma 3.5 implies that ∥RW |Fℓ ∥op ≤ 1 whenever the ℓ-th line is a MatMul instruction. Proof. We reason by induction on ℓ. Let k be the number of MatMul lines involving W up to line ℓ. If the ℓ-th line is NonLin, then there is nothing to prove. Hence, assume that the ℓ-th line is g(k+1) = Wh(k+1) . Write h = h(k+1) and g = g (k+1) for simplicity, and µ(ℓ) = µW|Fℓ , R(ℓ) = RW|Fℓ for each ℓ ≤ L. For the sake of clarity and since W is fixed, (ℓ)
(ℓ)
we will also omit the dependence of W in the block matrices HW and GW . We proceed to prove the claim for µ(ℓ) and R(ℓ) separately. Proof for R(ℓ) : We begin by studying the case h⊤ R(ℓ−1) h = 0. This implies that h is in the span of + {h(1) , . . . , h(k) } and h⊤ R(ℓ−1) h = 0, yielding +
R(ℓ) = R(ℓ−1) = Idm − H (ℓ−1) H (ℓ−1) . Note that since the column spaces of {h(1) , . . . , h(k) } and {h(1) , . . . , h(k) , h} coincide, we + + have H (ℓ) H (ℓ) = H (ℓ−1) H (ℓ−1) , so the claim holds in the degenerate case. We now consider the case h⊤ R(ℓ−1) h > 0. Denote by P the orthogonal complement to the span of {h(1) , . . . , h(k) } in Rm . Recalling the update (3.12), R(ℓ) = R(ℓ−1) Idm − hh⊤ R(ℓ−1) (h⊤ R(ℓ−1) h)−1 , is immediate under the induction hypothesis that R(ℓ) is symmetric and idempotent. Interpreting R(ℓ) as a linear operator on Rm , we have im R(ℓ) ⊆ im R(ℓ−1) = P . Moreover, for every x ∈ Rm , h⊤ R(ℓ) x = h⊤ R(ℓ−1) x − h⊤ R(ℓ−1) h(h⊤ R(ℓ−1) h)−1 h⊤ R(ℓ−1) x = 0,
20
ANDREA AGAZZI, ELOY MOSIG GARCÍA, AND DARIO TREVISAN
so that im R(ℓ) ⊆ P . Conversely, if y ∈ P , then R(ℓ−1) y = y and h⊤ y = 0, hence R(ℓ) y = R(ℓ−1) y − R(ℓ−1) h(h⊤ R(ℓ−1) h)−1 h⊤ R(ℓ−1) y = y. Thus P ⊆ im R(ℓ) . The uniqueness of the orthogonal projection implies the first claim. Proof for µ(ℓ) : Again, we study first the case h⊤ R(ℓ−1) h = 0. In this case, the update (3.10) reduces to µ(ℓ) = µ(ℓ−1) . In effect, note that G(i) = W H (i) for each i ≤ ℓ. Therefore, by reasoning as in the previous case: +
+
+
+
G(ℓ) H (ℓ) = W H (ℓ) H (ℓ) = W H (ℓ−1) H (ℓ−1) = G(ℓ−1) H (ℓ−1) . This proves the claim in the degenerate case. Let’s now assume that h⊤ R(ℓ−1) h > 0. The recursion (3.10) can be written as a rank-one correction µ(ℓ) = µ(ℓ−1) + g − µ(ℓ−1) h qℓ⊤ , where −1 qℓ⊤ := h⊤ R(ℓ−1) h h⊤ R(ℓ−1) . By induction hypothesis, +
+
µ(ℓ) = G(ℓ−1) H (ℓ−1) + (g − G(ℓ−1) H (ℓ−1) h)qℓ . On the other hand, since G(ℓ) = G(ℓ−1) , g and H (ℓ) = H (ℓ−1) , h , and since h is not + in the span of h(i) 1≤i≤k by assumption, we have (Idm − H (ℓ−1) H (ℓ−1) )h ̸= 0, so the first branch of Greville’s formula (3.3) can be applied with A = H (ℓ−1) and a = h. The corresponding row vector b ∈ R1×m is: +
b = (h − H (ℓ−1) H (ℓ−1) h)+ = (R(ℓ−1) h)+ . Since R(ℓ−1) is an orthogonal projector and R(ℓ−1) h ̸= 0, the pseudoinverse is given by: b=
(R(ℓ−1) h)⊤ = qℓ⊤ . h⊤ R(ℓ−1) h
Then, Greville’s formula (3.3) yields H
(ℓ) +
=
+ H (ℓ−1) (Idm − hqℓ ) , qℓ⊤
Therefore, (3.20)
+
+
G(ℓ) H (ℓ) = G(ℓ−1) H (ℓ−1) (Idm − hqℓ⊤ ) + gqℓ⊤ +
+
(3.21)
= G(ℓ−1) H (ℓ−1) + (g − G(ℓ−1) H (ℓ−1) h)qℓ⊤
(3.22)
= µ(ℓ) .
This concludes the proof of the second claim. Lastly, the formulas for the conditional mean and covariance of g follow directly from the first two claims and the general update formulas (3.13) and (3.14), substituting LW |Fℓ−1 = 1 □ m Idn in the latter.
QUANTITATIVE GAUSSIAN-PROCESS LIMITS OF TENSOR PROGRAMS
21
Infinite-width execution. In the infinite-width execution, A-vars are not realized as actual random matrices. Nevertheless, it is convenient to associate to each A-var W : A(n, m) a pair of barred parameters µW|F ℓ and ΣW|F ℓ playing the role of conditional mean and covariance. As in the finite-width execution, the covariance update preserves separability; we write ΣW|F ℓ = LW|F ℓ ⊗ RW|F ℓ . These quantities evolve by the same regression update maps as in the finite-width execution, but with all contractions computed using the infinite-width pairing ∗ defined in (1.4). Concretely, when line ℓ is MatMul g = Wh with executions ḡ, h̄ and F ℓ := σ(H̄ℓ ), we set + h̄∗ RW|F ℓ−1 h̄ h̄∗ RW|F ℓ−1 , µW|F ℓ = µW|F ℓ−1 + ḡ − µW|F ℓ−1 h̄ (3.23) and where LW|F ℓ = LW|F ℓ−1 = Idn , and + (3.24) RW|F ℓ = RW|F ℓ−1 − RW|F ℓ−1 h̄ h̄∗ RW|F ℓ−1 h̄ h̄∗ RW|F ℓ−1 . Such barred parameters can be thought as mere notational convenience, in order to describe the actual conditional law of G-vars: indeed, conditionally on F ℓ−1 , the new ḡ ∈ Rn is Gaussian with mean µḡ|F ℓ−1 = µW|F ℓ−1 h̄,
(3.25) and covariance
Σḡ|F ℓ−1 = LW|F ℓ−1 h̄∗ RW|F ℓ−1 h̄ .
(3.26)
We notice that also in the infinite-width execution, NonLin lines do not enlarge the σalgebras, hence F ℓ = F ℓ−1 for such lines, and conditional laws remain unchanged in such case. To provide a geometric description for the infinite-width conditional laws, we proceed analogously to the finite-width execution. We collect all the MatMul lines before (possibly including) line ℓ, in which W appears: g(1) = Wh(1) , . . . , g(k) = Wh(k) . We define the collections of their infinite-width executions ḡ (i) ∈ Rn and h̄(i) ∈ L2 (Rm ) via the formal operators: (ℓ) (ℓ) GW := ḡ (i) , H W := h̄(i) . 1≤i≤k
1≤i≤k
Analogously to (3.16), we define the Gram matrix in the infinite-width execution using the infinite-width pairing ∗ defined in Equation (1.4): (ℓ) (ℓ) H W ∗ H W := h̄(i)∗ h̄(j) . 1≤i,j≤k
In view of (1.3), the covariance matrix of the stacked vector (ḡ (i) )i=1,...,k is given by (ℓ)
(ℓ)
(H W ∗ H W ) ⊗ Idn . The orthogonal projection on the span of {h̄(i) }1≤i≤k in L2 (Rm ) is written compactly as: (ℓ) (ℓ) (ℓ) + (ℓ) (ℓ) (ℓ) + H W H W ∗H W HW ∗ = HW HW . A proof analogous to Lemma 3.5 yields the same structure for the infinite-width execution, with the only difference that all products are computed using the infinite-width pairing ∗ instead of the Euclidean inner product. In particular, Greville’s formula (3.3) still holds in this Hilbert space setting by choosing an orthonormal basis {ui }1≤i≤k of the span of h̄(i) 1≤i≤k in L2 (Rm ), as detailed at the beginning of the present section.
22
ANDREA AGAZZI, ELOY MOSIG GARCÍA, AND DARIO TREVISAN
Lemma 3.6 (Structure lemma for infinite-width execution). Consider an A-var W : A(n, m) in a Netsor program and let F ℓ be the σ-algebra generated by the H-vars defined up to line ℓ by the infinite-width execution of the program. With the notation introduced above, RW |F ℓ is the projector on the orthogonal subspace to the span of h̄(i) 1≤i≤k , (ℓ)
(ℓ),+
RW |F ℓ = Idm − H W H W
.
and (ℓ)
(ℓ),+
µW |F ℓ = GW H W
.
If moreover line ℓ is MatMul g = Wh, then, conditionally upon F ℓ−1 , ḡ is Gaussian with mean (ℓ−1)
µḡ|F ℓ−1 = GW
(ℓ−1),+
HW
h̄
and covariance (3.27)
(ℓ−1) (ℓ−1),+ Σḡ|F ℓ−1 = Idn h̄∗ h̄ − h̄∗ H W H W h̄ ,
where h̄ ∈ F ℓ−1 is the infinite-width execution of h. Remark 3.7 (Input parents in the basic Netsor language). In a basic Netsor program, if a MatMul line g = Wh has parent dimension m ∈ nin , then the finite execution h of h is obtained from deterministic input variables by ordinary NonLin operations. Hence its finite-width and infinite-width executions coincide. In particular, the corresponding Gram quantities are structural. This statement is used only in the proof for basic Netsor programs. In the extended NetsorK language, input-dimensional parents may also depend on scalar random variables, and the resulting finite-width fluctuations are controlled separately in the NonLin+ /Scal cases of the induction step. 3.3. Program equivalence and reduction to non-degeneracy. In this subsection we introduce a notion of equivalence between Netsor programs, and we use it to show that every program is equivalent to one in which all MatMul lines have nondegenerate conditional Gaussian outputs, in the infinite-width execution. This reduction is useful in the proof of the main theorem, since it allows us to work with densities and avoid technicalities related to degenerate Gaussians and non-invertible covariance matrices. Definition 3.8 (Program equivalence). Let T and T ′ be two Netsor programs of length M , and denote by h(1) , . . . , h(M )
and
h′(1) , . . . , h′(M )
their H-vars. We say that T and T ′ are infinite-width-equivalent if there exists a bijection π : {h(1) , . . . , h(M ) } −→ {h′(1) , . . . , h′(M ) } such that, writing π(h(i) ) = h′(π(i)) , the corresponding infinite-width executions satisfy law h̄(1) , . . . , h̄(M ) = h̄′(π(1)) , . . . , h̄′(π(M )) . We say that T and T ′ are equivalent if the same bijection also identifies the finite-width executions, namely law h(1) , . . . , h(M ) = h′(π(1)) , . . . , h′(π(M )) . We stress that in the definition above the bijection is not required to preserve the syntactic subtype: a G-var in one program may correspond to an H-var which is not a G-var in the other.
QUANTITATIVE GAUSSIAN-PROCESS LIMITS OF TENSOR PROGRAMS
23
Lemma 3.9 (Non-degenerate representative). Let T be a Netsor program whose NonLin functions are Lipschitz continuous. Then T is equivalent to a Netsor program T ′ such that, in the infinite-width execution of T ′ , every MatMul line has non-degenerate conditional Gaussian output: if g is generated at line ℓ of T ′ , then ḡ | F ℓ−1 is a non-degenerate Gaussian random vector. Remark 3.10 (Gram formulation of non-degeneracy). By the structure lemma for the infinite-width execution, the conclusion of Lemma 3.9 is equivalent to saying that no MatMul line contracts a parent which is already in the infinite-width span of the previous parents associated with the same A-var. More explicitly, if the line is g = Wh and the previous uses of the same A-var W have parents h(1) , . . . , h(k) , then the conditional covariance of ḡ is non-degenerate precisely when h̄∗ RW|F ℓ−1 h̄ > 0, or equivalently when the enlarged Gram matrix (ℓ)
(ℓ)
H W ∗H W has rank one more than
(ℓ−1)
(ℓ−1)
H W ∗H W . In particular, after applying Lemma 3.9, the Gram matrices for each A-var are invertible along the program. Before proving Lemma 3.9 we present an auxiliary result on the geometrical properties of the covariance tensor update (3.12). Proposition 3.11. Given a program with line ℓ of MatMul type g = Wh, the following are equivalent: (i) ḡ is Gaussian degenerate, conditionally upon F ℓ−1 ; (ii) h̄∗ RW|F ℓ−1 h̄ = 0; (iii) there exist (deterministic) (αi )i=1,...,k and previously defined G-vars (g(i) )i=1,...,k ∈ Hℓ−1 using the same A-var W such that, replacing the MatMul line g = Wh by a NonLin instruction: k X ′ g := αi g(i) i=1
yields an infinite-width-equivalent program under the correspondence g 7→ g′ and keeps the other variables fixed. Moreover, if all the G-vars in Hℓ−1 are non-degenerate in the infinite-width execution and non-linearities are Lipschitz continuous, then (iii) yields an equivalent program under the same correspondence. Remark 3.12. The equivalence (iii) changes the filtration of the Netsor program by substituting the G-var g with the H-var g′ . Informally, this provides a way to eliminate degenerate conditional infinite-width executions of G-vars by skipping “redundant” levels of the filtration and replacing the G-var with a linear combination of previous G-vars, which can be implemented as a NonLin line with a linear map. Proof of Lemma 3.9. We construct a finite sequence of programs T0 , T1 , . . . , TN , starting with T0 = T , where each program is equivalent to the previous one and has strictly fewer degenerate MatMul lines. Iterate through the lines of the current program T sequentially, from ℓ = 1 to M . Suppose line ℓ is the first MatMul instruction g(ℓ) = Wh
24
ANDREA AGAZZI, ELOY MOSIG GARCÍA, AND DARIO TREVISAN
that is degenerate in the sense of Proposition 3.11, i.e., ḡ (ℓ) | F ℓ−1 is a degenerate random variable. Because we proceed sequentially, all G-vars generated by MatMul lines prior to ℓ are already non-degenerate in the infinite-width execution conditionally upon F. We apply Proposition 3.11 to replace the instruction at line ℓ with the NonLin instruction P g′ := ki=1 αi g(i) , yielding a new program Tnew . By Proposition 3.11.(iii), Tnew is infinite-width-equivalent to the previous program. Since the nonlinearities are assumed to be continuous, the final step of Proposition 3.11 applies and the replacement yields an equivalent program. We update our current program to Tnew and continue the iteration. Since the original program T has a finite number of lines M , and each substitution strictly removes a single degenerate MatMul line without introducing new ones, this procedure terminates. The final program TN := T ′ is equivalent to T by transitivity, and every remaining MatMul line has strictly positive conditional variance. □ Proof of Proposition 3.11. (i)⇔(ii). By Equation (3.26) Σḡ|F ℓ−1 = LW|F ℓ−1 h̄∗ RW|F ℓ−1 h̄ , so conditional degeneracy of ḡ given F ℓ−1 is equivalent to vanishing of the scalar factor h̄∗ RW|F ℓ−1 h̄. (ii)⇒(iii). By Lemma 3.6, the identity h̄∗ RW|F ℓ−1 h̄ = 0 implies that h̄ belongs to the span of the previous parents {h̄i : g(i) = Wh(i) occurs before line ℓ}. Since only finitely many parents occur before ℓ, there exist deterministic coefficients (αi )ki=1 such that (3.28)
h̄ =
k X
αi h̄i
almost surely.
i=1
Using the bilinearity of the covariance, we can expand the covariance matrix of the P difference ḡ − ki=1 αi ḡ (i) . Applying (1.3) to each term yields: Σḡ−Pk
(i) i=1 αi ḡ
= Σḡ − 2
k X
αi Σḡ,ḡ(i) +
i=1
= (h̄∗ h̄)Idn − 2
k X
αi αj Σḡ(i) ,ḡ(j)
i,j=1 k X
αi (h̄∗ h̄i )Idn +
i=1
k X
αi αj (h̄i ∗ h̄j )Idn .
i,j=1
By the bilinearity of the infinite-width pairing ∗ , this expression factors into: ! k k X X Σ Pk h̄ − αi h̄i ∗ h̄ − αj h̄j Idn . (i) = ḡ−
i=1 αi ḡ
i=1
j=1
By (3.28), the term inside the brackets is zero, meaning this covariance matrix is identically zero. Since a centered Gaussian vector with zero variance must be zero almost surely, we Pk (i) almost surely, hence replacing the G-var g = Wh with the Hobtain ḡ = i=1 αi ḡ P var g′ = i αi g(i) yields an infinite-width-equivalent program. To show the last claim, we introduce continuous functions (ϕ(i) )i=1,...,k and ϕ such that ϕ(i) ((ḡ (j) )j=1,...,k′ ) = h̄(i) and similarly ϕ such that ϕ((ḡ (j) )j=1,...,k′ ) = h̄, where possibly the set of G-vars involved (g(j) )j=1,...,k′ is different than (ḡ (i) )i=1,...,k , but still contained in Hℓ−1 , hence they are ′ non-degenerate Gaussians. Consider the continuous function Φ : Rk → R given by Φ((zj )j=1,...,k′ ) := ϕ((zi )j=1,...,k′ ) −
k X i=1
αi ϕ(i) ((zj )j=1,...,k′ ).
QUANTITATIVE GAUSSIAN-PROCESS LIMITS OF TENSOR PROGRAMS
25
Then Φ((ḡ (j) )j=1,...,k′ ) = 0 (componentwise in Rn ) Lebesgue-almost surely. Since the ′ (ḡ (j) )j=1,...,k′ are non-degenerate, the set {z ∈ Rk : Φ(z) ̸= 0} has null Lebesgue measure, and by continuity of Φ it is an open set, hence empty. In particular, Φ((g (j) )j=1,...,k′ ) = 0 almost surely under finite-width execution as well, so the same linear relation (3.28) holds in finite-width execution. Multiplying by W gives g ′ = g almost surely. P (i) (iii)⇒(ii). Assume ḡ = i αi ḡ where each ḡ (i) = Wh̄(i) uses the same W and appears before ℓ. By Lemma 3.6, the conditional covariance of the stacked block (ḡ (1) , . . . , ḡ (k) , ḡ) equals a Kronecker product of In with the Gram matrix of the parents (h̄(1) , . . . , h̄(k) , h̄). The linear relation among the ḡ’s forces this Gram matrix to be singular, hence h̄ lies in the span of previous parents. Equivalently, its orthogonal projection coefficient vanishes: □ h̄∗ RW|F ℓ−1 h̄ = 0. 4. Proof of main results The proof of Theorem 1.2 is by induction on the program lines. A key input in the induction is a quantitative kernel LLN for pairs of H-vars already constructed before the current line. To avoid any circularity, we first prove this kernel estimate in the following conditional form: if the conclusion of Theorem 1.2 holds for the truncated program up to line r − 1, then the corresponding kernel LLN holds for all H-vars generated up to line r − 1. This conditional kernel estimate will be used in the induction step from r − 1 to r. Once the induction proving Theorem 1.2 is complete, the same argument gives Corollary 1.3 for the full program. We therefore begin with the proof of the kernel estimate, and then turn to the proof of Theorem 1.2. For the sake of completeness we recall some elementary inequalities involving the Wasserstein distance that will be used in the proofs of Theorem 1.2 and Corollary 1.3. The reader interested in a general treatment of the Wasserstein distance is referred to [Vil08]. Given random variables (X, X ′ ), (Y, Y ′ ) taking values in Rn × Rn , it holds (4.1) Wp X ⊤ X ′ , Y ⊤ Y ′ ≤ X ′ L2p + ∥Y ∥L2p W2p (X, X ′ ), (Y, Y ′ ) , with c = c(p) < ∞. Indeed, consider any coupling of (X, X ′ ) and (Y, Y ′ ). By adding and subtracting Y ⊤ X ′ , (4.2) X ⊤ X ′ − Y ⊤ Y ′ ≤ (X − Y )⊤ X ′ + Y ⊤ (X ′ − Y ′ ) ≤ ∥X − Y ∥ ∥X ′ ∥ + ∥Y ∥ ∥X ′ − Y ′ ∥. Taking the Lp norm, using the triangle inequality and Cauchy-Schwarz inequality, we obtain (4.1). Similarly, the following bound holds: (4.3) Wp X ⊤ X ′ , Y ⊤ Y ′ ≤ 1 + (Y, Y ′ ) L2p + W2p (X ′ , Y ′ ) W2p (X, X ′ ), (Y, Y ′ ) . This is obtained from using the triangle inequality ∥X ′ ∥ ≤ ∥X ′ −Y ′ ∥+∥Y ′ ∥ in (4.2) above. 4.1. Proof of Corollary 1.3. Throughout this proof, c < ∞ denotes a constant that may depend on p and on the structural constants of the program, but not on the hidden widths. Its value may change from line to line. Moreover, we can assume p ≥ 2 without loss of generality, and also that n ∈ / nin , otherwise the thesis is trivial because both variables are deterministic, noticing that the left hand side is zero in this case. Let h, h′ : H(n) be two H-vars of common length n, with finite-width executions (h, h′ ) and infinite-width executions (h̄, h̄′ ). Set 1 1 K̄n := h̄⊤ h̄′ , K := E[K̄n ] = h̄∗ h̄′ . Kn := h⊤ h′ , n n Since K is deterministic, ∥Kn − K∥Lp = Wp (Kn , K).
26
ANDREA AGAZZI, ELOY MOSIG GARCÍA, AND DARIO TREVISAN
By the triangle inequality for Wasserstein distances, (4.4)
∥Kn − K∥Lp = Wp (Kn , K) ≤ Wp (Kn , K̄n ) + Wp (K̄n , K).
We estimate the two terms separately. For the first term, apply (4.3) with h X := √ , n
h′ X ′ := √ , n
h̄′ Y ′ := √ . n
h̄ Y := √ , n
Then X ⊤ X ′ = Kn and Y ⊤ Y ′ = K̄n , and therefore Wp (Kn , K̄n ) ≤ cp 1 + ∥(Y, Y ′ )∥L2p + W2p (X ′ , Y ′ ) W2p (X, X ′ ), (Y, Y ′ ) . The moment factor ∥(Y, Y ′ )∥L2p is bounded uniformly in the widths by the infinite-width moment bounds. Moreover, ! ! (ℓ) (ℓ) h h̄ . , √ W2p (X ′ , Y ′ ) ≤ W2p (X, X ′ ), (Y, Y ′ ) ≤ W2p √ nℓ nℓ ℓ=1,...,M
ℓ=1,...,M
Using (1.5) with exponent 2p, we obtain therefore: X
W2p (X ′ , Y ′ ) + W2p (X, X ′ ), (Y, Y ′ ) ≤ c
m∈nhidd
1 √ . m
Thus (4.5)
Wp (Kn , K̄n ) ≤ c
X m∈nhidd
1 √ . m
It remains to estimate the sampling term Wp (K̄n , K). Since K is deterministic, Wp (K̄n , K) = ∥K̄n − K∥Lp . We write
n
K̄n =
1X Ui , n
Ui := h̄i h̄′i .
i=1
By definition of the infinite-width execution, the variables (Ui )ni=1 are i.i.d., and K = E[U1 ]. Let Ũi := Ui − E[Ui ] for i = 1, . . . , n. Then Rosenthal’s inequality (3.6) gives, for each p ≥ 2, 1/p !p/2 n n n p X X ∥Ũi ∥2L2 ∥ Ũ ∥ 1X i Lp Ũi p ≤ c(p) + n n2 np L i=1 i=1 i=1 ! ∥Ũ1 ∥L2 ∥Ũ1 ∥Lp √ ≤ c(p) + 1−1/p , n n where the constant c(p) depends on p only. Moreover, since U1 = h̄1 h̄′1 has finite moments by Lemma 3.2 and 1 − 1/p ≥ 1/2 for p ≥ 2, we can further simplify: n
1X c Ũi p ≤ √ , n n L i=1
where c depends on the p-th moment of the Ũi , hence it is a structural constant of the program. Therefore, (4.4) can be estimated with: X X 1 c 1 √ +√ ≤c √ . ∥Kn − K∥Lp ≤ c m n m m∈n m∈n ∪n hidd
hidd
out
QUANTITATIVE GAUSSIAN-PROCESS LIMITS OF TENSOR PROGRAMS
27
4.2. Proof of Theorem 1.2. We prove the theorem by constructing an explicit coupling between the finite-width and infinite-width executions, line by line. Throughout the proof, constants denoted by c may depend on p and on the structural constants of the program, but never on the hidden widths. By Lemma 3.9, we may assume without loss of generality that the program is nondegenerate in the infinite-width execution. In particular, the Gram matrices appearing in the structure formulas for the conditional laws stay in a fixed smooth stratum at the infinite-width point, so that the maps involving pseudoinverses and matrix square roots are differentiable there. We denote by L the number of lines in the program. For r ≤ L, define the normalized error vector ! h(ℓ) − h̄(ℓ) . ∆r := √ nℓ ℓ=1,...,r
We prove by induction on r that the coupling can be chosen so that X 1 √ , (4.6) ∥∆r ∥Lp ≤ c m m∈nhidd (r)
where nhidd (r) denotes the hidden widths of the program restricted to the first r lines. The theorem follows immediately from (4.6), since Wasserstein distance is bounded above by the Lp -distance under any coupling. The case r = 0 is trivial (T contains only input, hence deterministic, H-vars). Assume that the coupling has been constructed up to line r − 1 and that (4.6) holds at level r − 1. We distinguish the two possible operations at line r. The NonLin case. Assume that line r is h(r) = ϕ(h1 , . . . , hk ), where h1 , . . . , hk : H(nr ) are previously defined H-vars. We extend the coupling by setting h(r) = ϕ(h1 , . . . , hk ),
h̄(r) = ϕ(h̄1 , . . . , h̄k ).
Since ϕ is Lipschitz and acts componentwise, we have pointwise 1/2 k X ∥h(r) − h̄(r) ∥ ≤ Lip(ϕ) ∥hj − h̄j ∥2 . j=1
√ Dividing by nr gives k X hj − h̄j h(r) − h̄(r) ≤ Lip(ϕ) √ √ nr nr
2
1/2
.
j=1
Therefore, h(r) − h̄(r) √ nr
≤ Lip(ϕ) Lp
hj − h̄j √ nr
. j=1,...,k Lp
Since h1 , . . . , hk are previous H-vars of the same shape nr , the last term is bounded by the induction hypothesis: X h(r) − h̄(r) 1 √ . ≤c √ nr m p L
m∈nhidd (r−1)
Combining this estimate with the induction hypothesis for the previous coordinates yields X 1 √ , ∥∆r ∥Lp ≤ c m m∈nhidd (r)
noticing that nhidd (r − 1) = nhidd (r).
28
ANDREA AGAZZI, ELOY MOSIG GARCÍA, AND DARIO TREVISAN
The MatMul case. Assume that line r is g = Wh,
h : H(mr ),
W : A(nr , mr ),
where nr and mr do not necessarily belong to nhidd (r −1), because nr may be a completely new width, and mr could be an input dimension or an output dimension for the restricted program up to line r − 1. In the former case (as well as in the case that mr ∈ nhidd (r − 1) already) it holds nhidd (r) = nhidd (r − 1), while in the latter we have nhidd (r) = nhidd (r − 1)∪{mr }. Let g and ḡ be the finite-width and infinite-width executions of g. Conditionally on the past lines, both are Gaussian: g | Fr−1 ∼ N µg|Fr−1 , Σg|Fr−1 , and
ḡ | F r−1 ∼ N µḡ|F r−1 , Σḡ|F r−1 .
The conditional means and covariances are given by the structure lemmas for the finitewidth and infinite-width executions. We extend the coupling by taking a standard Gaussian vector Z ∼ N (0, Idnr ), independent of the past, and setting q q (4.7) h(r) = g := µg|Fr−1 + Σg|Fr−1 Z, h̄(r) = ḡ := µḡ|F r−1 + Σḡ|F r−1 Z. This gives the correct conditional laws, hence the correct marginal laws up to line r. We now prove the two stability estimates X µg|Fr−1 − µḡ|F r−1 1 √ , (4.8) ≤c √ nr m Lp m∈nhidd (r) q q X √ 1 √ . (4.9) Σg|Fr−1 − Σḡ|F r−1 ≤ c nr p m L m∈nhidd (r)
Assuming (4.8)-(4.9) for the moment, we conclude the MatMul step. Let q q Ar := Σg|Fr−1 − Σḡ|F r−1 . By (4.7), ∥µg|Fr−1 − µḡ|F r−1 ∥ ∥Ar Z∥ ∥g − ḡ∥ ≤ + √ . √ √ nr nr nr Since Z is standard Gaussian and independent of the past, the Gaussian moment bound gives ∥Ar Z∥Lp ≤ cp ∥Ar ∥Lp . Therefore, by (4.8) and (4.9), X g − ḡ 1 √ . ≤c √ nr Lp m m∈nhidd (r)
Combining this with the induction hypothesis gives ∥∆r ∥Lp ≤ ∥∆r−1 ∥Lp +
g − ḡ √ nr
≤c Lp
X m∈nhidd (r)
1 √ . m
It remains to prove the two stability estimates. We first dispose of the case mr ∈ nin . In this case the parent h : H(mr ) is an input variable, or a NonLin transformation of them, hence its finite-width and infinite-width executions coincide deterministically. Consequently the finite-width and infinite-width Gram data entering the structure formulas are identical. The conditional means and covariances of g and ḡ therefore coincide, and the left-hand sides of both stability estimates vanish. Hence, we may assume from now on that mr ∈ / nin . We prove separately the case of mean and covariance.
QUANTITATIVE GAUSSIAN-PROCESS LIMITS OF TENSOR PROGRAMS
29
It remains to prove the two stability estimates. We first consider the case mr ∈ nin . In the basic Netsor language, the current parent h : H(mr ) is obtained from deterministic input variables by ordinary NonLin operations. Hence its finite-width and infinite-width executions coincide deterministically, and no new kernel LLN over the input dimension mr is required. The Gram entries involving the current parent are structural. However, if the same A-var W has already been used earlier in the program, the conditional means and covariances need not coincide exactly: the only possible discrepancy comes from previous uses of W, and is already controlled by the induction hypothesis and by the kernel estimates for the previous hidden parents. Therefore the estimates below remain valid with nhidd (r) = nhidd (r − 1). We may then continue with the same mean and covariance stability arguments, noting −1/2 that no additional mr contribution is generated when mr ∈ nin . Covariance stability. By the structure lemmas, there exist scalar quantities ar , ār ≥ 0 such that Σḡ|F r−1 = ār Idnr . Σg|Fr−1 = ar Idnr , Moreover, ar and ār are obtained by applying the same finite-dimensional map to the finite-width and infinite-width Gram data associated with the previous parents of the weight W. More precisely, 1 (r−1),⊤ (r−1) (r−1) (r−1) ar = f , ār = f H W ∗ H W HW HW , mr for a map f depending only on the program. In particular, by the structure lemmas, f (X) is exactly the Schur complement X/Y , where Y is the square submatrix of X obtained (r−1),⊤ (r−1) 1 HW HW , the by removing the last row and column. In our case, when X = m (r−2),⊤
1 submatrix Y is given by m HW
(r−2)
HW
(r−1) ∗
and when X = H W
(r−1)
HW
, Y is given by
(r−2) (r−2) H W ∗H W .
p Since the program is non-degenerate in infinite-width execution, the map x 7→ f (x) is differentiable at the infinite-width Gram point. Lemma 3.1, applied to this map, together with the kernel LLN for the truncated program up to line r − 1, yields X X √ √ 1 1 √ =c √ , ∥ ar − ār ∥Lp ≤ c m m m∈nhidd (r−1)∪{mr }
m∈nhidd (r)
since mr is either a hidden width or an output dimension for the restricted program. Moreover, q q √ √ Σg|Fr−1 − Σḡ|F r−1 = ( ar − ār ) Idnr , we obtain q q Σg|Fr−1 − Σḡ|F r−1
Lp
=
√
√ √ √ nr ∥ ar − ār ∥Lp ≤ c nr
X m∈nhidd (r)
1 √ , m
which proves (4.9). Mean stability. We prove a slightly stronger estimate. Fix the parent h : H(m) appearing in line r, and write h, h̄ for its finite-width and infinite-width executions. For every s ≤ r − 1, we claim that X µW |Fs h − µ̄W |F s h̄ 1 √ . (4.10) ≤c √ nr m p L
m∈nhidd (r)
Taking s = r − 1 gives (4.8). We prove (4.10) by induction over s. For s = 0, both conditional means are zero. Assume the estimate holds at level s − 1. If line s is not a MatMul line using the weight
30
ANDREA AGAZZI, ELOY MOSIG GARCÍA, AND DARIO TREVISAN
W, then the conditional mean of W does not change and there is nothing to prove. Thus suppose line s is g′ = Wh′ , with executions (g ′ , h′ ) and (ḡ ′ , h̄′ ). The update formula for the conditional mean gives µW |Fs h = µW |Fs−1 h + Zs us , where ′
′
Zs := g − µW |Fs−1 h
q
(h′⊤ RW |Fs−1 h′
+
∈ Rnr ,
and us :=
q
h′⊤ RW |Fs−1 h′
+ ′⊤ h RW |Fs−1 h ∈ R.
Similarly, µ̄W |F s h̄ = µ̄W |F s−1 h̄ + Z̄s ūs , with the analogous barred definitions. By the coupling construction at line s, the standardized residuals are driven by the same Gaussian noise, hence Zs = Z̄s under our coupling. Moreover, 1 Idn . Zs | Fs−1 ∼ N 0, mr r Therefore, µW |Fs h − µ̄W |F s h̄ √ nr
≤ Lp
µW |Fs−1 h − µ̄W |F s−1 h̄ √ nr
1 + √ ∥Zs (us − ūs )∥Lp . nr p
L
By Hölder’s inequality and Gaussian moment estimates, c 1 1 √ ∥Zs (us − ūs )∥Lp ≤ √ ∥Zs ∥L2p ∥us − ūs ∥L2p ≤ √ ∥us − ūs ∥L2p . nr nr mr It remains to control us − ūs . By the projector representation of the right covariance, both us and ūs are obtained from the corresponding Gram matrices via the same finitedimensional map: 1 (s−1),⊤ (s−1) (s−1) (s−1) us = fs HW HW , ūs = fs H W ∗ H W . mr By non-degeneracy, fs is differentiable at the infinite-width Gram point. Moreover, q + 1/2 |us | = h′⊤ RW |Fs−1 h′ h′⊤ RW |Fs−1 h ≤ ∥RW |Fs−1 h∥ ≤ ∥h∥, since ∥RW |Fs−1 ∥op ≤ 1, and the same bound holds for ūs . Thus the uniform moment bounds imply √ ∥us ∥L4p + ∥ūs ∥L4p ≤ c mr . Applying Lemma 3.1 to fs , together with the kernel LLN for the truncated program, yields X X √ √ 1 1 √ = c mr √ . ∥us − ūs ∥L2p ≤ c mr m m m∈nhidd (r−1)∪{mr }
m∈nhidd (r)
Consequently, 1 √ ∥Zs (us − ūs )∥Lp ≤ c nr
X m∈nhidd (r)
Iterating over s ≤ r − 1 proves (4.10). Hence (4.8) follows.
1 √ . m
QUANTITATIVE GAUSSIAN-PROCESS LIMITS OF TENSOR PROGRAMS
31
The NonLin and MatMul cases close the induction. Taking r = L in (4.6) gives h(ℓ) − h̄(ℓ) √ nℓ
! X
≤c
m∈nhidd
ℓ=1,...,L Lp
1 √ . m
Therefore,
h(ℓ) Wp √ nℓ
!L , ℓ=1
h̄(ℓ) √ nℓ
!L ≤c ℓ=1
X m∈nhidd
1 √ , m
as claimed. 4.3. Proof of Theorem 2.4. The proof follows the line-by-line coupling used for Theorem 1.2, enlarged to include the scalar variables. The MatMul case is governed by the same conditional Gaussian regression formulas as in the basic Netsor proof. Indeed, scalar variables are measurable functions of the previously generated variables and do not introduce additional observations of the weights. Thus, conditionally on the previous lines, the finite-width execution g of a MatMul line, g = Wh, is still Gaussian with conditional mean and covariance given by the structure formulas of Section 3.2. There is, however, one minor difference with respect to the basic Netsor setting. In a NetsorK program, an H-var with structural input dimension m ∈ nin need not be one of the deterministic input variables: it may have been obtained from input variables through a NonLin+ instruction depending on scalar variables. Such a variable is therefore not necessarily deterministic in the finite-width execution, as opposed to the variables obtained from NonLin instructions applied to input variables. Its infinite-width execution is in any case deterministic, and its finite-width fluctuation is already controlled by the induction hypothesis on the preceding NonLin+ and Scal lines. Consequently, no new kernel LLN over an input dimension is required. For variables of hidden dimension m ∈ nhidd , the argument is identical to the Netsor case. The local Lipschitz assumption on NetsorK programs ensures that the finitedimensional maps appearing in the structure formulas are differentiable at the infinitewidth Gram point. Assuming that the r-th line is a Kernel line s = h⋆ h′ , with h, h′ : H(mr ), the two executions are 1 ⊤ ′ h h, s̄ = h̄∗ h̄′ . s= mr Thus the estimate ∥s − s̄∥Lp ≤ c
X m∈nhidd (r)
1 √ m
is exactly Corollary 1.3, applied to the truncated program up to the current line. Take into account that mr belongs to nhidd (r), as opposed to the case of Netsor. For a Scal line, s = ψ(s1 , . . . , su ), the estimate follows from the finite-dimensional local Lipschitz property of ψ, the induction hypothesis for the previous scalar variables, and the moment bounds for scalar variables, which are obtained inductively from the kernel LLN and the growth assumptions. It remains to discuss a NonLin+ line, h = ϕ(h1 , . . . , hk ; s1 , . . . , su ) : H(mr ).
32
ANDREA AGAZZI, ELOY MOSIG GARCÍA, AND DARIO TREVISAN
Let h, h̄ be the corresponding executions. Applying (2.5) componentwise gives |ha − h̄a | ≤ C(1 + ∥s∥ + ∥s̄∥)
k X
!1/2 2
|hi,a − h̄i,a |
i=1
+ C 1 +
k X
!1/2 |hi,a |2
+
i=1
k X
!1/2 |h̄i,a |2
|s − s̄|,
i=1
where s = (s1 , . . . , su ) and s̄ = (s̄1 , . . . , s̄u ). Hence, squaring and summing over a = 1, . . . , mr , we obtain eventually ! k k X X . ∥h − h̄∥ ≤ c(1 + ∥s∥ + ∥s̄∥) hi − h̄i + c ∥s − s̄∥ 1 + ∥hi ∥ + h̄i i=1
i=1
√
Dividing both sides by mr and taking Lp norms and using Hölder’s inequality, the induction hypothesis controls the differences, while the normalized moment bounds for the previously constructed H-vars and C-vars control the coefficients. This yields X h − h̄ √ ≤c mr Lp m∈n
1 √ , m
hidd
and induction closes as in the proof of Theorem 1.2.
5. Numerical Experiments In this section, we validate the quantitative wide-limit convergence of finite-width random neural networks toward their infinite-width limits. We focus on the convergence of the network’s final output law via an ensemble sampling paradigm. For fully connected architectures, this output-law experiment is designed to reveal the sharper O(n−1 ) convergence rate known under non-degeneracy hypotheses [Tre23; Fav+25]. For the other architectures, the experiments should be interpreted as evidence of the finite-width convergence predicted by the present theorem, rather than as a theorem-level claim of an O(n−1 ) rate. 5.1. Experimental Setup and Architectures. We evaluate a suite of architectures of increasing complexity, corresponding to the examples introduced in Section 2.1, to illustrate the finite-width convergence predicted by Theorem 1.2 across several architectures. In all of the following experiments, the infinite-width limit moments are computed using deterministic Gauss-Hermite quadrature to avoid integration noise. Multi-Layer Perceptrons: We evaluate both shallow and deep Multi-Layer Perceptrons using tanh and ReLU activations, mirroring the structures of Algorithm 1 and Algorithm 2, respectively. These setups capture standard layer-wise transitions where separate, independent weight matrices are sampled per layer. The deep MLP experiment (L = 4) demonstrates how compounding non-linearities propagate through the infinite-width recursion while maintaining the finite-width convergence limits. Time-unrolled Recurrent Neural Networks: We test recurrent networks following the tiedweight architecture of Algorithm 3. In these experiments, the network is unrolled over T = 5 time steps with a tanh activation to demonstrate convergence under weight-sharing constraints. Because the recurrent matrix Wh is shared across time, the infinite-width computation must track cross-time second moments Kt,s = E[ht h⊤ s ]. The limiting marginal laws at any given time step t are then extracted from the diagonal blocks Kt,t .
QUANTITATIVE GAUSSIAN-PROCESS LIMITS OF TENSOR PROGRAMS
33
Residual Networks: To further validate the framework, we evaluate deep networks with skip connections, formulated as h(ℓ) = h(ℓ−1) + ϕ(W(ℓ) h(ℓ−1) ) as detailed in Algorithm 6, where ϕ is the ReLU activation. Because a new, independent weight matrix W(ℓ) is drawn on each residual update, the samples for the current state and its update can be treated as independent when constructing marginal samples for the subsequent layer. Netsor program 6 Residual Network Input: x : H(d) Input: W(0) : A(n, d) Input: W(ℓ) : A(n, n) for ℓ = 1, . . . , L 1 g(0) = W(0) x : G(n) 2 h(0) = g(0) : H(n) for ℓ = 1, . . . , L do 3 g(ℓ) = W(ℓ) h(ℓ−1) : G(n) 4 h(ℓ) = h(ℓ−1) + ϕ(g(ℓ) ) : H(n) end for
▷ MatMul ▷ NonLin ▷ MatMul ▷ NonLin
5.2. Output-Law Convergence and Estimation Noise. To quantify the convergence of the output distributions, we utilize the Sliced Wasserstein-1 distance. For probability measures µ and ν on RM , it is defined as: Z (5.1) SW 1 (µ, ν) = W1 (θ♯ µ, θ♯ ν) dσ(θ), SM −1
where θ♯ denotes the pushforward measure onto the direction θ ∈ SM −1 , and σ is the uniform measure on the unit sphere. The bound SW 1 ≤ W1 holds inherently because the orthogonal projection onto θ is a 1-Lipschitz map. Moreover, SW 1 circumvents the computational complexity of W1 associated with the curse of dimensionality, making it preferred for empirical evaluations. For each architectural width n, we initialize an ensemble of N = 5000 independent networks and extract the scalar readouts. By evaluating these readouts across independent seeds, we decouple the empirical sample size from the network width n. The SW 1 distance is computed between these N outputs and the analytical distribution of the activations of the limiting Gaussian process. For an empirical measure µ̂N constructed from N i.i.d. samples of a distribution µ, Theorem 1 in [FG15] yields E[SW 1 (µ, µ̂N )] ≤ c1 N −1/2 . Let µn denote the exact finitewidth distribution of the network activations, µ̂n,N denote the empirical measure formed by the N ensemble samples, and µ∞ denote the limiting distribution of the Gaussian process activations. By the triangle inequality, the empirical distance we observe is bounded by the true architectural convergence plus this statistical noise: c2 c1 c1 +√ , (5.2) SW 1 (µ̂n,N , µ∞ ) ≤ SW 1 (µn , µ∞ ) + √ ≤ n N N where we substitute the expected theoretical rate O(n−1 ) from [Fav+25; Tre23] for MLP architectures. The sampling error N −1/2 does √ not dominate the architectural signal for network widths up to approximately n ≤ N ≈ 70; for larger n, the observed slope is expected to flatten toward −1/2 as the estimation noise becomes the leading term. This behavior is consistent with the fitted slopes reported in Table 1, which lie strictly between −1/2 and −1 across all architectures, reflecting a mixture of the true convergence signal and the sampling floor. In Figure 7 we plot the Sliced Wasserstein results for our four experiments in log-log scale.
34
ANDREA AGAZZI, ELOY MOSIG GARCÍA, AND DARIO TREVISAN
For each architecture, we estimate the empirical convergence exponent by a least-squares fit of log dn = α + β log n, where dn denotes the SW 1 distance between the finite-width execution at width n and the infinite-width execution of each selected architecture. The fitted value of β is reported in Table 1. Slopes strictly below −1/2 indicate that the output-law convergence is visible beyond the sampling noise floor.
(a) Shallow MLP
(b) Deep MLP (L = 4)
(c) Residual Network
(d) Time-unrolled RNN
Fig. 7. Sliced Wasserstein-1 (SW 1 ) distances between the empirical output laws (N = 5000 ensemble samples) and the law of the limiting Gaussian process activations across four architectures, plotted in log-log scale. Architecture SW 1 slope Shallow MLP −0.721 Deep MLP L = 4 −0.916 Residual network −0.797 Time-unrolled RNN −0.579 Table 1. Least-squares log-log slopes for the SW 1 curves shown in Figure 7. All slopes lie strictly below −1/2, confirming that the architectural convergence signal is visible above the empirical sampling floor.
REFERENCES
35
References [Aga+26]
A. Agazzi, G. Bruno, E. M. Garcı́a, S. Saviozzi, and M. Romito. Stochastic Scaling Limits and Synchronization by Noise in Deep Transformer Models. 2026. arXiv: 2604.26898 [math.PR] (cit. on p. 5). [And05] T. Ando. “Schur Complements and Matrix Inequalities: Operator-Theoretic Approach”. In: The Schur Complement and Its Applications. Ed. by F. Zhang. Boston, MA: Springer US, 2005, pp. 137–162 (cit. on p. 14). [Aro+19] S. Arora, S. S. Du, W. Hu, Z. Li, R. Salakhutdinov, and R. Wang. “On exact computation with an infinitely wide neural net”. In: Proceedings of the 33rd International Conference on Neural Information Processing Systems. Red Hook, NY, USA: Curran Associates Inc., 2019 (cit. on p. 5). [BFF25] A. Bordino, S. Favaro, and S. Fortini. Non-asymptotic approximations of Gaussian neural networks via second-order Poincaré inequalities. 2025. arXiv: 2304.04010 [cs.LG] (cit. on p. 5). [BT22] A. Basteri and D. Trevisan. Quantitative Gaussian Approximation of Randomly Initialized Deep Neural Networks. cvgmt preprint. 2022 (cit. on pp. 1, 5–7). [Cam+23] V. Cammarota, D. Marinucci, M. Salvi, and S. Vigogna. A Quantitative Functional Central Limit Theorem for Shallow Neural Networks. 2023. arXiv: 2306. 16932 [math.PR] (cit. on p. 5). [Cel26] L. Celli. “Wide neural networks with general weights: convergence rate and explicit dependence on the hyper-parameters”. In: arXiv preprint arXiv:2601.21539 (2026) (cit. on pp. 1, 6). [EMS21] R. Eldan, D. Mikulincer, and T. Schramm. “Non-asymptotic approximations of neural networks by Gaussian processes”. In: Proceedings of Thirty Fourth Conference on Learning Theory. Ed. by M. Belkin and S. Kpotufe. Vol. 134. Proceedings of Machine Learning Research. PMLR, 2021, pp. 1754–1775 (cit. on p. 5). [Fav+25] S. Favaro, B. Hanin, D. Marinucci, I. Nourdin, and G. Peccati. “Quantitative CLTs in deep neural networks”. In: Probability Theory and Related Fields 191.3 (2025), pp. 933–977 (cit. on pp. 1, 5, 32, 33). [FG15] N. Fournier and A. Guillin. “On the rate of convergence in Wasserstein distance of the empirical measure”. In: Probability Theory and Related Fields 162.3 (2015), pp. 707–738 (cit. on p. 33). [GDP25] F. Girardi and G. De Palma. “Trained quantum neural networks are gaussian processes”. In: Communications in Mathematical Physics 406.4 (2025), p. 92 (cit. on p. 6). [GKR26] F. Giovagnini, S. Kotitsas, and M. Romito. “Universality in Deep Neural Networks: An approach via the Lindeberg exchange principle”. In: arxiv preprintarxiv:2605.02771 (2026) (cit. on pp. 1, 6). [GM+18] A. G. de G. Matthews, M. Rowland, J. Hron, R. E. Turner, and Z. Ghahramani. Gaussian Process Behaviour in Wide Deep Neural Networks. 2018. arXiv: 1804.11271 [stat.ML] (cit. on p. 5). [Gre60] T. N. E. Greville. “Some Applications of the Pseudoinverse of a Matrix”. In: SIAM Review 2.1 (1960), pp. 15–22 (cit. on p. 14). [Hay24] S. Hayou. “Commutative Scaling of Width and Depth in Deep Neural Networks”. In: Journal of Machine Learning Research 25.299 (2024), pp. 1–41 (cit. on p. 5). [Hro+20] J. Hron, Y. Bahri, J. Sohl-Dickstein, and R. Novak. Infinite attention: NNGP and NTK for deep attention networks. 2020. arXiv: 2006.10540 [stat.ML] (cit. on p. 6).
36
[JGH20]
[KGR26] [Lee+20]
[MAT26]
[MH+25]
[Nea96] [PF20]
[Rao73] [Ros70]
[SKI26]
[SM50]
[Tre23]
[UK97]
[Vil08] [Yan20a]
[Yan20b] [Yan21]
[YH20]
REFERENCES
A. Jacot, F. Gabriel, and C. Hongler. Neural Tangent Kernel: Convergence and Generalization in Neural Networks. 2020. arXiv: 1806 . 07572 [cs.LG] (cit. on pp. 5, 6). H. Koubbi, B. Geshkovski, and P. Rigollet. Homogenized Transformers. 2026. arXiv: 2604.01978 [math.PR] (cit. on p. 5). J. Lee, L. Xiao, S. S. Schoenholz, Y. Bahri, R. Novak, J. Sohl-Dickstein, and J. Pennington. “Wide neural networks of any depth evolve as linear models under gradient descent *”. In: Journal of Statistical Mechanics: Theory and Experiment 2020.12 (Dec. 2020), p. 124002 (cit. on pp. 1, 6). E. Mosig, A. Agazzi, and D. Trevisan. Quantitative convergence of trained single layer neural networks to Gaussian processes. 2026. arXiv: 2509.24544 [stat.ML] (cit. on p. 6). A. Melchor Hernandez, F. Girardi, D. Pastorello, and G. De Palma. “Quantitative Convergence of Trained Quantum Neural Networks to a Gaussian Process: A. Melchor Hernandez et al.” In: Annales Henri Poincaré. Springer. 2025, pp. 1–57 (cit. on p. 6). R. M. Neal. Bayesian Learning for Neural Networks, Vol. 118 of Lecture Notes in Statistics. Springer-Verlag, 1996 (cit. on pp. 1, 5). S. Peluchetti and S. Favaro. “Infinitely deep neural networks as diffusion processes”. In: Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics. Ed. by S. Chiappa and R. Calandra. Vol. 108. Proceedings of Machine Learning Research. PMLR, 2020, pp. 1126–1136 (cit. on p. 5). C. R. Rao. Linear Statistical Inference and its Applications. 2nd. John Wiley & Sons, 1973 (cit. on p. 16). H. P. Rosenthal. “On the subspaces of Lp (p > 2) spanned by sequences of independent random variables”. In: Israel Journal of Mathematics 8.3 (1970), pp. 273–303 (cit. on p. 15). M. Sakai, R. Karakida, and M. Imaizumi. Infinite-Width Limit of a Single Attention Layer: Analysis via Tensor Programs. 2026. arXiv: 2506 . 00846 [cs.LG] (cit. on pp. 6, 13). J. Sherman and W. J. Morrison. “Adjustment of an Inverse Matrix Corresponding to a Change in One Element of a Given Matrix”. In: The Annals of Mathematical Statistics 21.1 (1950), pp. 124 –127 (cit. on p. 14). D. Trevisan. Wide Deep Neural Networks with Gaussian Weights are Very Close to Gaussian Processes. 2023. arXiv: 2312 . 11737 [math.ST] (cit. on pp. 1, 5, 7, 32, 33). F. E. Udwadia and R. E. Kalaba. “An Alternative Proof of the Greville Formula”. In: Journal of Optimization Theory and Applications 94.1 (1997), pp. 23–28 (cit. on p. 14). C. Villani. “Optimal transport – Old and new”. In: vol. 338. Jan. 2008, pp. xxii+973 (cit. on p. 25). G. Yang. Scaling Limits of Wide Neural Networks with Weight Sharing: Gaussian Process Behavior, Gradient Independence, and Neural Tangent Kernel Derivation. 2020. arXiv: 1902.04760 [cs.NE] (cit. on p. 5). G. Yang. Tensor Programs II: Neural Tangent Kernel for Any Architecture. 2020. arXiv: 2006.14548 [cs.LG] (cit. on pp. 1, 5, 10). G. Yang. Tensor Programs I: Wide Feedforward or Recurrent Neural Networks of Any Architecture are Gaussian Processes. 2021. arXiv: 1910.12478 [cs.NE] (cit. on pp. 1, 2, 4–6, 8, 9). G. Yang and E. J. Hu. Tensor Programs III: Neural Matrix Laws. 2020. arXiv: 2009.10685 [cs.LG] (cit. on pp. 1, 5).
REFERENCES
[YH21]
37
G. Yang and E. J. Hu. “Tensor Programs IV: Feature Learning in InfiniteWidth Neural Networks”. In: Proceedings of the 38th International Conference on Machine Learning. Ed. by M. Meila and T. Zhang. Vol. 139. Proceedings of Machine Learning Research. PMLR, 2021, pp. 11727–11737 (cit. on p. 5).
Institute for Mathematical Statistics and Actuarial Sciences, University of Bern, Switzerland Email address: [email protected] Department of Mathematics, University of Pisa, Italy Email address: [email protected] Department of Mathematics, University of Pisa, Italy Email address: [email protected]