Conceptio › Archive › arXiv CS

arXiv CSopen access

Universality of first-order methods on random and deterministic matrices

Unknown · 2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

machine learning, deep learning, neural networks

Universality of first-order methods on random and deterministic matrices Nicola Gorini∗

Chris Jones†

Dmitriy Kunisky‡

Lucas Pesenti§

arXiv:2604.11729v1 [math.PR] 13 Apr 2026

April 14, 2026

Abstract General first-order methods (GFOM) are a flexible class of iterative algorithms which update a state vector by matrix-vector multiplications and entrywise nonlinearities. A long line of work has sought to understand the large-n dynamics of GFOM, mostly focusing on “very random” input matrices and the approximate message passing (AMP) special case of GFOM whose state is asymptotically Gaussian. Yet, it has long remained unknown how to construct iterative algorithms that retain this Gaussianity for more structured inputs, or why existing AMP algorithms can be as effective for some deterministic matrices as they are for random matrices. We analyze diagrammatic expansions of GFOM via the limiting traffic distribution of the input matrix, the collection of all limiting values of permutation-invariant polynomials in the matrix entries, to obtain the following results: 1. We calculate the traffic distribution for the first non-trivial deterministic matrices, including (minor variants of) the Walsh–Hadamard and discrete sine and cosine transform matrices. This determines the limiting dynamics of GFOM on these inputs, resolving parts of longstanding conjectures of Marinari, Parisi, and Ritort (1994). 2. We design a new AMP iteration which unifies several previous AMP variants and generalizes to new input types, whose limiting dynamics are Gaussian conditional on some latent random variables. The asymptotic dynamics hold for a large and natural class of traffic distributions (encompassing both random and deterministic input matrices) and the algorithm’s analysis gives a simple combinatorial interpretation of the Onsager correction, answering questions posed recently by Wang, Zhong, and Fan (2022).

∗

Bocconi University. [email protected] UC Davis. [email protected] ‡ Johns Hopkins University. [email protected] § ETH Zürich. [email protected] †

Contents 1 Introduction

1.1 Approximate message passing and simple effective dynamics . . . . . . . . . . . . . .

1.2

Our contributions: Combinatorial method for GFOM . . . . . . . . . . . . . . . . . .

1.3

Related work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

1.4

Organization of the paper . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

1.5 Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2 Preliminaries

2.1

Matrix notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.2

Modes of convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.3

Matchings and Wick calculus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

3 Diagrams and the w- and z-Bases of Polynomials

3.1

Classes of diagrams . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

3.2

Graph polynomials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

3.3

Partitions, change of basis, and Möbius inversion . . . . . . . . . . . . . . . . . . . . 18

3.4 The example of cycles: Moments versus free cumulants . . . . . . . . . . . . . . . . . 19 3.5

Solving equations in the traffic distribution . . . . . . . . . . . . . . . . . . . . . . . 20

3.6

Products and concentration of traffic observables . . . . . . . . . . . . . . . . . . . . 23

4 Traffic Distributions of Random Matrices

4.1 Wigner random matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 4.2

Orthogonally invariant random matrices . . . . . . . . . . . . . . . . . . . . . . . . . 24

4.3

Block-structured random matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

5 Universality for Deterministic Matrices 5.1

Calculation of cactus diagrams and diagonal distribution . . . . . . . . . . . . . . . . 29

5.2 The fundamental theorem of graph polynomials . . . . . . . . . . . . . . . . . . . . . 32 5.3

Main structural lemma: Open cactus decomposition . . . . . . . . . . . . . . . . . . 33

5.4 The effect of puncturing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 5.5

Support of the z-basis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

5.6

Support of the w-basis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

5.7

Putting everything together: Proof of Theorem 5.3 . . . . . . . . . . . . . . . . . . . 45

6 From Diagrams to Asymptotic GFOM Dynamics

6.1 Asymptotic limit of the vector diagrams . . . . . . . . . . . . . . . . . . . . . . . . . 46

6.2 The treelike AMP algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 6.3

Examples of state evolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

A Traffic Distributions via Feynman Diagrams

A.1 Calculation of the free energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 A.2 Calculation of general observables: Argument for Theorem 4.2

. . . . . . . . . . . . 73

A.3 Mathematical comments on the Feynman diagram method . . . . . . . . . . . . . . . 75 B Traffic Distributions via Weingarten Calculus

B.1 Weingarten formula for orthogonal matrices . . . . . . . . . . . . . . . . . . . . . . . 75 B.2 Möbius inversion on non-crossing partitions . . . . . . . . . . . . . . . . . . . . . . . 77 B.3 Tracial moments concentration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 B.4 Traffic distribution of orthogonally invariant matrices . . . . . . . . . . . . . . . . . . 79 B.5 Concentration of traffic observables . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 B.6 Traffic distribution of punctured orthogonally invariant matrices . . . . . . . . . . . 84 C Convergence of Stochastic Processes

C.1 Connection with convergence of the empirical distribution . . . . . . . . . . . . . . . 87 D Omitted Proofs

D.1 Combinatorial lemmas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 D.2 Handling empirical averages . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 D.3 Proof of Lemma 6.30 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 D.4 Proof of Lemma 6.31 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96

Introduction

Complex systems with a large number of simply interacting pieces underlie many natural processes and, more recently, have been studied in computer science in an effort to make sense of how simple machine learning algorithms can learn complex structures latent in large, semi-random input data. Iterative optimization algorithms making sequential updates can be viewed as dynamical systems, with the main task being to understand how the algorithm evolves over time and what properties the eventual output will have. When the size of these systems grows very large, a key insight from statistical physics is that the macroscopic properties of the system can simplify dramatically: As the size of a random, smoothly-interacting dynamical system grows, the effect of individual particles averages out, and the dynamical system’s trajectory approximately follows an asymptotic distributional equation. We refer to these distributional equations as (asymptotic) effective dynamics. We seek to prove this kind of theorem for discrete-time nonlinear iterative algorithms such as those used in modern optimization, statistics, and machine learning. Concretely, we study general first-order methods (GFOM) [CMW20, MW24] which take as input a symmetric matrix A ∈ Rn×n , maintain a vector state x ∈ Rn , and at each step can perform one of two possible operations: 1. either multiply the state by A:

xt+1 = Axt ,

2. or apply a function ft : Rt+1 → R componentwise to the previous states: xt+1 = ft (xt , . . . , x0 ), i.e., xt+1 [i] = ft (xt [i], . . . , x0 [i]) for each i ∈ [n] . The initial state will be either the deterministic all-ones vector x0 = 1, or a random Gaussian vector x0 ∼ N (0, I) independent of A. Without loss of generality, we may assume that these operations alternate, giving an iteration of the form xt+1 = Aft (xt , . . . , x0 ) . We fix some number of iterations t and view xt = xt (A) as the output of the algorithm. GFOM is a flexible computational model which is expressive enough to capture many types of gradient descent [CMW20, GTM+ 24] and message passing algorithms [FVRS22]. It may be viewed as a nonlinear version of the power method for estimating top eigenvectors. The alternation of linear and nonlinear steps also closely matches the structure of a feedforward neural network [CHS24]. One may view the structural restriction on GFOM as forcing xt viewed as a function of A to be permutation-equivariant: if we apply the same permutation to the rows and columns of A, then xt undergoes the same permutation, a natural condition of an algorithm’s not depending on the particular indexing of its inputs. GFOM and their special case of approximate message passing (AMP) are very popular algorithms for many statistical inference tasks and are known to perform optimally in various such settings [DMM09, Ran11, Mon12, RFSK16, BM11, FVRS22]. In these cases, an algorithm takes as 1

input not an arbitrary matrix A, but one that contains a corrupted observation of a signal (in a common example, the input A is a low-rank yy ⊤ plus independent random noise). GFOM have also been used as optimization algorithms in average-case settings without any such planted structures. For instance, they are the best known algorithms for optimizing quadratic forms with random coefficients over the non-negative orthant [MR15] (the non-negative PCA objective function), other convex cones [DMR14], and the hypercube [Mon19] (the Sherrington–Kirkpatrick Hamiltonian), all of which are NP-hard problems in the worst case. This situation is the main target of our analysis. We receive an input matrix A without any particular “signal” and wish to output x approximately solving an optimization problem parametrized by A, such as maximize ⟨x, Ax⟩ subject to x ∈ S

(1)

studied in the above references for various choices of the constraint set S ⊆ Rn . To view GFOM as an instance of the physical setting sketched above, we consider a growing sequence of matrices A = A(n) ∈ Rn×n , and think of the “particles” as being the coordinates xt [i] of xt ∈ Rn . To keep notation reasonable, while all of these objects depend on n, we omit the (n) superscript whenever possible. We analyze the empirical distribution of our particles, accessed by sampling a random coordinate of a vector: samp(x) := x[i] ∈ R for i ∼ Unif([n]) . In order to study a particle’s entire trajectory more generally, we may “stack” several vectors and define samp((x0 , . . . , xt )) := samp(x0 , . . . , xt ) := (x0 [i], . . . , xt [i]) ∈ Rt+1 for i ∼ Unif([n]). The analysis of GFOM hinges on the observation that these random variables often converge in distribution to certain limiting distributions. That is, for suitably nice test functions φ : Rt+1 → R, (n)

(n)

lim E φ samp(x0 , . . . , xt ) = lim E

n→∞

n 1X (n) (n) φ(x0 [i], . . . , xt [i]) = n i=1

∞ φ dν≤t ,

∞ . For example, we can analyze the objective function of a problem for some probability measures ν≤t like Eq. (1) in this way: given a GFOM to run for t iterations producing xt = xt (A), we extend it to xt+1 = Axt so that n 1 1X 1 E ⟨xt , Axt ⟩ = E ⟨xt , xt+1 ⟩ = E xt [i]xt+1 [i] , n n n i=1

a quantity accessible in the above formalism by a suitable choice of φ. We can also study the algorithm’s convergence by expanding n1 ∥xt − xt−1 ∥22 in the same way. ∞. The goal of an asymptotic effective dynamics is then to identify the asymptotic measures ν≤t Such a description is a natural first step to designing optimal GFOM for optimization problems: given an explicit description of the limiting performance of any GFOM, we then optimize the performance over all GFOM [CMW20, AMS21, MW25, Pes26].

The goal of this paper is to study the following three questions regarding effective dynamics: 1. Existence: What are minimal assumptions on the input matrices and the algorithm that ensure the existence of asymptotic effective dynamics? 2

2. Universality: What properties of the sequence of input matrices A(n) determine the asymptotic effective dynamics? In particular, how can we show that two sequences of A(n) share the same dynamics? 3. Explicit Calculation: What are the effective dynamics? In particular, for a given algorithm, ∞ for each fixed t ∈ N? how can one describe ν≤t

1.1

Approximate message passing and simple effective dynamics

The majority of results to date on effective dynamics for GFOM, including ours, are most useful for Approximate Message Passing (AMP) algorithms. Originating from physicists’ work on mean-field spin glass models [MPV87, DMM09], AMP algorithms are a special case of GFOM with very ∞ above on the last simple effective dynamics: each distribution νt∞ (the marginal distribution of ν≤t coordinate) is a Gaussian distribution, νt∞ = N (µt , σt2 ) , 2 ) in terms of (µ , σ 2 ), . . . , (µ , σ 2 ) via a formula known and the effective dynamics gives (µt+1 , σt+1 t t 0 0 as the state evolution equation. This gives a simple yet complete description of the leading-order behavior of an algorithm as n → ∞. In part due to the power afforded by such a description, AMP (and the closely related belief propagation, of which AMP is a limit in a suitable sense) has taken on an indispensable role in statistical physics [MPV87, MM09, CMP+ 23] and, more recently, in computational statistics [ZK16, FVRS22].

In fact, while the original appearances of AMP in statistical physics were intrinsically motivated, for statistics applications the simplicity of state evolution is so useful that a line of work has emerged trying to design GFOM that have Gaussian νt∞ and effective dynamics given by state evolution [JM13, BSK15, VSR+ 15, Fan22, ZWF24, LWF25]. The term “AMP” is now often used to describe any choice of GFOM for a given family of inputs A(n) that has these properties. While it is not clear that this should be the case a priori, a common fortuitous coincidence is that, for various problems, the best GFOM algorithms (in the sense of achieving optimal rates in estimation or inference tasks) happen to be in the special class of AMP. That is, in many cases, the GFOM with the simplest asymptotic effective dynamics are also the most useful in applications. Given the successes of AMP, it is a longstanding goal in the literature to identify AMP-like algorithms for as many different choices of inputs and input distributions as possible. Yet, even to go slightly beyond the simplest choices of matrices A(n) has proved challenging and subtle (e.g., random matrices with i.i.d. entries [JM13, BLM15], orthogonally invariant distributions [Fan22], or semi-random ensembles [DLS23, WZF22]). Constructing AMP algorithms in such settings involves carefully inserting so-called Onsager correction terms into the nonlinearities ft in ways that remain somewhat mysterious yet are crucial to obtain Gaussian limiting behavior. Here, we will present an approach to the analysis of GFOM that re-derives different existing variants of AMP in a unified way, derives AMP algorithms for new inputs (both random and deterministic), and offers new conceptual insights into the design of these algorithms and into the proof of their asymptotic effective dynamics, in particular giving a clear combinatorial explanation for the Onsager corrections mentioned above.

1.2

Our contributions: Combinatorial method for GFOM

We study GFOM by expressing them as vectors of polynomials in the entries of the input matrix. For this reason we focus on polynomial ft ; it is likely possible to treat more general nonlinearities by approximating them by polynomials (see Section 1.3 for some discussion). Definition 1.1. We call a GFOM as described above a polynomial GFOM (pGFOM) if all nonlinearities ft : Rt+1 → R are polynomials. Our approach is divided into two parts. The first is a “static” analysis of certain symmetric polynomials in the entries of the input A. The second translates this to “dynamic” information about vector-valued functions, allowing us to calculate effective dynamics for O(1) iterations of GFOM in a general way. 1.2.1

Statics of graph polynomials: Traffic distributions and universality

The basic objects of study for our static analysis are the following graph polynomials. Definition 1.2 (Diagram classes). We write A = A0 for the set of finite, undirected, connected (multi)graphs. We also write E = E0 ⊆ A0 for the set of 2-edge-connected (multi)graphs (ones that cannot be disconnected by removing any single edge) and C = C0 ⊆ E0 ⊆ A0 for the set of cactus graphs, ones where every edge belongs to exactly one simple cycle.1 See Figure 1. The optional subscript “0” of the diagram classes refers to the outputs of the polynomials being 0-dimensional, i.e., scalars, which will be useful to distinguish them from vector- and matrix-valued polynomials to be defined later (with subscript “1” and “2”, respectively).

Figure 1: A cactus graph in C. Intuitively, a cactus is a “tree of cycles”. Definition 1.3 (Scalar graph polynomials). Given α ∈ A and A ∈ Rn×n sym , define polynomials wα (A), zα (A) ∈ R[A] by: wα (A) =

A[i(u), i(v)] ,

i:V (α)→[n] {u,v}∈E(α)

zα (A) =

A[i(u), i(v)] .

i:V (α),→[n] {u,v}∈E(α) 1

This notion is sometimes more specifically called a bridgeless cactus; in this paper we take this to be part of the definition of a cactus.

That is, wα (A) and zα (A) are each multivariate polynomials in the n(n+1) entries on and above 2 the diagonal of the matrix A obtained by summing over all labelings of the vertices of α by [n] = {1, 2, . . . , n} and with each edge corresponding to an entry of A. The only difference between wα (A) and zα (A) is that the vertex labeling for zα (A) is restricted to be injective by the notation i : V (α) ,→ [n] whereas labels in wα (A) are allowed to repeat. Each monomial in the entries of A can be represented as a multigraph on {1, 2, . . . , n}. By summing all monomials with the same “shape”, the wα (A) and zα (A) give two different spanning sets for a subspace of the Sn -invariant polynomials in the entries of A, where Sn acts on A by permuting the rows and columns simultaneously. There are only a few possible distinct shapes for monomials with low degree, so analysis on the w or z polynomials is a highly compressed way to analyze Sn -invariant low-degree polynomial functions of A. The limiting values of the graph polynomials are a basic set of parameters for the sequence of matrices A(n) , introduced in random matrix theory by Male [Mal20], who termed them the traffic distribution. Definition 1.4 (Traffic distribution). For a sequence of random2 matrices A = A(n) ∈ Rn×n sym we say that D : A → R is the (limiting) traffic distribution of A if lim

E wα (A) = D(α) for all α ∈ A.

n→∞ n A

(2)

We say the (limiting) traffic distribution exists if the limit exists for all α ∈ A.3 When the limiting traffic distribution exists, it is easy to show that it determines the asymptotic behavior of all constant-time GFOM algorithms with input A: Claim 1.5. Assume that A = A(n) ∈ Rn×n sym have traffic distribution D, and that a pGFOM defines xt = xt (A) with x0 = 1. Then, for any fixed t and polynomial φ ∈ R[x], n 1X φ(xt [i]) = C, lim E n→∞ A n i=1

where C is a constant depending only on D, (fs )1≤s≤t , and φ. Because of this observation, the traffic distribution is a natural way both to show existence of effective dynamics for constant-time GFOM (when the traffic distribution exists then so do effective dynamics) and to characterize the universality class of GFOM (when two sequences of matrices have the same traffic distribution then they have the same effective dynamics). We now reach our first main contribution: by calculating their limiting traffic distributions, we obtain the first analysis of GFOM on non-trivial completely deterministic inputs. Namely, we prove that any delocalized orthogonal matrix, after a slight modification, has the same traffic distribution as a corresponding random matrix model, the regular random orthogonal model (r-ROM; see Definition 2.5). Theorem 1.6 (See Theorem 5.1). Let Π = Π(n) = I − n1 11⊤ and H = H (n) ∈ Rn×n sym be a sequence 2 3

Deterministic matrices are also allowed just by taking a constant distribution. Note that the diagram α cannot depend on n. It has constant size as n → ∞.

of orthogonal matrices such that 1

max |H[i, j]| ≤ n− 2 +o(1) .

1≤i,j≤n

(3)

Then, the traffic distribution of ΠHΠ exists and equals that of the r-ROM. The motivating examples for Theorem 1.6 are “Fourier transform matrices” such as the Walsh– Hadamard matrix (Definition 2.3), discrete sine transform matrix, or discrete cosine transform matrix (Definition 2.4). We call conjugating by the projection matrix Π puncturing the matrix. Theorem 1.6 implies that, after puncturing, the effective dynamics of GFOM on these matrices are the same as those for the r-ROM, which itself is a punctured version of the random orthogonal model (ROM) of [MPR94a]. Explicit state evolution equations for these dynamics are given in Theorem 6.29. Puncturing is necessary in Theorem 1.6 and is natural for Fourier transform matrices. For the Walsh–Hadamard matrix, puncturing removes the first row and column, all of whose entries are √ identically 1/ n. This row/column makes H1 have a single large entry; because of that imbalance, without puncturing the traffic distribution of H does not exist4 and some GFOMs do not have well-defined asymptotic dynamics. This phenomenon has also been observed experimentally: [Sch20] writes that “structured matrices (e.g., DCT, Hadamard, Fourier) should work as well as i.i.d. random ones. But, in practice, AMP often diverges with such structured matrices.” We propose, and our results corroborate, that it is precisely alignment with the all-ones vector that causes this behavior. Showing that Fourier transform matrices are pseudorandom orthogonal matrices has been a longstanding folklore open problem in the statistical physics and AMP literature. It seems to originate in the work of [MPR94a, MPR94b, PP95] in statistical physics, who proposed these matrices as couplings for spin glass models. Recently (nearly 30 years later), [DLS23] summarized the situation as follows: More generally, numerical studies reported in the literature [...] suggest that AMP algorithms exhibit universality properties as long as the eigenvectors are generic. Formalizing this conjecture remains squarely beyond existing techniques, and presents a fascinating challenge. Similar comments have been made in [JM12, RSFS19, BSK15], and relevant numerical experiments can be found in [ÇO19, ABKZ20, DLS23]. Fourier transform matrices are also favored in compressed sensing applications since they admit fast multiplications via the Fast Fourier Transform [WZF22, Example 2.26]. Although Theorem 1.6 concerns orthogonal matrices, we also prove generally that after puncturing, any sequence of delocalized matrices has the same traffic distribution as the orthogonally invariant ensemble with the same eigenvalue distribution, assuming stronger delocalization properties than Eq. (3). See Theorem 5.3 for the formal statement. 1.2.2

Cactus properties: conditions for simple traffic distributions

The traffic distribution is a complicated object in general, just because its indexing set A is very large. Fortunately, traffic distributions of many common matrices are much simpler. Specifically, 4

For example, when H is the Walsh–Hadamard matrix, the degree-D star diagram σD satisfies n1 |wσD (H)| = Θ(nD/2−1 ), which diverges for D > 2.

they often satisfy a cactus property: almost all of the graph polynomials zα (A) are asymptotically negligible as n → ∞, with the only exceptions being the cactus graphs α ∈ C ⊊ A (in the z basis, but not in the w basis). Definition 1.7 (Cactus properties and cactus type). For a sequence of symmetric matrices A = A(n) ∈ Rn×n , we say that: (i) A has the strong cactus property if limn→∞ n1 EA zα (A) = 0 for all α ∈ A \ C. (ii) A has the weak cactus property if limn→∞ n1 EA zα (A) = 0 for all α ∈ E \ C. (iii) A has the factorizing (strong or weak) cactus property if it has the (strong or weak) cactus Q property, and for each σ ∈ C we have limn→∞ n1 EA zσ (A) = ρ∈cyc(σ) κ|ρ| for some real numbers κq , where cyc(σ) is the set of cycles of a cactus and |ρ| is the length of a cycle.5 The idea that the non-negligible diagrams for many random matrix models are cactuses appeared in the physics literature as early as the 1990s [PP95, MFC+ 19] and we will show in Section A how it can be derived from the Feynman diagram expansion widely used in physics. More recent mathematical work [Mal20, CDM24] reviewed in Section 4 has rigorously established the strong cactus property for Wigner matrices and unitarily invariant matrices whose eigenvalue distributions converge weakly. In fact, the factorizing strong cactus property is essentially equivalent to A having the same limiting traffic distribution as some orthogonally invariant random matrix model. The strong cactus property implies that the traffic distribution is specified only by the limiting values associated to σ ∈ C, a much smaller set of graphs than A. Another way to say this is that, under the strong cactus property, the traffic distribution contains no extra information beyond the considerably simpler diagonal distribution, introduced by [WZF22]. Definition 1.8 (Diagonal distribution). For a sequence of symmetric matrices A = A(n) ∈ Rn×n , we say that D : C → R is the limiting diagonal distribution of A if 1 E wσ (A) = D(σ) for all σ ∈ C. n→∞ n A lim

We say the diagonal distribution exists if the limit exists for all σ ∈ C. Let us make several important observations about the definitions of the traffic distribution, the diagonal distribution, and the cactus properties. First, note that Definition 1.7 is stated in the z-polynomial basis, whereas Definitions 1.4 and 1.8 are stated in the w-polynomial basis. Throughout the paper, it will be helpful to move back and forth between these bases, since some properties are most natural (or even are only true) in one basis or the other. This can be done via Möbius inversion, as described in Section 3.3. Second, neither the diagonal distribution nor the traffic distribution is an actual probability distribution. Instead, they should be interpreted as specifying limiting moments of certain empirical distributions, namely, the empirical distributions of the entries of vector graph polynomials.6 5

In the traffic probability literature, the factorizing strong cactus property has been referred to as a traffic distribution being of cactus type [CDM24]. The parameters κq are the free cumulants appearing in free probability theory. 6 The reason for the name of the diagonal distribution D is that it can also be interpreted as specifying the moments of the empirical distribution over the diagonal of certain matrices, namely those that can be formed from A by matrix multiplication and the operation of zeroing out the off-diagonal entries of a matrix [WZF22].

Third, one can view the diagonal and traffic distributions as generalizations of the limiting spectral distribution of a sequence of matrices. The spectral moments are n1 Tr(Aq ) = n1 wα (A), where α is the q-cycle diagram, so they are included in both the diagonal and traffic distributions: “ spectral distribution ⊊ diagonal distribution ⊊ traffic distribution ” Just as the empirical spectral distribution characterizes the limiting behavior of all polynomials in A that are invariant under the action of the orthogonal group O(n) (acting by Q · A = QAQ⊤ ), the traffic distribution characterizes the limiting behavior of the larger space of polynomials invariant under the smaller symmetric group Sn , i.e., where Q is restricted to be a permutation matrix. Finally, the strong cactus properties describe when these inclusions can be reversed: if the strong cactus property holds, then the traffic distribution contains no more information than the diagonal distribution. If the factorizing strong cactus property holds, then the diagonal distribution, in turn, contains no more information than the spectral distribution. Due to the effect of the puncturing operation, the strong cactus property actually is not satisfied by the pseudorandom matrices or r-ROM matrices appearing in our Theorem 1.6. But, these matrices satisfy the weak cactus property, and establishing this is a key step in the analysis of these matrices (in fact, the weak cactus property holds for the Fourier transform matrices without puncturing, as we show in Part 1 of Theorem 5.3). 1.2.3

Dynamics of graph polynomials: asymptotic GFOM state and treelike AMP

Recall that our final goal is to describe the state xt = xt (A) of a GFOM. Since xt ∈ Rn we use vector diagrams for this task. Compared to the scalar diagrams in A0 , the only extra information in these diagrams is that one of the vertices is specially marked as the “root”, whose label specifies the coordinate of the vector output. Definition 1.9 (Vector diagram classes). We write A1 and C1 for the set of graphs in A and C respectively, further decorated with a distinguished root vertex. For α ∈ A1 , we write root(α) ∈ V (α) for its root vertex. Definition 1.10 (Vector graph polynomials). Given α ∈ A1 and A ∈ Rn×n sym , define vectors of n polynomials wα (A), zα (A) ∈ (R[A]) by, wα (A)[i] :=

A[j(u), j(v)] ,

j:V (α)→[n] {u,v}∈E(α) j(root(α))=i

zα (A)[i] :=

A[j(u), j(v)] ,

j:V (α),→[n] {u,v}∈E(α) j(root(α))=i

for all i ∈ [n]. To analyze the vector graph polynomials, we compute the moments of the empirical distribution of their entries. We will see that these are matched (asymptotically) by a family of scalar random variables Zα∞ , so the empirical distribution of the entries of zα (A) converges in a suitable sense to Zα∞ as n → ∞. Further, when A has the strong cactus property, an analogous property is inherited by these limiting distributions, only a small number of α ∈ A1 having a non-negligible limit. 8

Definition 1.11 (Treelike diagrams). We say that α ∈ A1 is treelike if it is a tree with hanging cactuses attached to the leaves of the tree (see Figure 2). We denote the set of treelike diagrams by T1 , and denote by G1 ⊆ T1 the set of treelike diagrams in which, after removing hanging cactuses, the root has degree exactly 1.

(a) A treelike diagram in T1

(b) A Gaussian diagram in G1

Figure 2: Examples of treelike and Gaussian diagrams. The root vertex is circled. Theorem 1.12 (Vector polynomial limits; see Theorem 6.2). Assume that A = A(n) has the strong cactus property with limiting diagonal distribution D. Assume also that the sequence of random variables (∥A(n) ∥)n≥1 is tight,7 i.e., that

for all ε > 0 there exists K > 0 such that sup Pr ∥A(n) ∥ > K ≤ ε . n≥1

(4)

Write zA1 (A) ∈ (RA1 )n for the stacking of values of all zα (A) for α ∈ A1 . Then, (d)

samp(zA1 (A)) −−−→ (Zα∞ )α∈A1 , n→∞

for a family of (partially dependent) random variables (Zα∞ )α∈A1 such that Zα∞ = 0 for all α not treelike, and which can be sampled as follows for α ∈ T1 : 1. Draw (Zσ∞ )σ∈C1 from a distribution determined by D. 2. Draw (Zγ∞ )γ∈G1 ∼ N (0, Σ∞ ) from a centered Gaussian distribution with countably infinite covariance matrix Σ∞ depending on (Zσ∞ )σ∈C1 . 3. Set (Zα∞ )α∈T1 \(G1 ∪C1 ) to be certain deterministic polynomial functions of (Zα∞ )α∈G1 ∪C1 . We note that samp(zA1 (A)) is a random variable taking values in RA1 , a countable product space. Thus, its convergence in distribution is the same as convergence in distribution of any finite-dimensional projection; see Section C. The application to pGFOM is as follows. Analogously to Claim 1.5, it is easy to see that the 7

If the matrices A(n) are deterministic, this should be understood as (∥A(n) ∥)n≥1 being bounded.

iterates xt (A) of a pGFOM admit a diagrammatic expansion of the form xt (A) =

ct,α zα (A) ,

(5)

α∈A1

for finitely supported coefficients (ct,α )α∈A1 . Given the limits of the individual diagrams above, for a given GFOM, number of iterations t, and coefficients as in Eq. (5), we write 

(X0∞ , . . . , Xt∞ ) := 

 X

c0,α Zα∞ , . . . ,

α∈A1

ct,α Zα∞  ,

α∈A1

a random variable in Rt+1 that describes the joint empirical distribution of the first t steps of the GFOM. We call this the asymptotic state of a GFOM (Definition 6.16). By Theorem 1.12, the asymptotic state describes limiting empirical averages over the GFOM states, in the sense that n 1X φ(x0 [i], . . . , xt [i]) = E φ(X0∞ , . . . , Xt∞ ) n→∞ A n i=1

lim E

for any φ : Rt+1 → R either a polynomial or a bounded continuous function (Lemma 6.17). In particular, if the only nonzero ct,α in Eq. (5) are non-treelike α or treelike α ∈ G1 , then the GFOM has an asymptotic state that is Gaussian conditional on (Zσ∞ )σ∈C1 . This observation leads to our second main contribution: a new family of treelike AMP algorithms simultaneously generalizing Orthogonal Approximate Message Passing (OAMP) algorithms [RSF19, Fan22] for orthogonally invariant matrices, and Generalized Approximate Message Passing (GAMP) algorithms [Ran11, JM13] for matrices with independent entries that are not necessarily identically distributed.8 Theorem 1.13 (Treelike AMP; see Theorem 6.18). Assume that A = A(n) satisfies the assumptions of Theorem 1.12. Given polynomial functions ft : R → R, define the pGFOM: x0 := 1 ,

xt := Aft−1 −

t−1 X

bs,t · fs ,

(The product bs,t · fs is entrywise.)

s=0

ft := ft (xt ) , bs,t [i] :=

ft′ := ft′ (xt ) .

n X

′ ′ ′ A[is , it−1 ]ft−1 [it−1 ]A[it−1 , it−2 ]ft−2 [it−2 ] · · · fs+1 [is+1 ]A[is+1 , is ] .

is ,...,it−1 =1 distinct is =i

Then, for any fixed t as n → ∞, the asymptotic state (X1∞ , . . . , Xt∞ ), conditional on (Zσ∞ )σ∈C1 , is a centered Gaussian vector. A formula for its covariance is given in Proposition 6.26. The subtracted terms t−1 s=0 bs,t · fs generalize the “Onsager correction terms” appearing in different variants of AMP. Theorem 1.13 and its proof address two questions posed in [WZF22], namely (1) to obtain a combinatorial interpretation of the Onsager correction for OAMP algorithms, and (2) to identify a more general class of AMP algorithms whose state evolution is characterized by the diagonal distribution of the input matrix. Theorem 1.13 shows that (2) is possible for arbitrary P

The second comparison is with the caveat that GAMP uses a certain class of “non-separable” nonlinearities (applying a different function ft to each coordinate of xt ) which are not directly covered by our result [Ran11, JM13].

matrices satisfying the strong cactus property, and explicitly describes such an algorithm and its conditionally Gaussian asymptotic states. We show in Section 6.3 how the treelike AMP iteration simultaneously generalizes several variants of AMP introduced in prior work. We emphasize that, in contrast to all existing state evolution results we are aware of, we derive an Onsager correction and state evolution formula without assuming an explicit random model for A. The iteration in Theorem 1.13 is the same regardless of the limiting diagonal distribution of A, provided that these matrices (random or deterministic) satisfy the strong cactus property and have some limiting diagonal distribution (which will affect the covariance formula in Proposition 6.26). Note that the matrices in our universality result (Theorem 1.6) and their random counterparts (the r-ROM), satisfy the weak cactus property instead of the strong cactus one. Nevertheless, the Onsager correction and the state evolution can still be determined by a reduction to the strong-cactus-property setting, as we explain in Section 6.3.2.

1.3

Related work

Moment method for AMP. Our overall approach to graph polynomials generalizes prior work for the case of Wigner matrices [JP25]. Similar techniques have also appeared in prior works using the moment method to study AMP algorithms [BLM15, WZF22, MW25, DLS23, IS24, DSL24]. The w and z polynomials are rather fundamental objects which, along with their vector, matrix, and tensor generalizations, have variously been called “graph monomials” or “traffics” in free probability, “graph matrices” in computer science, “graph homomorphism polynomials” in combinatorics, and are also related to “tensor networks” and “Feynman diagrams” in physics. Polynomial vs. non-polynomial GFOM. In random and semi-random models, general firstorder methods with a constant number of iterations using (1) only polynomial nonlinearities or (2) arbitrary Lipschitz nonlinearities are generally expected to have the same computational power. Using polynomial approximation arguments, this has been made precise in several previous works [MW25, IS24, WZF22]. For example, [WZF22, Lemma 2.12] gives an abstract reduction showing that if state evolution for AMP on rotationally-invariant matrices holds for polynomial nonlinearities, then it also holds for arbitrary Lipschitz nonlinearities. While we study more general matrix models, we expect the assumption of polynomial nonlinearities is not essential. AMP vs. GFOM. A simple reduction shows that every algorithm in the GFOM class can be expressed as a certain post-processing of an AMP algorithm (allowing “memory terms”) [CMW20]. Therefore, these two classes of algorithms are equivalent from the standpoint of computational power. In our analysis, this is mirrored by the fact that, in Theorem 1.12, all possible non-Gaussian limits after conditioning on the draw of (Zσ∞ )σ∈C1 are deterministic functions of the possible Gaussian limits. GFOM on independent entry matrices. The analysis of GFOM and AMP on Wigner matrices or inhomogeneous versions thereof was the first case widely considered in the literature, and goes back to the origins of the mathematical analysis of AMP in the statistical physics literature on spin glasses [Bol14, DMM09, BM11, Mon12, BSK15, RV18, LW22]. See [FVRS22] for a survey of many of these works. Further, see [BLM15, CL21] for universality results over such models allowing for different entry distributions (but still requiring entrywise independence), [DJM13, JM13] for results 11

on block-structured variance profiles along the lines of our block GOE model, and [GHN26, BHX25] for recent progress on more general variance profiles. GFOM on orthogonally invariant matrices. The correct form of AMP (to ensure Gaussian limiting distributions) in orthogonally invariant models was first predicted non-rigorously for physics applications by [OÇW16] using dynamical mean-field theory (DMFT), and then proved by [Fan22]. Precursors for special “divergence-free” forms of AMP were also obtained by [ÇO19, MP17, RSF19, Tak19] under the names of Vector AMP and Orthogonal AMP. Related calculations for a more general statistical physics framework subsuming these AMP variants are carried out in [MFC+ 19]; in particular, this work includes special cases of and discusses the more general form of the calculations we detail in Section B. See the discussion in [Fan22] for a more thorough overview of these distinctions. Universality principles for GFOM. Beyond the above results, the main ones we are aware of that reduce the amount of randomness required for AMP are the recent works [WZF22, DLS23], which, modulo technical differences, both prove universality results over random matrices whose distribution is invariant under signed permutations. In other words, they treat broad classes of matrices provided that these are conjugated by random signed permutations, a considerable reduction in randomness from, e.g., conjugating by random Haar-distributed orthogonal or unitary matrices as in OAMP. Numerous experimental works have found universality phenomena for “sufficiently pseudorandom” deterministic matrices, but we are not aware of any rigorous results for completely deterministic matrices prior to our work. See discussion in [ÇO19, Sch20, ABKZ20, DLS23].

1.4

Organization of the paper

We give preliminaries on the matrices considered in this work and modes of convergence for our limiting theorem in Section 2. We introduce our definitions of diagrams and consequences of Möbius inversions for the traffic distribution in Section 3. In Section 4, to build intuition on traffic distributions, we describe them for several random matrix ensembles. Section 5 is dedicated to the proof of our first main result, the polynomial universality of delocalized deterministic matrices (Theorem 1.6). Section 6 details and proves the effective dynamics of GFOM under the strong cactus property (Theorems 1.12 and 1.13). We illustrate two viable approaches to computing the traffic distribution of orthogonally invariant matrix models: Section A is based on Feynman diagrams and Section B relies on Weingarten calculus. Section C provides background on convergence of stochastic processes, and Section D contains omitted proofs.

1.5

Acknowledgments

Thanks to Zhou Fan, Cynthia Rush, and Subhabrata Sen for helpful discussions over the course of this project. CJ was supported in part by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 101019547). LP’s work was supported by the Swiss National Science Foundation (SNSF), grant no. 10004947.

Preliminaries

2.1

Matrix notation

Given matrices A, B ∈ Rn×n , we will use: • A ∈ Rn×n sym to specify that A is symmetric. • A ∈ O(n) ⊆ Rn×n to specify that A is orthogonal. • A[i, j] to denote its (i, j)-th entry for i, j ∈ [n] := {1, . . . , n}. • ∥A∥ := max∥x∥2 =1 ∥Ax∥2 to denote its spectral or operator norm. • ∥A∥2F :=

2 i,j=1 A[i, j] to denote its Frobenius norm.

• Tr(A) :=

i=1 A[i, i] to denote its trace.

• λ1 (A) ≥ . . . ≥ λn (A) to denote its eigenvalues when A is symmetric. • A ⊙ B to denote the entrywise or Hadamard product with entries (A[i, j]B[i, j])i,j∈[n] . 1 ⊤ Definition 2.1 (Puncturing). Let H ∈ Rn×n sym and Π := I − n 11 be the projection orthogonal to the all-ones direction. The puncturing of H is the matrix A = ΠHΠ.

Definition 2.2 (GOE). The (normalized) Gaussian Orthogonal Ensemble GOE is the distribution of random matrices A ∈ Rn×n sym with A[i, j] = A[j, i] ∼ N (0, 1/n) independently for all 1 ≤ i < j ≤ n, and A[i, i] ∼ N (0, 2/n) independently for all i ∈ [n]. Definition 2.3 (Hadamard matrices). When n is a power of 2, the (normalized) Walsh–Hadamard (n) matrix Hhad ∈ Rn×n sym is defined recursively by (1) Hhad =

1 (2n) Hhad := √

(n)

Hhad Hhad (n) (n) . 2 Hhad −Hhad

h i

1 ,

√ (n) Hhad is a symmetric orthogonal matrix with entries in ±1/ n. (n)

Definition 2.4 (DST and DCT matrices). The discrete sine transform matrices Hsin ∈ Rn×n sym are s (n) Hsin [i, j] :=

2 πij sin n+1 n+1

∀i, j ∈ [n] .

(n)

The discrete cosine transform matrices Hcos ∈ Rn×n sym are (n) Hcos [i, j] :=

π(i − 12 )(j − 21 ) 2 cos n n

∀i, j ∈ [n] .

√ (n) (n) Hcos and Hsin are symmetric orthogonal matrices with entries at most O(1/ n) in magnitude. Definition 2.5 (ROM and r-ROM). The Random Orthogonal Model ROM is the distribution of random matrices H = QDQ⊤ , where Q ∈ O(n) is Haar-distributed, and D is a diagonal matrix 13

with i.i.d. Unif({−1, 1}) entries, independent from Q. The Regular Random Orthogonal Model r-ROM is the distribution of the puncturing of H, when H is sampled from the ROM. Random matrices from the ROM are symmetric orthogonal matrices, satisfying H 2 = I. They are a special case of the orthogonally invariant models we discuss in Section 4.2.

2.2

Modes of convergence

We will use a few standard modes of convergence from scalar-valued probability theory. Definition 2.6 (Modes of convergence: scalars). For a sequence of random variables x(n) ∈ R, we say that: • x(n) converge in expectation if, for some c ∈ R, limn→∞ Ex(n) = c. • x(n) converge in probability if, for some c ∈ R, for all ε > 0, limn→∞ P[|x(n) − c| > ε] = 0. • x(n) converge in L2 if they converge in expectation and limn→∞ E(x(n) −c)2 = 0, or equivalently if they converge in expectation and limn→∞ Var x(n) = 0. We write a symbol M ∈ {E, P, L2 } to indicate these modes of convergence, and in this notation say that the x(n) converge in M. Moreover, we say a sequence of random vectors x(n) ∈ Rd in fixed dimension d ≥ 1 converges in distribution to a random vector x ∈ Rd if for every bounded continuous function φ : Rd → R, E φ(x(n) ) −→ E φ(x) , n→∞

(d)

in which case we write x(n) −→ x. See Section C for a generalization to random variables indexed by a countably infinite index set. Definition 2.7 (Modes of convergence: tracial moments). For a mode of convergence M, we say that a sequence of random matrices A ∈ Rn×n converges in tracial moments in M if, for every k ≥ 1, n1 Tr Ak converges in M. We say that it converges in tracial moments in M to a probability measure µ over R if Z 1 Tr Ak → xk dµ(x) n in the mode of convergence M.

2.3

Matchings and Wick calculus

Given a set S, let M(S) denote the set of matchings on S. Let Mperf (S) denote the subset of perfect matchings. The elements of M ∈ M(S) are written as pairs {i, j} ⊆ S. For several sets S1 , . . . , Sk , denote by M(S1 , . . . , Sk ) the set of matchings on the disjoint union S1 ⊔ · · · ⊔ Sk that do not match any two elements of the same Si . For two sets S1 , S2 of the same size, denote by Mperf (S1 , S2 ) the bipartite perfect matchings of S1 ⊔ S2 that only match elements of S1 to ones of S2 . We will abbreviate M({1, 2, . . . , q}) as M(q).

Lemma 2.8 (Wick lemma). Let X1 , . . . , Xq be jointly Gaussian random variables with mean zero. Then: X Y E[X1 · · · Xq ] = E[Xi Xj ] . M ∈Mperf (q) ij∈M

The Wick products are the multivariate generalization of the Hermite polynomials to correlated Gaussians [Jan97, Chapter 3]. Definition 2.9 (Wick product). Let I be an index set, X = (Xi )i∈I be formal variables, and I Σ ∈ RI×I sym . The Wick products are defined by, for each finitely supported α ∈ N , Heα (X ; Σ) :=

(−1)|M |

Y uv∈M

M ∈M(α)

Σ[u, v]

Xu ,

u∈M /

where M(α) denotes the set of matchings on a collection consisting of αi copies of each i ∈ I. When |I| = 1, X ∼ N (0, 1), and Σ = 1, then He(p) (X ; Σ) equals the pth Hermite polynomial. When the Xi are mean-zero Gaussian random variables and Σ is their covariance matrix, the Wick products satisfy the (partial) orthogonality property that for each finitely supported α, β ∈ NI P P with i αi ̸= i βi , E [Heα (X ; Σ) Heβ (X ; Σ)] = 0 . In general, we have E [Heα (X ; Σ) Heβ (X ; Σ)] =

Σ[u, v] .

M ∈Mperf (α,β) uv∈M

Since by the Wick lemma E i∈α Xi · j∈β Xj equals the same sum over all matchings of α ⊔ β, the Wick products achieve a general “partial orthogonalization” that removes all terms from this covariance where any pairs within α or within β are matched. Q

For each choice of Σ ∈ RI×I sym , the Wick products are a basis for polynomials in the Xi . Multiplication of polynomials gives an algebra structure to this space which we call the Wick algebra of X. Below is a combinatorial formula for multiplication in the Wick algebra. Proposition 2.10 ([Jan97, Theorem 3.15]). Let I be an index set, X = (Xi )i∈I be formal variables, 1 k I and Σ ∈ RI×I sym . Let α , . . . , α ∈ N . Then: k Y

Heαj (X; Σ) =

j=1

Σ[u, v] HeU (M ) (X ; Σ) ,

M ∈M(α1 ,...,αk ) uv∈M

where αj is a multiset of size |αj | with αij copies of each i ∈ I. Here U (M ) ∈ NI for M a matching of α1 ⊔ · · · ⊔ αk counts the number of unmatched elements of each type. In the special case where each group αj consists of a single element, we obtain: Corollary 2.11. For every i1 , . . . , ik ∈ I, k Y j=1

Xij =

Σ[iu , iv ] HeU (M ) (X ; Σ) .

M ∈M(k) uv∈M

Diagrams and the w- and z-Bases of Polynomials

All graphs considered in this paper are multigraphs (loops and multiedges are allowed) and will be denoted by Greek letters (α, β, γ, . . .). We use the terms graphs and diagrams interchangeably in this paper. Given a diagram α, we use V (α) to denote its vertex set and E(α) to denote its edge set. We denote by α[S] the subgraph of α induced by S ⊆ V (α). We count self-loops as contributing 2 to the degree of a vertex.

3.1

Classes of diagrams

Each diagram can have either 0, 1, or an ordered pair of 2 special vertices called its root(s). With the exception of the class of graphs defined in Definition 5.4, the roots of a graph can be arbitrary vertices (in particular, they might be equal if there are two of them). Notation 3.1. Let A = A0 (resp. A1 or A2 ) be the set of all connected graphs with no root (resp. 1 root or 2 roots). We also refer to such graphs as scalar (resp. vector or matrix) diagrams. Given α ∈ A, an edge e ∈ E(α) is a bridge of α if deleting e would disconnect the graph. α ∈ A is 2-edge-connected if it contains no bridge. In general, α ∈ A can be decomposed into a tree of 2-edge-connected components connected by bridges. Notation 3.2. Let E = E0 ⊆ A (resp. E1 ⊆ A1 or E2 ⊆ A2 ) be the set of all 2-edge-connected scalar (resp. vector or matrix) diagrams. Given α ∈ A, a vertex u ∈ V (α) is an articulation point of α if removing u and its incident edges disconnects the graph. α is 2-vertex-connected if it has no articulation point. Any α ∈ A decomposes into its 2-vertex-connected components (blocks), which refine the 2-edge-connected components. The block-cut graph (whose vertices are the articulation points and the blocks, with edges for incidence) is a tree. A connected graph is a cactus if every edge lies on exactly one simple cycle. Thus, cactuses are in a sense the minimal 2-edge-connected graphs. Notation 3.3. Let C = C0 ⊆ A (resp. C1 ⊆ A1 ) be the set of all scalar (resp. vector) cactus diagrams. For a cactus σ, we will denote by cyc(σ) the set of (unrooted) cycles of σ. Finally, as in Definition 1.11, we will denote the treelike diagrams by T1 and the treelike diagrams such that the root has degree 1 after deleting all hanging cactuses by G1 .

3.2

Graph polynomials

Each diagram represents different scalar-, vector-, or matrix-valued polynomials in a matrix input, depending on whether it is viewed in the w-basis or the z-basis. In the following definitions, we fix A ∈ Rn×n sym , α to be a scalar, vector, or matrix diagram, and i, j ∈ [n].

Definition 3.4. Define wα (A) ∈ R, wα (A) ∈ Rn , and Wα (A) ∈ Rn×n by wα (A) =

A[φ(u), φ(v)]

if α is a scalar diagram,

A[φ(u), φ(v)]

if α is a vector diagram with root r,

φ : V (α)→[n] {u,v}∈E(α)

wα (A)[i] =

φ : V (α)→[n] {u,v}∈E(α) φ(r)=i

Wα (A)[i, j] =

A[φ(u), φ(v)]

if α is a matrix diagram with roots (r1 , r2 ).

φ : V (α)→[n] {u,v}∈E(α) φ(r1 )=i,φ(r2 )=j

Definition 3.5. Define zα (A) ∈ R, zα (A) ∈ Rn , and Zα (A) ∈ Rn×n by zα (A) =

A[φ(u), φ(v)]

if α is a scalar diagram,

A[φ(u), φ(v)]

if α is a vector diagram with root r,

φ : V (α),→[n] {u,v}∈E(α)

zα (A)[i] =

φ : V (α),→[n] {u,v}∈E(α) φ(r)=i

Zα (A)[i, j] =

A[φ(u), φ(v)]

if α is a matrix diagram with roots (r1 , r2 ).

φ : V (α),→[n] {u,v}∈E(α) φ(r1 )=i,φ(r2 )=j

The only difference between the w- and z-bases is the summation domain: Definition 3.5 sums over injective embeddings φ, whereas Definition 3.4 sums over all embeddings. Finally, we define two extensions of Definition 3.4 that we will need in the proofs. The following allows us to use a different matrix on each edge of the graph: Definition 3.6. Let α be a matrix diagram with roots (r1 , r2 ) and A = (Ae )e∈E(α) be such that n×n by Ae ∈ Rn×n sym for all e ∈ E(α). Define Wα (A) ∈ R Wα (A)[i, j] =

Ae [φ(u), φ(v)] .

φ : V (α)→[n] e={u,v}∈E(α) φ(r1 )=i,φ(r2 )=j

The following is an intermediate quantity between Definition 3.4 and Definition 3.5 which only restricts the sum over injective labelings on two vertices: Definition 3.7. Let A ∈ Rn×n sym , α be a scalar/vector/matrix diagram, i, j ∈ [n], and s, t ∈ V (α).

Define wαs̸=t ∈ R, wαs̸=t ∈ Rn , and Wαs̸=t ∈ Rn×n by wαs̸=t (A) =

A[φ(u), φ(v)]

if α is a scalar diagram,

A[φ(u), φ(v)]

if α is a vector diagram with root r,

A[φ(u), φ(v)]

if α is a matrix diagram with roots (r1 , r2 ).

φ : V (α)→[n] {u,v}∈E(α) φ(s)̸=φ(t)

wαs̸=t (A)[i] =

φ : V (α)→[n] {u,v}∈E(α) φ(s)̸=φ(t) φ(r)=i

Wαs̸=t (A)[i, j] =

φ : V (α)→[n] {u,v}∈E(α) φ(s)̸=φ(t) φ(r1 )=i φ(r2 )=j

3.3

Partitions, change of basis, and Möbius inversion

While (zα (A))α∈A and (wα (A))α∈A span the same space of Sn -invariant polynomials in the entries of A, some properties are better expressed in one basis than the other. Here we take a closer look at these bases and derive change-of-basis formulas. Given a set S, let P(S) denote the set of all partitions of S, sets of non-empty disjoint subsets of S whose union is all of S. We call the parts of a partition blocks. Each block is a set, and P is the set of blocks, so we denote the blocks by b ∈ P . For a (scalar, vector, or matrix) diagram α and a partition P ∈ P(V (α)), we define a new diagram αP by identifying the vertices within each block of P into a single vertex. The vertices of αP may thus be identified with the blocks of P . αP retains all edges of α, which may become multiedges or self-loops. The status of being one of the (0, 1, or 2) roots of α is inherited by the block containing that root. To change from the w- to the z-basis, we then simply sum over all partitions: Claim 3.8. For all (scalar, vector, or matrix) diagrams α, wα (A) =

zαP (A) .

X P ∈P(V (α))

Define the relation α ⪯ β on scalar diagrams if there exists a partition P ∈ P(V (β)) such that α = βP . It is easy to check that this relation gives a partial ordering, inherited from the standard partial ordering on partitions. We write α ≺ β as a shorthand for α ⪯ β and α ̸= β. Lemma 3.9. There exist (cα,β )α,β∈A and (c′α,β )α,β∈A not depending on n such that cα,β ∈ N, c′α,β ∈ Z and for any α, β ∈ A, wβ (A) =

cα,β zα (A) ,

zβ (A) =

α⪯β

c′α,β wα (A) .

α⪯β

Proof. The coefficients in the left equation count symmetries in Claim 3.8, i.e., cα,β equals the number of ways to choose a partition P ∈ P(V (β)) such that βP is isomorphic to α. Reciprocally, since ⪯ is a partial ordering, this transformation can be inverted using Möbius inversion [Rot64] on

this poset. Although an explicit formula for c′α,β is available in terms of the combinatorial structure of the graphs, we will not need it in this paper.

3.4

The example of cycles: Moments versus free cumulants

The difference between the w- and z-bases is illustrated nicely by the special case of the diagrams σq which are cycles of length q ≥ 1. In this case, n1 wσq (A) and n1 zσq (A) are versions of the limiting spectral moments and free cumulants, respectively, for finite-dimensional matrices. Let P(q) denote the set of partitions of {1, 2, . . . , q} and let NC(q) denote the subset of noncrossing partitions (partitions such that there does not exist i < j < k < ℓ with i, k in the same block and j, ℓ in the same block, different from the one i, k are in). It is convenient to view these as partitions of the vertices of the q-cycle so that the term non-crossing may be interpreted visually: in a non-crossing partition, the blocks do not intersect one another when drawn as “blobs” inside the cycle. In the w-basis, we have n n 1 1X 1 X 1 q wσ (A) = A[i1 , i2 ]A[i2 , i3 ] . . . A[iq , i1 ] = Tr(A ) = λi (A)q . n q n i ,...,i =1 n n i=1 1

(6)

Suppose that the expression in Eq. (6) converges as n → ∞ to the qth moment mq ∈ R of a limiting R spectral distribution, mq = λq dµ(λ). The free cumulants are defined from the moments by a formula similar to the classical cumulants vis-à-vis the moments of a random variable: Definition 3.10 (Free cumulant). The free cumulants (κq )q≥1 corresponding to (mq )q≥1 are defined implicitly by: X Y mq = κ|b| . (7) σ∈NC(q) b∈σ

The κq can be computed explicitly in terms of the mq by applying Möbius inversion to Eq. (7); see Eq. (59). Analogous to Eq. (6) which is in the w-basis, it appears to be folklore9 that if A is drawn from an orthogonally invariant matrix ensemble with free cumulants (κq )q≥1 , then 1 E zσq (A) −→ κq . n→∞ n

(8)

The quantity n1 zσq (A) has also been called the qth injective trace of A. Below in Lemma 3.12, we prove Eq. (8) using a change of basis from w to z. For example, below are the parameters mq and κq for the GOE and the ROM, whose limiting empirical spectral distribution are the Wigner semicircle distribution and the Rademacher distribution, respectively. 1 2k Claim 3.11. Let Cat(k) := k+1 k be the kth Catalan number. For the GOE, the limiting spectral

This is for example explicitly stated in [MFC+ 19, Theorem 1 and Appendix D.1].

moments and free cumulants are: (

mq =

)

Cat(q/2) if q is even 0 if q is odd

(

κq =

1 if q = 2 0 otherwise

)

For the ROM, the limiting spectral moments and free cumulants are: (

mq =

3.5

1 0

if q is even if q is odd

)

(

κq =

(−1)q/2−1 Cat(q/2 − 1) 0

if q is even if q is odd

)

(9)

Solving equations in the traffic distribution

The traffic distribution is defined as the limiting values of all w-basis polynomials, but we show now how it can be derived from various combinations of limits of w- and z-basis polynomials. In our other arguments, we will also find it convenient to describe the traffic distribution of sequences of matrices (random or deterministic) using the two bases simultaneously. While Lemma 3.9 shows that we could in principle express all these results in a single basis, this would involve precisely tracking very complicated combinatorial coefficients (in fact, this was a major technical obstacle in previous diagrammatic analyses of AMP). As we have discussed, when a matrix satisfies the strong cactus property, its traffic distribution is determined by its values on the cactus diagrams (equivalently, by the diagonal distribution), and when it satisfies the factorizing strong cactus property, its traffic distribution is determined by the spectral distribution. We show that one can use either the w-basis or z-basis for these determinations. Lemma 3.12. Suppose that A = A(n) satisfies the weak cactus property, i.e., for all α ∈ E \ C, 1 E zα (A) −→ 0 . n→∞ nA Then the following are equivalent: (i) For all σ ∈ C there exists mσ ∈ R such that n1 EA wσ (A) −→ mσ . n→∞

(ii) For all σ ∈ C there exists kσ ∈ R such that n1 EA zσ (A) −→ kσ . n→∞

Furthermore, when they exist, (mσ )σ∈C and (kσ )σ∈C determine each other. The following are also equivalent: (i) There exist real numbers (mq )q∈N such that for all σ ∈ C, n1 EA wσ (A) −→

(ii) There exist real numbers (κq )q∈N such that for all σ ∈ C, n1 EA zσ (A) −→

ρ∈cyc(σ) κ|ρ| .

ρ∈cyc(σ) m|ρ| .

n→∞

Furthermore, when they exist, (mq )q∈N and (κq )q∈N are related by Eq. (7). We use the following observation which will be used repeatedly in Section 5: Lemma 3.13. If α ∈ E and β ⪯ α, then β ∈ E.

Proof of Lemma 3.13. By Menger’s theorem, a graph is 2-edge-connected if and only if there exist two edge-disjoint paths between every pair of distinct vertices. These paths are maintained when α is contracted into β. Proof of Lemma 3.12. (ii) =⇒ (i). Using Claim 3.8, 1 1 X E wσ (A) = cβ,σ E zβ (A) . A nA n β⪯σ Every diagram β ⪯ σ remains 2-edge-connected by Lemma 3.13. There are only finitely many terms in the sum, so we can directly take the n → ∞ limit and use the assumptions to obtain that P 1 β⪯σ cβ,σ kβ . n EA wσ (A) converges to Note that by the weak cactus property, the only asymptotically nonzero β ⪯ σ are when β is a Q cactus. Assuming furthermore that kβ = ρ∈cyc(β) κ|ρ| factors over the cycles of each cactus β we will derive the second part of the lemma. Using the more specific result of Claim 3.8, we have X 1 1 lim E zσP (A) E wσ (A) = n→∞ n A n→∞ n A P ∈P(V (σ))

lim

Since A has the weak cactus property and σ is a cactus, only the terms where σP is a cactus contribute. These are precisely the terms where P restricted to each cycle of σ is non-crossing. Given Pρ ∈ NC(V (ρ)) for each ρ ∈ cyc(σ), let us write P (Pρ : ρ ∈ cyc(σ)) for the partition obtained by composing these partitions of each cycle, and let us write, following our previous notation, cyc(ρPρ ) for the set of cycles created when the single cycle ρ is contracted according to Pρ . Then, we have =

lim

X Pρ ∈NC(V (ρ)) for each ρ∈cyc(σ)

E zσP (Pρ :ρ∈cyc(σ)) (A)

n→∞ n A

κ|π|

Pρ ∈NC(V (ρ)) ρ∈cyc(σ) π∈cyc(ρPρ ) for each ρ∈cyc(σ)



 ρ∈cyc(σ)



Y Y

κ|π| 

P ∈NC(V (ρ)) π∈cyc(ρP )

m|ρ| .

ρ∈cyc(σ)

Thus we have the claimed factorization. Further, the coefficients mq and κq indeed have the relation between moments and free cumulants from Eq. (7): mq =

κ|b| .

σ∈NC(q) b∈σ

(i) =⇒ (ii). This direction uses a recursive change of basis technique that will be very useful

in Section 5. Using Lemma 3.9 in both directions, we get 1 1 X ′ 1 X ′ E zσ (A) = cβ,σ E wβ (A) + c E wβ (A) A nA n β⪯σ n β≺σ β,σ A β∈C

β∈E\C

1 X ′ 1 X ′ X = cβ,σ E wβ (A) + c cα,β E zα (A) A A n β⪯σ n β≺σ β,σ α⪯β β∈C

β∈E\C





 1 X ′ 1 X   X ′ cβ,σ cα,β  E zα (A) cβ,σ E wβ (A) +  A A n β⪯σ n α≺σ  β∈E\C β∈C

α⪯β≺σ

Note that every diagram in this expansion remains 2-edge-connected by Lemma 3.13. Every contraction identifying a non-empty subset of vertices decreases the number of vertices in the graph, and the w and z bases coincide for 1-vertex graphs. Therefore, we can apply the same steps inductively on terms for which α ∈ C to finally obtain 1 X ′′ 1 1 X ′′ cβ,σ E wβ (A) + E zσ (A) = c E zβ (A) . A nA n β⪯σ n β⪯σ β,σ A β∈C

β∈E\C

for some coefficients {c′′α,β } independent of n. Take the n → ∞ limit to obtain X 1 c′′β,σ mβ , E zσ (A) = n→∞ n A β⪯σ

lim

β∈C

which finishes the proof of the first equivalence. Assuming furthermore that mβ factors over the cycles of each cactus β, then n1 EA zσ (A) also asymptotically factors over its cycles: n1 EA zσ (A) −→ Q ρ∈cyc(σ) κ|ρ| for some numbers κq . This is because the cactuses β ⪯ σ still only arise by contracting a separate non-crossing partition for each cycle of σ, and so we can perform the above recursive analysis separately inside each cycle. The following lemma shows that the properties of graph polynomials we will establish for delocalized deterministic matrices in Section 5 characterize their traffic distribution. We emphasize our use of a combination of assumptions on limits of the w- and z-bases that makes this formulation convenient. Lemma 3.14. Suppose that A = A(n) satisfies: 1. The weak cactus property, i.e., that for all α ∈ E \ C, n1 EA zα (A) −→ 0. n→∞

2. For all α ∈ A \ E, n1 EA wα (A) −→ 0. n→∞

3. For all σ ∈ C, there exists mσ ∈ R such that n1 EA wσ (A) −→ mσ . n→∞

Then the traffic distribution of A exists and only depends on {mσ : σ ∈ C}.

Proof. We want to show that for every α ∈ A, limn→∞ n1 EA wα (A) exists and only depends on {mσ : σ ∈ C}. By assumption, it suffices to prove it for α ∈ E \ C. By Lemma 3.9, 1 1 X E wα (A) = cβ,α E zβ (A) . A nA n β⪯α By Lemma 3.13, every β in the support of the sum is 2-edge-connected. If β ∈ C, then the value of limn→∞ n1 EA zβ (A) exists and only depends on {mσ : σ ∈ C} by Lemma 3.12. Otherwise, β ∈ E \ C, and limn→∞ n1 EA zβ (A) = 0 by assumption. This implies that limn→∞ n1 EA wα (A) exists and only depends on {mσ : σ ∈ C}, which concludes the proof. Note that, more generally, by Lemma 3.12, the same statement will hold with Condition 3 of Lemma 3.14 taken in terms of either the w- or z-basis.

3.6

Products and concentration of traffic observables

Recall that the traffic distribution specifies the limits of n1 EA wα (A) for all α ∈ A. In all of the random matrix models we consider, these expectations are highly concentrated. We say that the traffic distribution concentrates for A if the following property holds, studied in [Mal20]. Definition 3.15. Let A = A(n) ∈ Rn×n sym and assume that the traffic distribution of A exists. We say that the traffic distribution concentrates for A if for all k ≥ 2 and α1 , . . . , αk ∈ A, 

lim E 

n→∞ A

k Y 1

j=1



wαj (A) =

k Y

1 E wαj (A) . n→∞ n A j=1 lim

The case k = 2 and α1 = α2 = α of the definition specializes to the statement: Lemma 3.16. Let A = A(n) ∈ Rn×n sym have traffic distribution D. If the traffic distribution concentrates for A, then n1 EA wα (A) converges to D(α) in L2 . The full condition may be viewed as a strengthening of this straightforward notion of concentration. We note that the product of several w-basis polynomials is equivalent to taking the disjoint union of their diagrams: wα1 (A) · · · wαk (A) = wα1 ⊔···⊔αk (A). Therefore, Definition 3.15 says that the values of disconnected diagrams asymptotically factor over the components. This justifies defining A and the traffic distribution to include only connected diagrams. The following shows that concentration may equally well be considered in the z-basis. Lemma 3.17 ([Mal20, Lemma 2.9]). Let A = A(n) ∈ Rn×n sym and assume that the traffic distribution of A exists. The traffic distribution concentrates for A if and only if, for all k ≥ 2 and α1 , . . . , αk ∈ A, 

lim E 

n→∞ A

k Y 1

j=1



zαj (A) =

k Y

1 E zαj (A) . n→∞ n A j=1 lim

For vector diagrams, the componentwise or Hadamard product is wα1 (A) · · · wαk (A) = wα1 ⊕···⊕αk (A) , 23

where α1 ⊕ · · · ⊕ αk is the diagram formed by taking the disjoint union of α1 through αk and then identifying the roots together into a single root. We sometimes refer to this operation as grafting α1 , . . . , αk at the root.

Traffic Distributions of Random Matrices

As both a technical preliminary for our results and useful background, this section describes the traffic distributions of several common random matrix ensembles. A common theme is that all of these classical models satisfy the strong cactus property. Most of these results have appeared previously in the literature, though we provide some extensions and new interpretations.

4.1

Wigner random matrices

A Wigner matrix is a random symmetric matrix with i.i.d. entries on and above the diagonal. Changes to the diagonal entries such as setting them to zero (which is the convention used in some works), or taking the diagonal variances to be twice the off-diagonal ones (as in the GOE model), do not affect the results. The limiting traffic distribution of a sequence of Wigner matrices was derived by Male [Mal20], by generalizing the combinatorial proof of the semicircle limit theorem for the limiting spectral distribution [AGZ10]. The same result was re-discovered in [JP25] in the context of analyzing pGFOM on such matrices. Theorem 4.1 (Traffic distribution of Wigner matrices). Let ν be a probability measure on R with e(n) ∈ Rn×n have entries on and all moments finite, mean 0, and variance 1. For all n ≥ 1, let A sym e(n) . Then, for all α ∈ A, above the diagonal drawn i.i.d. from ν. Define A(n) := √1n A 

1 1 lim E zα (A(n) ) = n→∞ n 0

if α is a cactus of 2-cycles, otherwise.

The same result holds for normalized GOE matrices. Note that a cactus of 2-cycles may equivalently be viewed as a “doubled tree”, a tree where every edge is repeated exactly twice, which is the formulation used in the previous works [Mal20, JP25]. Thus, sequences of Wigner matrices have the factorizing strong cactus property, with the especially simple sequence of free cumulants κ2 = 1 and κq = 0 for all q = ̸ 2. These are also the free cumulants of the semicircle law, which is the limiting eigenvalue distribution of A(n) .

4.2

Orthogonally invariant random matrices

⊤ Let the orthogonal group O(n) act on Rn×n sym by conjugation, with Q ∈ O(n) acting as Q·A := Q AQ. n×n Let µ denote a probability measure on Rsym that is invariant under this action of O(n). In this case, we call A ∼ µ an orthogonally invariant random matrix. n×n If µ has a density on Rn×n sym , an equivalent condition is that the density at A ∈ Rsym depends only on the unordered multiset of eigenvalues of A. An important class of examples in physics is given by matrix models with potential V : R → R, whose density is proportional to exp(− Tr V (A)).

For example, the GOE model corresponds to V (t) = t2 /2. We will come back to these examples in Section A. For the complex-valued variant where O(n) is replaced by the unitary group U (n), the limiting traffic distribution of such unitarily invariant random matrices is described in [CDM24, Theorem 1.1]. The same description holds in the orthogonal case. The proof is a straightforward generalization of the unitarily invariant case, but for the sake of completeness we present it in detail in Section B. Theorem 4.2 (Traffic distribution of orthogonally invariant random matrices). Let A(n) ∈ Rn×n sym be a sequence of orthogonally invariant random matrices that converges in tracial moments in L2 to a probability measure µ. Then, for all α ∈ A,   

1 E zα (A(n) ) = σ∈cyc(α) n→∞ n   0 lim

κ|σ|

if α ∈ C, (10) otherwise.

where κq is the qth free cumulant of µ (Definition 3.10), and |σ| denotes the length of the cycle. Eq. (10) shows that the factorizing strong cactus property holds for orthogonally invariant random matrices, and in particular their limiting traffic distribution is supported only on cactus diagrams in the z-basis. Actually, in this case the strong cactus property is non-trivial only for the Eulerian diagrams, since the non-Eulerian ones have identically zero expectation for each fixed dimension n: Claim 4.3. Let A(n) ∈ Rn×n sym be an orthogonally invariant random matrix. Then for all α ∈ A which are not Eulerian, E zα (A(n) ) = 0. We show this at the beginning of our proof in Section B. Both the proof of [CDM24, Theorem 1.1] and our proof of Theorem 4.2 are based on the Weingarten calculus, a combinatorial description of the entrywise moments of Haar-distributed matrices from a matrix group. In Section A, we present an alternative (albeit non-rigorous) derivation of Theorem 4.2 using the Feynman diagram method from physics. Arguably, the combinatorics of the Feynman diagram method is simpler than that of the Weingarten calculus proof.

4.3

Block-structured random matrices

Wigner random matrices and orthogonally invariant random matrices both extend the GOE in different directions, while still satisfying the factorizing strong cactus property. We now consider a third generalization, block matrices, which typically do not satisfy the factorizing property. (n)

n/q×n/q

Fix q ∈ N. For r, c ∈ [q], let Ar,c = Ar,c ∈ Rsym be a sequence of random matrices with Ar,c = Ac,r . The corresponding block matrix model is the symmetric n-by-n matrix whose rows and columns are partitioned into blocks of sizes n/q which has blocks (Ar,c )r,c∈[q] . We let block(i) ∈ [q] denote the block label of i ∈ [n]. The simplest example of a block matrix model is the block GOE model, which has previously been studied in the context of the Generalized AMP algorithm [JM13].10 10

In this paper, we study a slightly more symmetric variant, in which the blocks themselves are symmetric. This modification is made purely for technical reasons, since we work in our other definitions only with symmetric matrices.

Definition 4.4 (Block GOE model). Let q ∈ N and let Σ ∈ Rq×q be a symmetric with nonnegative n/q×n/q entries. For 1 ≤ r ≤ c ≤ q, let Ar,c ∈ Rsym be a symmetric random matrix whose entries on and above the diagonal are independent Gaussians with mean 0 and variance Σ[r, c]/n, and let Ar,c = Ac,r for q ≥ r > c ≥ 1. The block GOE model A ∼ BlockGOE(n, Σ) is the block matrix with blocks (Ar,c )r,c∈[q] . Following the arguments of [Mal20, JP25], one can prove that the block GOE model with fixed parameter Σ satisfies the strong cactus property. Indeed, as in Theorem 4.1, it is still only the doubled trees or cactuses of 2-cycles that have non-zero value in the traffic distribution. However, these values depend non-trivially on Σ, and in general the block GOE model does not satisfy the factorizing strong cactus property.11 Traffic independence. We study block models through the notion of traffic independence. Traffic independence was introduced by Male [Mal20] as a generalization of free independence of matrices. Free independence is a property of the mixed traces of several random matrices (in our notation, these traces are represented by cycle diagrams), whereas traffic independence is a property of all diagrams. Using this concept, below we prove a general result that block-structured matrices have the strong cactus property provided that (i) each of the blocks separately has the strong cactus property, and (ii) those blocks are asymptotically traffic independent. k For a sequence of symmetric matrices (A1 , . . . , Ak ) ∈ (Rn×n sym ) , we generalize the graph polynomials to wα (A1 , . . . , Ak ) and zα (A1 , . . . , Ak ), where α is a multigraph whose edges are additionally colored by A1 , . . . , Ak . The graph polynomial defined by α uses the entries of Ai on each edge whose color is Ai , as in Definition 3.6.

Define a colored component to be a maximal connected subgraph of α whose edges all have the same label Ai . Let CC(α) denote the set of colored components. Define the graph of colored components GCC(α) to be the bipartite graph χ with: V (χ) = CC(α) ∪ {u ∈ V (α) : u belongs to at least two colored components}, E(χ) = {(C, u) : u belongs to the colored component C}. (n)

(n)

k Definition 4.5 (Traffic independence). Let (A1 , . . . , Ak ) = (A1 , . . . , Ak ) ∈ (Rn×n sym ) be sequences of symmetric random matrices, with respective limiting traffic distributions D1 , . . . , Dk . We say that A1 , . . . , Ak are asymptotically traffic independent if, for all connected undirected multigraphs α with edges labeled by A1 , . . . , Ak ,

  

1 lim E zα (A1 , . . . , Ak ) = C∈CC(α) n→∞ n A1 ,...,Ak   0

Di(C) (C)

if GCC(α) is a tree otherwise

Here, i(C) denotes the matrix label associated with the colored component C. Next, we prove that traffic independence of the blocks preserves the strong cactus property: (n)

n/q×n/q

Proposition 4.6. Let q ∈ N. For r, c ∈ [q], let Ar,c = Ar,c ∈ Rsym be a sequence of symmetric random matrices such that Ar,c = Ac,r . Assume that each Ar,c has a limiting traffic distribution 11

If the row sums of Σ are constant, yielding what is sometimes called a generalized Wigner matrix, then up to rescaling the traffic distribution is again that of the GOE and the factorizing property does hold.

that satisfies the strong cactus property and (Ar,c )1≤r≤c≤q are asymptotically traffic independent. Then, the block matrix A ∈ Rn×n sym with blocks (Ar,c )r,c∈[q] also has a limiting traffic distribution that satisfies the strong cactus property. Proof. Let α ∈ A. In the graph polynomial zα (A) we partition the sum based on the block of each vertex: X X Y 1 1 zα (A) = A [i(u), i(v)] . n n χ:V (α)→[q] i:V (α)→[ n ] uv∈E(α) χ(u),χ(v) q

We can interpret the inner summation as a generalized graph polynomial whose edges are labeled by the matrices Ar,c . Call this diagram αχ and write: X 1 1 zα (A) = zαχ ((Ar,c )r,c∈[q] ) . n n χ:V (α)→[q]

Taking the expectation and the limit n → ∞, by traffic independence, all limits exist (so the block matrix has a limiting traffic distribution), and the nonzero terms on the right-hand side are those for which GCC(αχ ) is a tree. By the strong cactus property for each Arc , each colored component must be a cactus. Therefore, any nonzero α is formed by gluing several cactuses along a tree, which forms a bigger cactus. Finally, traffic independence is shown in [Mal20] to hold quite generally for independent random matrices Ai , each of which has a permutation-invariant distribution. Theorem 4.7 ([Mal20, Theorem 1.8]). Let A1 , . . . , Ak ∈ Rn×n sym be independent random matrices such that for each i ∈ [k], (i) The law of Ai ∈ Rn×n sym is Sn -invariant (i.e., invariant under the simultaneous action of Sn on the rows and columns of Ai ). (ii) The limiting traffic distribution of Ai exists. (iii) The traffic distribution concentrates for Ai (Definition 3.15). Then A1 , . . . , Ak are asymptotically traffic independent. Together with Proposition 4.6, Theorem 4.7 implies that block-structured matrices with independent blocks, each satisfying the strong cactus property and Conditions (i), (ii), (iii) also satisfy the strong cactus property (such as the block GOE matrix). We note that Condition (i) can be ensured by applying an independent random permutation to the rows and columns of each Ai . Condition (iii) is proven for orthogonally invariant random matrices in Lemma B.7.

Universality for Deterministic Matrices

Recall the definition of puncturing (Definition 2.1) and of the r-ROM (Definition 2.5). Our main theorem in this section is:

Theorem 5.1. Let H = H (n) ∈ Rn×n sym be a sequence of symmetric orthogonal matrices such that 1

max |H[i, j]| ≤ n− 2 +o(1) .

(11)

1≤i≤j≤n

Then, the limiting traffic distribution of the puncturing of H exists and equals that of the r-ROM. Theorem 5.1 directly applies to H being the sequence of Walsh-Hadamard matrices, discrete sine transform matrices, or discrete cosine transform matrices. Theorem 5.1 follows from the more general Theorem 5.3 below, which applies to symmetric matrices that are not necessarily orthogonal, but have a limiting diagonal distribution and satisfy a generalized delocalization assumption. (n) > 0. We introduce the assumptions: Assumption 5.2. Let H = H (n) ∈ Rn×n sym and ε = ε

∥H∥ ≤ 1, max |Wα (H)[i, j]| ≤ ε

1≤i<j≤n

1 √ ∥Πwσ (H)∥2 ≤ ε n

(12) for each open cactus α (Definition 5.4),

(13)

for all σ ∈ C1 ,

(14)

where Π = Π(n) = I − n1 11⊤ denotes the projection orthogonal to the all-ones direction. For example, one of the constraints of Eq. (13) is that |H k [i, j]| ≤ ε uniformly for all k, n ∈ N and distinct i, j ∈ [n] (a bound which is uniform in n, i, j but may depend on k would also be sufficient, but we omit this for simplicity). (n) > 0. Theorem 5.3 (Universality). Let H = H (n) ∈ Rn×n sym , A be the puncturing of H, and ε

1. If H satisfies Eqs. (12) and (13), then for all α ∈ E \ C, 1 1 |zα (H)| ≤ Oα ε(n) + √ n n

1 1 |zα (A)| ≤ Oα ε(n) + √ n n

and

In particular, if ε(n) = o(1), then both H and A satisfy the weak cactus property. 2. If H satisfies Eqs. (12) to (14), then for all α ∈ A \ E, √ Oα (1) 1 1 |wα (A)| ≤ √ · 1 + ε(n) n . n n 1

In particular, if ε(n) = n− 2 +o(1) , then the right-hand side is n− 2 +oα (1) . 1

Hence, if H satisfies Eqs. (12) to (14) with ε(n) = n− 2 +o(1) , and the diagonal distribution of H exists, then the traffic distribution of A exists and is determined by the diagonal distribution of H. We emphasized in the statement that all constants in the O notations depend on α. We will drop this dependency in the rest of the section. Comparison with prior work. In [WZF22, Theorem 2.8], the authors assume (i) delocalization of open cactuses (Eq. (13)) and (ii) the existence of a limiting diagonal distribution. They show that, after conjugation by a randomly signed permutation matrix, the resulting “semi-random” matrix 28

lies in the same universality class (in the sense of AMP dynamics) as an orthogonally invariant matrix with the same diagonal distribution. Theorem 5.3 shows that the same conclusion holds for deterministic matrices, if we replace random conjugation with puncturing. The universality result of [WZF22] can also be extended in a black-box way to deterministic matrices, but only for GFOM with odd nonlinearities [DLS23, ZWF24]. This assumption lets one only consider the limiting traffic distribution evaluated on Eulerian diagrams. Under the same assumption, our proof would also significantly simplify. Indeed, in Theorem 5.1, the number of monomials appearing in wα (H) is O(n|V (α)| ), and each term has magnitude maxi,j∈[n] |H[i, j]||E(α)| ≤ n−|E(α)|/2+o(1) , giving the upper bound |wα (H)| ≤ no(1) if α has minimum degree 4. It only remains to incorporate paths of degree-2 vertices, which simply compute H k ∈ {I, H} for some k ≥ 1.

5.1

Calculation of cactus diagrams and diagonal distribution

To apply Theorem 5.3, one needs to compute the diagonal distribution of H and small strengthenings of it in order to verify Assumption 5.2. Notice that the only diagrams involved in the assumptions are cactuses, so this is a much simpler task than calculating the entire traffic distribution. In this subsection, we do this calculation directly to prove Theorem 5.1 assuming Theorem 5.3. Let H be a delocalized orthogonal matrix satisfying the assumption of Theorem 5.1. Note that it satisfies H 2 = I. Hence, Eq. (12) is automatic. Next, we define the notion of open cactus appearing in Eq. (13). An open cactus is a matrix diagram with two roots such that merging the roots yields a cactus. Definition 5.4. An open cactus is a graph obtained from a simple path by attaching vertex-disjoint cactuses to each vertex of the path. Formally, α = (V (α), E(α)) is an open cactus if there exist k ≥ 2, vertex-disjoint cactuses β1 , . . . , βk , and distinct vertices u1 ∈ V (β1 ), . . . , uk ∈ V (βk ) with V (α) =

k [

V (βi ) ,

E(α) = {{ui , ui+1 } : i ∈ {1, . . . , k − 1}} ∪

i=1

k [

E(βi ) .

i=1

We call (u1 , uk ) the endpoints of α, and (u1 , . . . , uk ) the base path of α. Unless specified otherwise, we will view an open cactus α ∈ A2 as a matrix diagram rooted at its two ordered endpoints. In general, if α is a matrix diagram and α′ is the scalar diagram formed by merging the roots of α, then Tr(Wα (A)) = wα′ (A). For an open cactus α, this α′ is a cactus, and so wα′ (A) is one of the quantities whose limit is included in the diagonal distribution of A; further, all values of the diagonal distribution can be obtained in this way from the diagonal entries of open cactus matrices. From this perspective, Eq. (13) is a natural counterpart to the diagonal distribution since it concerns all of the off-diagonal entries of the open cactus matrices. We compute the open cactus matrices for H in the following lemma. Lemma 5.5. Let σ be an open cactus and let H satisfy Eq. (11). If all cycles in all of the hanging cactuses have even length, then Wσ (H) = I if the base path has even length and Wσ (H) = H if 1 the base path has odd length. Otherwise, ∥Wσ (H)∥ ≤ n− 2 +o(1) . Proof. First, the leaf 2-vertex-connected components of σ consisting of cycles of even length can be iteratively removed without changing the value of Wσ (H). This is because a hanging cycle of even length k contributes diag(H k ) = diag(I) = 1 in the definition of Wσ . Therefore, if all cycles in all 29

hanging cactuses have even length, then Wσ (H) = H ℓ ∈ {I, H} where ℓ is the length of the base path. In the remaining case where σ has an odd cycle, we use induction. Let β1 , . . . , βk be the hanging cactuses of σ. We convert each βi into an open cactus diagram βi′ by splitting the vertex at which βi meets σ. With this notation, we have the matrix factorization: Wσ (H) = diag(Wβ1′ (H))H diag(Wβ2′ (H))H . . . Hdiag(Wβk′ (H)) . The odd cycle in σ has either become an odd-length base path in some βi′ or it continues to be an odd cycle in some βi′ . In the second case, by sub-multiplicativity of the spectral norm, 1

∥Wσ (H)∥ ≤ diag(Wβi′ (H)) ≤ Wβi′ (H) ≤ n− 2 +o(1) with the last inequality by induction. In the first case, we have Wβi′ (H) = H. Then 1

diag(Wβi′ (H)) = ∥diag(H)∥ ≤ n− 2 +o(1) by the delocalization assumption, and this case is also complete. We use the lemma to complete the proof of Theorem 5.1. Proof of Theorem 5.1 from Theorem 5.3. Eq. (12) holds automatically for H a symmetric orthogonal matrix. Verifying Eq. (13), Lemma 5.5 implies that the off-diagonal entries of all open cactus matrices satisfy 1 max |Wσ (H)[i, j]| ≤ ∥Wσ (H)∥ ≤ n− 2 +o(1) 1≤i<j≤n

when σ has an odd cycle, and the remaining cases Wσ (H) = H or Wσ = I are easily checked. Next, each vector cactus diagram σ ∈ C1 satisfies wσ (H) = diag(Wσ′ (H)) where σ ′ is an open cactus obtained by splitting the root of σ. By Lemma 5.5 the diagonal of an open cactus matrix is either 1 (in which case Eq. (14) is satisfied with ε = 0) or it satisfies 1 1 √ ∥diag(Wσ′ (H))∥2 ≤ ∥diag(Wσ′ (H))∥∞ ≤ n− 2 +o(1) , n 1

in which case Eq. (14) is satisfied with ε = n− 2 +o(1) . The diagonal distribution is computed by averaging the diagonal entries of open cactus matrices: 

n 1 1 1X lim wσ (H) = lim Wσ′ [i, i] = n→∞ n n→∞ n 0 i=1

if all cycles in σ have even length otherwise

where on the left-hand side, we convert σ ∈ C0 to an open cactus diagram σ ′ by rooting it arbitrarily and splitting the root. The right-hand side is by Lemma 5.5. That is, the diagonal distribution of H is just the indicator function that all cycles of the cactus are even. Thus, we showed that Eqs. (12) to (14) hold and the diagonal distribution converges to the same fixed limit for any orthogonal matrix with delocalized entries. By Theorem 5.3, the traffic distribution of such matrices exists and is always the same. 30

Finally, we show that the r-ROM is also in this class, by showing that, after conditioning on a suitable high-probability event, the above argument applies to an r-ROM matrix as well. Let HROM = QDQ⊤ , where Q is Haar-distributed and D is diagonal with i.i.d. ±1 entries, independent of Q. Claim 5.6. There exists c > 0 such that for any t > 0, 1

max |HROM [i, j]| ≤ t2 n− 2

(15)

i,j∈[n] 2

holds with probability at least 1 − n2 e−ct . Proof. Since every entry of Q is O(n−1/2 )-subgaussian, by a union bound 1

max |Q[i, j]| ≤ tn− 2

i,j∈[n] 2

holds with probability at least 1 − n2 e−Ω(t ) . Next, we have HROM [i, j] = nk=1 D[k, k]Q[i, k]Q[j, k], which, conditioned on Q, is a sum of independent random variables. By Hoeffding’s bound, any fixed entry of HROM is O(σ)-subgaussian with parameter σ 2 :=

n X

Q[i, k]2 Q[j, k]2 ≤ max Q[i, j]2 , i,j∈[n]

k=1

since every row of Q has ℓ2 -norm 1. The conclusion follows from a union bound over all entries. Fix α ∈ A. Let En denote the event Eq. (15), with t = no(1) . By the law of total expectation, we decompose 1 1 1 E wα (ΠHROM Π) = E[wα (ΠHROM Π) | En ] Pr(En ) + E[wα (ΠHROM Π) | Enc ] Pr(Enc ) . n n n The left-hand side converges to the traffic distribution of the r-ROM evaluated at α. Moreover, since ∥ΠHROM Π∥ ≤ 1, we may crudely bound the second term by 1 E [wα (ΠHROM Π) | Enc ] · Pr(Enc ) ≤ n|V (α)|−1 Pr(Enc ) −→ 0 . n→∞ n Since Pr(En ) −→ 1, we deduce that n→∞

lim

n→∞ n

E [wα (ΠHROM Π) | En ] = lim

n→∞ n

E wα (ΠHROM Π) .

Finally, on the event En , the matrix HROM satisfies the assumptions of Theorem 5.1. Consequently, the traffic distribution of punctured delocalized orthogonal matrices coincides with that of the r-ROM, as desired. As a consequence of the above argument, the traffic distribution of the r-ROM is specified implicitly as the solution to the following equations: 1. For every α ∈ A \ E, n1 E wα (A) −→ 0. n→∞

2. For every α ∈ E \ C, n1 E zα (A) −→ 0. n→∞

3. For every σ ∈ C, n1 E wσ (A) −→ 1 if all cycles of σ are even and 0 otherwise. n→∞

These equations determine a unique traffic distribution by Lemma 3.14. It is possible to give an explicit but much more complicated description using the Weingarten calculus, which we do in Section B.6. However, the above characterization is arguably the conceptually clearer one, and we emphasize that it involves both the w- and z-bases. We note also as a point of reference that the last part, the limiting values of cactuses in the w-basis, are the same as those for the (unpunctured) ROM, as follows from combining Claim 3.11 with Lemma 3.12, and corresponds simply to the moments of the Rademacher distribution being 1 for moments of even order and 0 for ones of odd order.

5.2

The fundamental theorem of graph polynomials

The main proof of Theorem 5.3 throughout the rest of the section relies on the “fundamental theorem of graph polynomials” of Bai and Silverstein [BS10]. This result can be used to easily bound 2-edge-connected graph polynomials expressed in the w-basis, which is one reason that it is convenient to restrict to such diagrams in our definition of the weak cactus property. The proof of the fundamental theorem uses a spectral bound on tensor powers of A; see [MS12] for another related result. Theorem 5.7 ([BS10, Theorems A.31 and A.32]). For every n ≥ 1, α ∈ E ∪ E1 ∪ E2 and collection of n × n symmetric matrices A = (Ae )e∈E(α) , Y 1 |wα (A)| ≤ ∥Ae ∥ n e∈E(α)

∥wα (A)∥∞ ≤

if α ∈ E,

∥Ae ∥

if α ∈ E1 ,

∥Ae ∥

if α ∈ E2 .

e∈E(α)

∥Wα (A)∥ ≤

Y e∈E(α)

The result of [BS10] only covers scalar and matrix diagrams, but we provide a quick reduction of the vector case to the scalar case. Proof of vector case of Theorem 5.7. For all q ≥ 1, we can diagrammatically express ∥wα (A)∥2q 2q as the diagram formed by merging 2q copies of α at the root, and then forgetting the identity of the root to obtain a scalar diagram. Let α2q = α⊕2q denote this diagram. The graph α2q remains 2-edge-connected, therefore by the scalar case of the result we have: 2q

  ∥wα (A)∥2q 2q = wα2q (A) ≤ n ·

∥Ae ∥

e∈E(α)

Taking q → ∞ with n fixed, we obtain ∥wα (A)∥∞ ≤

e∈E(α) ∥Ae ∥ .

We will apply the fundamental theorem by decomposing a general graph into its 2-edge-connected components, which are joined together by a tree of bridge edges. Decomposing diagrams into their 2-edge-connected components is also a fundamental idea in physics, where a 2-edge-connected Feynman diagram is called a “1-particle-irreducible diagram”.

5.3

Main structural lemma: Open cactus decomposition

To prove the weak cactus property of Theorem 5.3, we begin by observing that any 2-edge-connected non-cactus graph contains three edge-disjoint paths between some pair of vertices. How can we quantify that such a graph is a cactus plus excess edges? We answer this question by introducing the open cactus decomposition. Our main structural result is that one can identify an “extra” open cactus subgraph inside any 2-edge-connected graph which is not a cactus, in the sense that the subgraph can be removed without spoiling 2-edge-connectedness. Proposition 5.8. For any α ∈ E1 \ C1 , there exist distinct s, t ∈ V (α) and an induced subgraph β of α such that 1. β is an open cactus with endpoints {s, t}. 2. α [V (α) \ (V (β) \ {s, t})] is 2-edge-connected. 3. root(α) ∈ / V (β) \ {s, t}.

Figure 3: Example for Proposition 5.8 of a 2-edge-connected graph which is not a cactus. If the open cactus in red is removed, the graph remains 2-edge-connected. To prove Proposition 5.8, we will consider the last ear in an ear decomposition of α. We prove a small variant of the classical ear decomposition (see [Rob39] or [BM08, §5.3]) which lets us exclude a specified vertex from the internal vertices of the last ear. Lemma 5.9. Let α ∈ E1 be 2-edge-connected with at least 2 vertices. There exists a path π = (u1 , . . . , uk ) in α with k ≥ 2 such that: 1. Each internal vertex u2 , . . . , uk−1 has degree 2 in α. 2. Each internal vertex u2 , . . . , uk−1 satisfies ui ̸= root(α). 3. u1 , . . . , uk are pairwise distinct, except possibly u1 = uk . 4. Removing internal vertices and edges of π from α leaves α 2-edge-connected. Proof of Lemma 5.9. Consider the following sequence (αt )t≥0 of 2-edge connected subgraphs of α: 33

1. Start from α0 being any cycle of α containing root(α). 2. Let t ≥ 0. If αt spans all vertices of α, then stop. / V (αt ). Since α is 3. Otherwise, there exists {u1 , u2 } ∈ E(α) such that u1 ∈ V (αt ) and u2 ∈ 2-edge-connected, there exists a simple path (u2 , . . . , uk ) in α \ {{u1 , u2 }} such that ui ∈ / V (αt ) for all 2 ≤ i ≤ k − 1, and uk ∈ V (αt ). Set αt+1 = (V (αt ) ∪ {u2 , . . . , uk−1 }, E(αt ) ∪ {{ui , ui+1 } : 1 ≤ i < k}) . For any t ≥ 0, αt is 2-edge-connected. Therefore, if at the end of the algorithm V (αt ) = V (α) but E(αt ) ̸= E(α), then any edge in E(α) \ E(αt ) is a length-1 path that satisfies the conclusion of the lemma. Otherwise, this means that α is obtained from αt−1 (which is 2-edge-connected) by adding a path of internal degree-2 vertices in α which must all be distinct from root(α) ∈ V (α0 ) ⊆ V (αt−1 ). This concludes the proof. Proof of Proposition 5.8. Starting with the graph α, consider the following procedure: 1. Delete all self-loops in α. 2. If no leaf 2-vertex-connected component (i.e., a 2-vertex-connected component meeting the rest of the graph at a single articulation point) consists of a single cycle, then stop. 3. Otherwise, choose an arbitrary such component. Let v be the articulation point connecting this component to the rest of graph. Delete all edges of this component from the graph. 4. Delete newly isolated vertices; exactly one vertex of the component remains, namely v. Since α∈ / C1 , the procedure does not delete the entire graph. 5. If the root was removed in Step 4, set v as the new root of the diagram. 6. Return to Step 1. Call β ∈ A1 the resulting rooted graph. Note that β is still 2-edge-connected, so by Lemma 5.9, we can find a path π = (u1 , . . . , uk ) in β with internal degree-2 vertices. π cannot be a cycle because of our initial step of removing cyclic 2-vertex-connected components. Therefore, π is a simple path and the root of β is not an internal vertex of π. Observation 5.10. For 2 ≤ i < k, let σi be the connected component of ui in α \ E(π). Then α′ := π ∪ σ2 ∪ . . . ∪ σk−1 is an open cactus in α with endpoints u1 , uk . Moreover, root(α) is not an internal vertex of the open cactus. Proof. π is a simple path in β, and adding back loops and cyclic 2-vertex-connected components we removed from α, we obtain an open cactus. The recursive pruning procedure we used to transfer the root ensures that root(α) is not in any of the cyclic 2-vertex-connected components that are added to π. Observation 5.11. α [V (α) \ (V (α′ ) \ {u1 , uk })] is 2-edge-connected. Proof. By Lemma 5.9, β [V (β) \ {u2 , . . . , uk−1 }] is 2-edge-connected. Adding 2-vertex-connected cyclic components to this graph preserves 2-edge-connectivity. 34

Observation 5.10 and Observation 5.11 conclude the proof of Proposition 5.8.

5.4

The effect of puncturing

The main result of this subsection is: n Proposition 5.12. Let H ∈ Rn×n sym such that ∥H∥ ≤ 1 and u ∈ R be a unit vector. Denote by ⊤ ⊤ A = (I − uu )H(I − uu ). Then for any open cactus α ∈ A2 ,

∥Wα (A) − Wα (H)∥F ≤ |E(α)| · ∥A − H∥F ≤ 3|E(α)| . We deduce in the following that puncturing does not change the diagonal distribution. In particular, matrices such as the ROM and the r-ROM have the same diagonal distribution. Corollary 5.13. Let H and A be as in Proposition 5.12. Then for any σ ∈ C1 ∥wσ (H) − wσ (A)∥2 ≤ O(1) , and for any σ ∈ C,

1 1 |wσ (H) − wσ (A)| ≤ O √ n n

Proof of Corollary 5.13 from Proposition 5.12. Except for the case where σ ∈ C1 has one vertex (in which case the statement holds because the diagonal entries are bounded), root(σ) has degree ≥ 2. Create two copies r1 , r2 of root(σ) and re-assign the edges incident to root(σ) to r1 or r2 in such a way that r1 and r2 have degree at least 1. The resulting graph is an open cactus α with endpoints r1 and r2 such that merging these endpoints yields back σ. Hence, ∥wσ (H) − wσ (A)∥2 = ∥diag(Wα (H)) − diag(Wα (A))∥F ≤ O(1) . The second statement then follows from Cauchy-Schwarz: |wσ (H) − wσ (A)| = |⟨1, wσ (H) − wσ (A)⟩| ≤

√

√ n · ∥wσ (H) − wσ (A)∥2 ≤ O( n) .

This concludes the proof. However, H and its punctured version A may not have the same traffic distribution, even on scalar open cactuses. Thus, the diagonal distribution (i.e., the values of cactus diagrams) is not sensitive to the behavior of H in any single direction u, while some diagrams in the traffic distribution are sensitive to the behavior in the 1 direction. Example 5.14 (Puncturing of the Walsh-Hadamard matrix). Let H (n) be the normalized WalshHadamard matrices (Definition 2.3). Then for the 2-path diagram α (which is an open cactus), 1 1 (wα (H) − wα (A)) = ⟨1, (H 2 − A2 )1⟩ −→ 1 . n→∞ n n

This does not contradict Proposition 5.12: E = Wα (H) − Wα (A) indeed satisfies n X

E[i, j]2 ≤ O(1)

and

n X

E[i, j] = Ω(n) .

i,j=1

In general, as the following example demonstrates, the off-diagonal structure of the error matrix E = Wα (H) − Wα (A) in Proposition 5.12 may be intricate. In the following example, E has entries of magnitude Ω(1), even though its Frobenius norm remains bounded. Example 5.15 (Puncturing of the DST matrix). Let H (n) be the discrete sine transform matrices (Definition 2.4). Then for any fixed odd i ≥ 1, the normalized sum of the ith row of H (n) is √ Z 1 n √ 2 2 1 X √ sin(πit) dt −→ H[i, j] = ( 2 + o(1)) . n→∞ iπ n j=1 0 Consider the 2-path diagram α. While the off-diagonal entries of Wα (H) = H 2 vanish (since H is a symmetric orthogonal matrix), on the other hand, for any fixed distinct odd numbers i, j ≥ 1, Wα (A)[i, j] = (A2 )[i, j] −→ − n→∞

8 , ijπ 2

which is Ω(1) for constant i ̸= j. The proof of Proposition 5.12 relies on expanding A in terms of uu⊤ and H. All rank-1 terms can be neglected thanks to the following lemma: Lemma 5.16. Let α be an open cactus, e∗ ∈ E(α), and A = (Ae )e∈E(α) be a collection of matrices such that ∥Ae ∥ ≤ 1 for all e ∈ E(α) \ {e∗ }. Then, ∥Wα (A)∥F ≤ ∥Ae∗ ∥F . Proof. We first run a pruning procedure that iteratively removes parts of α not containing e∗ , without decreasing the Frobenius norm of Wα (A) during the procedure. To this end, we use repeatedly the standard inequalities: Claim 5.17. ∥M1 M2 ∥F ≤ ∥M1 ∥F ∥M2 ∥. Claim 5.18. ∥M1 ⊙ M2 ∥F ≤ ∥M1 ∥F · max1≤i,j≤n |M2 [i, j]| ≤ ∥M1 ∥F ∥M2 ∥, where ⊙ denotes entrywise or Hadamard product. Initially, let uL be one of the endpoints of α. 1. If e∗ belongs to a cactus hanging from uL , then stop. 2. Otherwise, remove the cactus hanging from uL from the diagram. Using Claim 5.18 and Theorem 5.7 (the spectral norm of the cactus matrix diagram with a double root at uL is at most 1), this does not decrease the Frobenius norm. 3. At this point, uL must have degree equal to 1 in the current graph. If e is the edge adjacent to uL , then stop.

4. Otherwise, remove the edge adjacent to uL . By Claim 5.17 and the assumption, this does not decrease the Frobenius norm. Set uL to be the vertex that was adjacent to uL , and go back to the first step. Then, apply the symmetric procedure from the other endpoint uR of α. At this point, there are two cases. If uL ̸= uR , then the resulting graph must consist of the single edge e∗ = {uL , uR }, so we get the desired upper bound on the Frobenius norm. Therefore, we assume from now on that uL = uR . The resulting graph must be a cactus rooted at uL = uR , and e∗ is one of the edges of this cactus. If there are several cycles incident to uL , we use again Claim 5.18 and Theorem 5.7 to remove all such cycles not containing e∗ without decreasing the Frobenius norm. Finally, we bound the Frobenius norm of the diagonal cactus matrix rooted at uL by the Frobenius norm of an open cactus obtained by creating two copies of the root and turning the unique cycle hanging at uL into a simple path between these two copies (we used a similar procedure in Corollary 5.13). We claim that this open cactus has strictly less edges than the one we started with before running the pruning procedure. Indeed, the base path had at least one edge, which was removed during the pruning stage when uL = uR at the end. We conclude by induction on the number of edges of the open cactus. Proof of Proposition 5.12. We replace iteratively H by A in the graph polynomial Wα (H): let e1 , . . . , e|E(α)| be the edges of α, and write |E(α)|

Wα (A) − Wα (H) =

Wα (Ai ) ,

i=1

where Ai [ej ] = H if j < i, Ai [ej ] = A if j > i, and Ai [ei ] = A − H. For each i ∈ [|E(α)|], we apply Lemma 5.16 with e∗ = ei . We have ∥A∥ ≤ 1 and ∥H∥ ≤ 1 so the assumptions of the lemma are satisfied, and we deduce ∥Wα (Ai )∥F ≤ ∥A − H∥F , and by the triangle inequality ∥Wα (A) − Wα (H)∥F ≤ |E(α)| · ∥A − H∥F . Finally, we have

A − H = ⟨u, Hu⟩uu⊤ − (Huu⊤ + uu⊤ H) .

Since ∥H∥ ≤ 1 and u is a unit vector, we have |⟨u, Hu⟩| ≤ 1 and ∥Hu∥2 ≤ 1, so ∥A−H∥F ≤ 3.

5.5

Support of the z-basis

Let H = H (n) be a family of matrices satisfying Eqs. (12) and (13) and A = A(n) be their puncturing. The main result of this subsection is that A and H satisfy the weak cactus property, that is, their traffic distribution in the z-basis is supported on cactuses and graphs with bridges. Proposition 5.19. For any α ∈ E \ C, 1 1 |zα (H)| ≤ O ε + √ n n

1 1 |zα (A)| ≤ O ε + √ n n

and

The fundamental theorem of graph polynomials can be used to show that these quantities are O(1) (after converting between the z and w-bases). The idea of Proposition 5.19 is to isolate an open cactus in α by Proposition 5.8 and apply Assumption 5.2 to gain an additional ε factor. We emphasize that analogous bounds in the w-basis are false in general; summation over some distinct indices is necessary to prove Proposition 5.19. We prove that, using the notation in Section 3.2: Lemma 5.20. Let α ∈ E1 \ C1 and let s, t be the endpoints of an open cactus in α satisfying the guarantees of Proposition 5.8. Then 1 1 √ ∥wαs̸=t (A)∥2 ≤ O ε + √ n n

1 1 √ ∥wαs̸=t (H)∥2 ≤ O ε + √ n n

and

(16)

The constraint s = ̸ t in Eq. (16) ensures that we only use off-diagonal entries of the open cactuses in the graph polynomial. These are the only entries assumed to be small in Assumption 5.2 (and indeed, the diagonal entries of Wα (H) can be large, for example, in the 2-path diagram). Proof of Proposition 5.19 from Lemma 5.20. Let M ∈ {A, H} and s, t be two distinct vertices of α to be fixed later. Using Möbius inversion (Lemma 3.9) recursively, we can expand zα (M ) = cα wαs̸=t (M ) +

cβ zβ (M ) ,

β≺α

for some constant coefficients cβ ∈ R. Since all β ≺ α are 2-edge-connected by Lemma 3.13 and have strictly less vertices than α, by induction on the number of vertices of α, it suffices to prove: 1 1 s̸=t |wα (M )| ≤ O ε + √ n n

(17)

But Eq. (17) follows from Lemma 5.20: pick s, t to be the endpoints of an open cactus decomposition provided by Proposition 5.8, so that by Cauchy-Schwarz 1 s̸=t 1 1 1 |w (M )| = ⟨wαs̸=t (M ), 1⟩ ≤ √ ∥wαs̸=t (M )∥2 ≤ O ε + √ n α n n n

which concludes the proof. We now move to the proof of Lemma 5.20. A useful concept will be the following graphical interpretation of squaring the polynomial expressed by a diagram: Definition 5.21 (Lift). Let α ∈ A and T ⊆ V (α). Let S1 and S2 be two new disjoint sets of size |V (α)| − |T | (also disjoint from V (α)). For i ∈ {1, 2}, let pi be a bijection between V (α) \ T and Si , which is extended to V (α) by pi (u) = u for all u ∈ T . The lift of α with respect to T is the graph LiftT (α) with V (LiftT (α)) = T ∪ S1 ∪ S2 ,

E(LiftT (α)) = {{pi (u), pi (v)} : i ∈ {1, 2}, {u, v} ∈ E(α)} .

Claim 5.22. Let α ∈ A2 with roots (s, t), and T ⊆ V (α) be such that {s, t} ⊆ T . Then for any

M ∈ Rn×n sym , 2



WLiftT (α) (M )[i, j] =

 φ : T →[n] φ(s)=i,φ(t)=j

M [φ(u), φ(v)] .

φ : V (α)\T →[n] {u,v}∈E(α)

Lemma 5.23. Let α ∈ E, let π be a connected subgraph of α, let α1 be any connected component of α \ E(π), and let α2 the graph spanned by E(α) \ E(α1 ). Then for all j ∈ {1, 2}, LiftV (α1 )∩V (α2 ) (αj ) is 2-edge-connected. Proof. First, α1 is connected by definition. Since α is connected, every connected component in (V (α), E(α) \ E(π)) must be connected to π. Together with the fact that π itself is connected, we get that α2 is connected. In particular, the lifts of α1 and α2 are connected. Fix j ∈ {1, 2} and an edge e′ in the lift of αj . We need to show that e′ belongs to at least one simple cycle in the lift of αj . There exist i ∈ {1, 2} and e = {x, y} ∈ E(αj ) such that e′ = {pi (x), pi (y)} (where p1 , p2 are the lift maps from Definition 5.21). Since α is 2-edge-connected, e belongs to a simple cycle in α. Consider the longest subpath of this cycle containing e and consisting only of vertices in V (αj ). If this subpath is the entire cycle, then we have found a cycle containing e in αj , and so a cycle containing e′ in its lift. Otherwise, the endpoints of this path must be in V (α1 ) ∩ V (α2 ). The images of this path through the lift maps p1 and p2 are disjoint, so their union forms a cycle in the lift of αj containing e′ . Lemma 5.24. Let α ∈ E2 have two distinct roots. Let β be a leaf 2-vertex-connected component of α (i.e., removing internal vertices of β leaves α connected) that does not contain the roots of α. We view β ∈ E1 as a vector diagram rooted at the articulation point connecting β to the rest of α. For any distinct s′ , t′ ∈ V (β) and M ∈ Rn×n sym such that ∥M ∥ ≤ 1, n X

′

Wαs ̸=t (M )[i, j] ≤

√

′

n · ∥wβs ̸=t (M )∥2 .

i,j=1

Proof. Let (s, t) be the roots of α. Since α is 2-edge-connected, there exist two edge-disjoint simple paths between s and t. Let π be one of them. Let α1 be the connected component of s in (V (α), E(α) \ E(π)), and α2 be the graph spanned by E(α) \ E(α1 ) (including only the vertices incident with one of these edges). Finally, let S = V (α1 ) ∩ V (α2 ). Claim 5.25. {s, t} ⊆ S. Proof. On the one hand, E(π) ⊆ E(α2 ) and {s, t} are the endpoints of π, so {s, t} ⊆ V (α2 ). On the other hand, s ∈ V (α1 ) by definition, and there is an s–t path in α \ E(π), so t ∈ V (α1 ). Claim 5.26. For any {u, v} ∈ E(α) with u ∈ V (α1 ) and v ∈ V (α2 ), we have u ∈ S or v ∈ S. Proof. Suppose that v ∈ / V (α1 ). Since u ∈ V (α1 ), u is connected to s by edges of E(α) \ E(π), and since v ∈ / V (α1 ), v is not connected to s by these edges. But, {u, v} ∈ E(α), so it must be that {u, v} ∈ E(π). And, E(π) ⊆ E(α2 ), so u ∈ V (α2 ).

As π is a simple path and β is connected to the rest of α at an articulation vertex, π does not contain any edge of β, so it must be that either E(β) ⊆ E(α1 ) or E(β) ⊆ E(α2 ). Assume without loss of generality that this holds for α1 (the argument will be exactly symmetric for α2 , as we will only use the fact that these subgraphs satisfy the conclusion of Lemma 5.23). In particular, we then have s′ , t′ ∈ V (α1 ). We first use the triangle inequality to push the absolute value inside the sum over labelings of vertices in S: n X

′

Wαs ̸=t (M )[φ(s), φ(t)]

φ(s),φ(t)=1

≤

M [φ(u), φ(v)]

(18)

φ : S→[n] φ : V (α)\S→[n] {u,v}∈E(α) φ(s′ )̸=φ(t′ )

2 Y

M [φ(u), φ(v)]

φ : S→[n] j=1 φ : V (αj )\S→[n] {u,v}∈E(αj ) φ(s′ )̸=φ(t′ ) if j = 1



2  12



2 Y  

≤

j=1 φ : S→[n]

   

     

(19)

M [φ(u), φ(v)]  ,

φ : V (αj )\S→[n] {u,v}∈E(αj ) φ(s′ )̸=φ(t′ ) if j = 1

where we applied Cauchy-Schwarz in the second inequality. Note that Eq. (18) is well-defined by Claim 5.26. By Lemma 5.23 and Claim 5.22, the term for j = 2 in Eq. (19) is a 2-edge-connected graph polynomial, so by Theorem 5.7 and the assumption ∥M ∥ ≤ 1, this term is bounded by 2

 X

M [φ(u), φ(v)] ≤ n .

 φ : S→[n]

φ : V (α2 )\S→[n] {u,v}∈E(α2 )

We now switch to the term j = 1 in Eq. (19). This graph polynomial can be interpreted as 2

 X φ : S→[n]

   

  

′

M [φ(u), φ(v)] = ⟨wβs ̸=t (M ), Wα′ (M )wβs ̸=t (M )⟩ ,

φ : V (α1 )\S→[n] {u,v}∈E(α1 ) φ(s′ )̸=φ(t′ )

where α′ is the lift of α [V (α) \ (V (β) \ {r})] with respect to S (here r denotes the root of β, the articulation vertex connecting β to the rest of α), and we add two roots in α′ at the two copies of r created during the lift operation.

Hence, 2

 X φ : S→[n]

   

′

  

′

M [φ(u), φ(v)] ≤ ∥Wα′ (M )∥ · ∥wβs ̸=t (M )∥22 .

φ : V (α1 )\S→[n] {u,v}∈E(αj ) φ(s′ )̸=φ(t′ )

Note that α′ is 2-edge-connected by Lemma 5.23, so that ∥Wα′ (M )∥ ≤ 1 by Theorem 5.7. Putting everything together, we obtain n X

′

Wαs ̸=t (M )[i, j] ≤

√

′

n · ∥wβs ̸=t (M )∥2 ,

i,j=1

as desired. Proof of Lemma 5.20. Let M ∈ {A, H}. Consider β ∈ A2 defined by: 1. Start from the lift of α with respect to its root. Let p1 and p2 be the lift maps. 2. Delete the edges and internal vertices of the image under p1 of the open cactus in α. 3. Root the resulting graph at p1 (s) and p1 (t). Recall that s and t are the endpoints of the “extra” open cactus in α. Thus, β is, in short, α grafted to its mirror image at the roots, with just one of the copies of that extra open cactus deleted except for its endpoints, and those endpoints made the roots of the matrix diagram β. See Figure 4 for an illustration of this and the rest of the proof. Let σ be the image of the open cactus in α under the lift map p2 , and let s′ and t′ be the images of the endpoints of this open cactus through the lift map p2 . Thus s′ and t′ are the mirror images of the vertices chosen to be the roots of β above. We can then rewrite ∥wαs̸=t (M )∥22 =

n X i,j=1 i̸=j n X

′

Wβs ̸=t (M )[i, j]Wσ (M )[i, j]

(20) n X

Wβs ̸=t (M )[i, j]Wσ (H)[i, j] +

i,j=1 i̸=j

′

Wβs ̸=t (M )[i, j](Wσ (M ) − Wσ (H))[i, j]

i,j=1 i̸=j

≤ max |Wσ (H)[i, j]| 1≤i<j≤n

n X

′

Wβs ̸=t (M )[i, j] + ∥Wσ (M ) − Wσ (H)∥F ∥Wβs ̸=t (M )∥F ,

(21)

i,j=1

using Hölder on the first term and Cauchy-Schwarz on the second. We further bound the first term with Assumption 5.2 and Lemma 5.24: max |Wσ (H)[i, j]| ·

1≤i<j≤n

n X

√ ′ ′ Wβs ̸=t (M )[i, j] ≤ ε n · ∥wαs̸=t (H)∥2 .

i,j=1

For the second term, observe that by Proposition 5.12, we know that the change due to puncturing is small in Frobenius norm, i.e., ∥Wσ (M ) − Wσ (H)∥F ≤ O(1) for M ∈ {A, H}. Moreover, in the ′ ′ other factor, ∥Wβs ̸=t (M )∥2F is nothing but the lift of β with respect to {p1 (s), p1 (t)}. This lift can be interpreted as: ′

′

∥Wβs ̸=t (M )∥2F = ⟨wαs̸=t (M ), Wβ ′ (M )wαs̸=t (M )⟩ ,

(22)

where β ′ is the lift of α [V (α) \ (V (σ) \ {s, t})] with respect to {s, t}. By the guarantees of Proposition 5.8, α [V (α) \ (V (σ) \ {s, t})] is already 2-edge-connected, and therefore so is β ′ . As a result, by Theorem 5.7, ′

′

∥Wσ (M ) − Wσ (H)∥F · ∥Wβs ̸=t (M )∥F ≤ O(1) · ∥Wβ ′ (M )∥ 2 ∥wαs̸=t (M )∥2 ≤ O(1) · ∥wαs̸=t (M )∥2 . We obtain

√ ∥wαs̸=t (M )∥22 ≤ O(1 + ε n) · ∥wαs̸=t (M )∥2 ,

and the result follows after rearranging the inequality.

(a) Example of vector diagram α with extra open cactus in red

(b) Lift of α at the root (Eq. (20))

(d) Lift again. The left and right sides are copies of α while the inner part is 2-edge-connected (Eq. (22)).

Figure 4: Illustration of the main diagrammatic manipulations in the proof of Lemma 5.20.

5.6

Support of the w-basis

In this subsection we prove the second part of Theorem 5.3: Proposition 5.27. Suppose that H satisfies Eqs. (12) to (14). Then for any α ∈ A \ E, √ 1 1 |wα (A)| ≤ √ · (1 + ε n)O(1) . n n

These calculations are simpler than the previous ones, but this is also the point in the proof of Theorem 5.3 where puncturing is essential (note that it was not used at all in the previous section, and indeed those results apply equally well to the original H or the puncturing A). But, without puncturing, the values of graph polynomials that contain bridges can fail to be universal. For instance, when α is the degree-d star, Walsh–Hadamard matrices H = H (n) satisfy wα (H) = Θ(nd/2 ), so the limiting traffic distribution does not even exist when d ≥ 3. As Proposition 5.27 shows, puncturing effectively forces all such diagrams to vanish in the traffic distribution. To prove Proposition 5.27, we will isolate a bridge edge in the graph, and show by induction over the tree of 2-edge-connected components that: Lemma 5.28. For all α ∈ A1 , √ ∥Awα (A)∥2 ≤ (1 + ε n)O(1) . Proof of Proposition 5.27 from Lemma 5.28. Decompose α = α1 ⊔ α2 ⊔ {u, v}, where {u, v} ∈ E(α) is a bridge edge, α1 ∈ E1 is rooted at u, and α2 ∈ A1 is rooted at v. Then, |wα (A)| = |⟨wα1 (A), Awα2 (A)⟩| ≤ ∥wα1 (A)∥2 ∥Awα2 (A)∥2 ≤

√

√ n · (1 + ε n)O(1) ,

using Theorem 5.7 on the first term and Lemma 5.28 on the second. We prove Lemma 5.28 by first treating the cactus special case (Lemma 5.29), then the 2-edgeconnected special case (Lemma 5.30), and then finally the general case by the induction mentioned above. Lemma 5.29. For any α ∈ C1 , √ ∥Awα (A)∥2 ≤ O(1 + ε n) . Proof. We first decompose wα (A) as: wα (A) = (wα (A) − wα (H)) + Πwα (H) +

1 ⟨1, wα (H)⟩1 . n

Since A1 = 0 and ∥A∥ ≤ 1 by assumption, by the triangle inequality we have ∥Awα (A)∥2 ≤ ∥wα (A) − wα (H)∥2 + ∥Πwα (H)∥2 . By Corollary 5.13, the first term is O(1), and by our assumption in Eq. (14), the second term is at √ most ε n. Lemma 5.30. For all α ∈ E1 , √ ∥Awα (A)∥2 ≤ O(1 + ε n) . Proof. We proceed by induction on |V (α)|. For α ∈ C1 (in particular, if α has only one vertex, which is our base case), the claim follows from Lemma 5.29. For α ∈ E1 \ C1 , we apply Proposition 5.8: there is an open cactus σ induced in α such that removing the internal vertices and edges from σ

leaves α rooted and 2-edge-connected. Let {s, t} be the endpoints of σ, and β be the graph obtained from α by merging s and t. Then, we can decompose: wα (A) = wαs̸=t (A) + wβ (A) . On the one hand, β ∈ E1 by Lemma 3.13 and has strictly less vertices than α, so by induction √ ∥Awβ (A)∥2 ≤ O(1 + ε n) . On the other hand, by Lemma 5.20, √ ∥wαs̸=t (A)∥2 ≤ O(1 + ε n) . Putting everything together and using ∥A∥ ≤ 1 and the triangle inequality, we obtain √ ∥Awα (A)∥2 ≤ O(1 + ε n) , which concludes the induction. Lemma 5.31. Let α ∈ E, v ∈ V (α), and S the set of edges adjacent to v in α. Then there exists β ∈ E and v1 , v2 ∈ V (β) such that V (β) = (V (α) \ {v}) ∪ {v1 , v2 } , E(β) = (E(α) \ S) ∪ ϕ(S) ∪ {{v1 , v2 }} , where ϕ(e) ∈ {{v1 , u}, {v2 , u}} for all e = {v, u} ∈ S. Proof. We use the ear decomposition construction from the proof of Lemma 5.9. Consider the step of the ear decomposition which adds v. During this step, v is a new interior vertex of a path or cycle added to α. We define β by splitting v into two vertices v1 , v2 with a new edge between them. When other ears attach to v in α, we can attach them to either v1 or v2 in β. This process yields an ear decomposition for β, hence β is also 2-edge-connected. Proof of Lemma 5.28. We proceed by induction on the number of 2-edge-connected components in α. If α is 2-edge-connected, then the result follows by Lemma 5.30. We assume from now on that α is not 2-edge-connected. Let C be the 2-edge-connected component of the root of α. Let β1 , . . . , βk (k ≥ 1) be the connected components disjoint from β in the graph obtained after removing E(C) and all bridges incident to C. We root βi at the (unique) vertex of V (βi ) adjacent to β. We also consider u1 ∈ V (α), the unique vertex in V (C) that is adjacent to V (β1 ). Let β ∈ A2 be the graph obtained from α by adding a second root at u1 , and deleting V (β1 ), E(β1 ), and the bridge between u1 and β1 . Then, for i = 2, . . . , k, we iteratively apply the graph transformation from Lemma 5.31, label the new edge e = {v1 , v2 } by Ae = diag(Awβi (A)), and transfer the old labels for all other edges. In this way, we obtain a 2-edge-connected graph β ′ ∈ E2 and a family of matrices A = (Ae )e∈E(β ′ ) such that Wβ ′ (A) = Wβ (A) .

All involved matrices Ae ∈ A are either A or of the form diag(Awβi (A)) for some i ∈ {2, . . . , k}, √ so they satisfy ∥Ae ∥ ≤ (1 + ε n)O(1) by induction. Next, applying Theorem 5.7, we get √ ∥Wβ (A)∥ = ∥Wβ ′ (A)∥ ≤ (1 + ε n)O(1) . As a result, √ ∥Awα (A)∥2 = ∥AWβ (A)Awβ1 (A)∥2 ≤ ∥A∥ · ∥Wβ (A)∥ · ∥Awβ1 (A)∥2 ≤ (1 + ε n)O(1) , using again induction on ∥Awβ1 (A)∥2 . This concludes the induction.

5.7

Putting everything together: Proof of Theorem 5.3

Proof of Theorem 5.3. The first part follows from Proposition 5.19, and the second part follows from Proposition 5.27. For the third part, suppose that H satisfies Eqs. (12) to (14) with ε(n) = 1 n− 2 +o(1) . Summarizing, we know that: 1. For all σ ∈ C, n1 wσ (A) −→ mσ ∈ R by assumption. n→∞

2. For all α ∈ E \ C, n1 zα (A) −→ 0 by the first part. n→∞

3. For all α ∈ A \ E, n1 wα (A) −→ 0 by the second part. n→∞

By Lemma 3.14, the traffic distribution of A then exists and is uniquely determined by {mσ : σ ∈ C}, completing the proof.

From Diagrams to Asymptotic GFOM Dynamics

The traffic distribution captures the limiting behavior of all scalar-valued, permutation-invariant polynomials. In this section, we show how to leverage this information to derive the limiting empirical laws of vector-valued, permutation-invariant polynomials. Our main application is a description of the limiting dynamics of GFOM. We will mostly work under the assumption that the input matrices A = A(n) satisfy the strong cactus property, which we recall is the statement that n1 EA zα (A) → 0 as n → ∞ for all non-cactus α (i.e., all α ∈ A \ C, a statement about scalar graph polynomials). In Section 6.3.2 we will briefly suspend this assumption to discuss punctured matrices, so as to connect to the setting of Section 5. We tackle two tasks in this section: 1. First, we study the joint asymptotic limit of the empirical distributions of the vector diagrams zα (A) over α ∈ A1 . Assuming the strong cactus property, we show that only the small subset of treelike α ∈ T1 are asymptotically nonzero in the z-basis, in a sense to be made precise below. We then show that the asymptotic algebra of the treelike diagrams is isomorphic to a Wick algebra, an algebra defined by a family of Gaussian random variables. This will give a precise version of Theorem 1.12.

2. Second, we work with the asymptotic limit of treelike diagrams to identify a generalized Onsager correction, derive a treelike Approximate Message Passing algorithm, and prove its state evolution over arbitrary input matrices having the strong cactus property and a limiting diagonal distribution. This will give a precise version of Theorem 1.13.

6.1

Asymptotic limit of the vector diagrams

In this section, given a family (Xi )i∈I and J ⊆ I, we will write as a shortcut XJ = (Xj )j∈J . Recall that C1 denotes the set of rooted cactuses and T1 ⊆ A1 denotes the set of rooted trees with hanging cactuses. We call the diagrams in T1 treelike, and we call Gaussian trees the subset of diagrams G1 ⊆ T1 such that the root has degree exactly 1 after removing hanging cactuses. Definition 6.1 (Type). For each τ ∈ T1 , let type(τ ) ∈ NG1 ∪C1 , where type(τ )α count the number of copies of α ∈ G1 ∪ C1 attached to the root of τ , with the additional convention that type(τ )α = 0 for all α ∈ G1 that has cactuses hanging at the root. The following theorem identifies the limiting distribution of zA1 (A) under the strong cactus property. We refer the reader to Section C for the definition of convergence in distribution for random elements indexed by countably infinite index sets. Theorem 6.2. Assume that A = A(n) satisfies Eq. (4), has the strong cactus property, and a limiting diagonal distribution. Then, (d)

∞ samp(zA1 (A)) −→ ZA , 1 ∞ ∈ RA1 is a random variable satisfying the following properties: where ZA 1

1. Zα∞ = 0 for all non-treelike α. 2. Conditioned on ZC∞1 , ZG∞1 is a centered Gaussian process with covariance Σ∞ from Eq. (31). 3. Let He denote the Wick product (Definition 2.9). Then for every τ ∈ T1 , Zτ∞ = Hetype(τ ) (ZG∞1 ; Σ∞ ) ·

(Zσ∞ )type(τ )σ .

σ∈C1 ∞ of permutation-invariant, vector-valued polyTheorem 6.2 shows how the limiting algebra ZA 1 ∞ nomials in A can be derived from ZC1 . Although we have not specified the description of the law of ZC∞1 , it is fully determined by the limiting diagonal distribution of A. For example, when A further satisfies the factorizing strong cactus property, ZC∞1 is deterministic:

Proposition 6.3. If A satisfies the factorizing strong cactus property and Eq. (4), then the conclusion of Theorem 6.2 holds with the additional property that for every σ ∈ C1 , Zσ∞ =

Y ρ∈cyc(σ)

1 lim E zρ (A) n→∞ n

Proof. Let σ ∈ C1 . The first moment of samp(zσ (A)) is E samp(zσ (A)) =

n 1X 1 E zσ (A)[i] = E zσ0 (A) , n i=1 A nA

(23)

where σ0 is the unrooted version of σ. As n → ∞, Eq. (23) converges to the deterministic constant Y 1 E zσ0 (A) = n→∞ n A ρ∈cyc(σ )

κσ0 := lim

1 E zρ (A) n→∞ n lim

by the factorizing cactus property. We now switch to the second moment, E samp(zσ (A))2 =

n 1 X E zσ (A)[i]2 . n A i=1

Expand the scalar polynomial ni=1 zσ (A)[i]2 in the z-basis. The support of that expansion is the set of diagrams that can be obtained by grafting two copies of σ at the root and merging pairs of vertices across the two different copies. By the strong cactus property, it suffices to find which cactuses can be obtained in this way. By Lemma D.1, the only cactus that can occur in this way has no merging, and it contributes κ2σ0 by the factorizing cactus property. Thus, P

E samp(zσ (A))2 = κ2σ0 + o(1) = (E samp(zσ (A))2 + o(1) . We showed that samp(zσ (A)) converges to the desired deterministic quantity in expectation, and furthermore that its variance converges to 0. This implies that it converges to the constant in distribution. By unicity of the limit in distribution, Zσ∞ equals that constant almost surely. However, if we drop the factorizing cactus property assumption, the variables ZC∞1 may no longer be deterministic. For example, this can be the case when A is a block-structured matrix as in Section 4.3: (n)

(n)

Example 6.4. Let A1 and A2 Define the 2n × 2n matrix,

be two n × n matrices satisfying the assumptions of Theorem 6.2. " (n)

(2n)

A1 = 0

0 (n) A2

From the block-diagonal structure, for any α ∈ A1 , zα (A)[i] =

 z (A )[i]

if i ∈ [n]

zα (A2 )[i − n]

if i ∈ [2n] \ [n]

Hence, the law of ZC∞1 (A) is a uniform mixture of the law of ZC∞1 (A1 ) and that of ZC∞1 (A2 ). We will prove a generalization of Example 6.4 later; see Lemma 6.31. In Example 6.4, the randomness of ZC∞1 may be viewed as coming solely from the samp(·) operator, but this is not always the case. For instance, our model also captures orthogonally invariant distributions that do not satisfy the traffic concentration property: 47

Example 6.5. Let (λn )n≥1 be an exchangeable sequence of random variables in [−1, 1] and consider A(n) = (Q(n) )⊤ diag(λ1 , . . . , λn )Q(n) , for Haar-distributed matrices Q(n) ∈ O(n), independent from (λn ). By de Finetti’s theorem, there exists a latent random probability measure µ almost surely supported on [−1, 1] such that conditionally on µ, λ1 , λ2 , . . . are i.i.d. with common law µ. By Theorem 4.2, A(n) satisfies the strong cactus property conditionally on µ, so it also satisfies the strong cactus property unconditionally. Applying Theorem 6.2 and Proposition 6.3, we get that conditionally on µ, samp(zC1 (A)) converges in distribution to 



κ|ρ| (µ)

Y  ρ∈cyc(σ)

(24)

σ∈C1

where (κq (µ))q≥1 are the free cumulants of µ. Therefore, unconditionally, samp(zC1 (A)) converges in distribution to the random quantity Eq. (24). Note that Examples 6.4 and 6.5 do not contradict Proposition 6.3 because in these examples, A(n) typically does not satisfy the factorizing cactus property. 6.1.1

Non-treelike diagrams are asymptotically negligible

The remainder of Section 6.1 is dedicated to the proof of Theorem 6.2. In the whole proof, we drop the dependence of zα and wα on A to lighten notation. We start by proving that non-treelike diagrams are negligible. Lemma 6.6. Suppose that A satisfies the strong cactus property. Then for each non-treelike α, L2

samp(zα ) −→ 0 . Proof. By definition, we have E samp(zα )2 =

n h i 1X E zα2 [i] n i=1

By Lemma D.3, we can expand zα2 in the z-basis to obtain, for some constant coefficients cβ =

n X 1X cβ E zβ [i] n i=1 β∈A \T 1

1 X c′ E zβ = n β∈A \T β 0

for some other constant coefficients c′β . Since no diagram in A0 \ T0 is a cactus, by the strong cactus property, we get E samp(zα )2 −→ 0, as desired. n→∞

6.1.2

Asymptotic limit of the treelike diagrams

Next, we analyze the treelike diagrams. All results in Section 6.1.2 are purely combinatorial, meaning that they hold for arbitrary A ∈ Rn×n sym . The covariance of treelike diagrams is defined in terms of homeomorphic matchings between them. We start by defining this new concept. Definition 6.7 (Core). Let τ ∈ T1 . Define core(τ ) to be the rooted tree obtained from τ by 1. Removing all hanging cactuses. 2. Removing all non-root degree-2 vertices and the two edges they are incident with, and adding back a new edge between their two neighbors. Note that the vertex set V (core(τ )) may be identified with a subset of V (τ ), even though the second rule may lead to edges being present in core(τ ) that do not exist in τ . Definition 6.8 (Homeomorphic matchings). Let τ1 , τ2 ∈ T1 . We say that a partial matching P ⊆ V (τ1 ) × V (τ2 ) of τ1 and τ2 is homeomorphic if 1. (root(τ1 ), root(τ2 )) ∈ P . 2. Restricted to V (core(τ1 )) × V (core(τ2 )), P is a rooted graph isomorphism between core(τ1 ) and core(τ2 ). 3. Let {u, u′ } ∈ E(core(τ1 )), let (u = u1 , . . . , uk = u′ ) be the path between u and u′ in τ1 . Let v = P (u), v ′ = P (u′ ), and (v = v1 , . . . , vℓ = v ′ ) be the path between v and v ′ in τ2 . Then there is no matching edge between {u1 , . . . , uk , v1 , . . . , vℓ } and its complement. Moreover, for all (ui , vj ) ∈ P and (ui′ , vj ′ ) ∈ P , we have i ≤ i′ ⇐⇒ j ≤ j ′ (the matching restricted to the vertices in the paths is non-crossing). 4. No inner vertices from the hanging cactuses are matched. We denote by H(τ1 , τ2 ) the set of homeomorphic matchings between τ1 and τ2 . This definition is motivated by the following lemma stating that, when computing the covariance of two treelike diagrams, the matchings giving rise to cactuses are precisely the homeomorphic ones. Lemma 6.9. Let τ1 , τ2 ∈ T1 and τ = τ1 ⊔ τ2 . For any matching P ⊆ V (τ1 ) × V (τ2 ) such that (root(τ1 ), root(τ2 )) ∈ P , we have τP ∈ C1 if and only if P ∈ H(τ1 , τ2 ). In particular, if τ1 , τ2 ∈ C1 , only the matching P = {(root(τ1 ), root(τ2 ))} creates a cactus τP . We are now ready to describe the algebra of treelike diagrams: Lemma 6.10. For all γ1 , . . . , γℓ ∈ G1 ∪ C1 , ℓ Y j=1

zγj −

M ∈M(ℓ) Puv ∈H(γu ,γv ) ∀uv∈M

zL uv∈M

where ⊕ denotes the grafting at the root.

γPu,v ⊕

L u∈M /

γu ∈ span(zA1 \T1 ) ,

(25)

The proofs of Lemmas 6.9 and 6.10 are deferred to Section D.1. Note that the error in Lemma 6.10 is measured in terms of non-treelike diagrams. By inverting Eq. (25), we can formulate the algebra of treelike diagrams in the language of Wick products (Definition 2.9). Corollary 6.11. For all τ ∈ T1 , zτ − Hetype(τ ) (zG1 ; Σ)

(Zσ∞ )type(τ )σ ∈ span(zA1 \T1 ) ,

(26)

σ∈C1

where for all γ, γ ′ ∈ G1 , we defined the “finite-n” covariance matrix Σ[γ, γ ′ ] :=

(27)

zγP .

P ∈H(γ,γ ′ )

Proof. We proceed by induction on the number of vertices of τ . First, Eq. (26) trivially holds if τ has one vertex, which proves the base case. Now, suppose that τ = γ1 ⊕ . . . ⊕ γℓ is the grafting at the root of γ1 , . . . , γℓ ∈ G1 ∪ C1 . By Lemma 6.10, zτ +

zL uv∈M

M ∈M(ℓ) Puv ∈H(γu ,γv ) M ̸=∅ ∀uv∈M

γPu,v ⊕

L u∈M /

γu −

ℓ Y

zγj ∈ span(zA1 \T1 ) .

j=1

Applying the induction hypothesis and using additivity of types, we have: zL

γ ⊕ uv∈M Pu,v

γu − u∈M /

Y uv∈M

zγu · HeP

zγPuv

type(γu ) u∈M / γu ∈G1

u∈M / γu ∈C1

(zG1 ; Σ) ∈ span(zA1 \T1 ) .

Since cactuses are not matched by homeomorphic matchings by definition, the product over cactuses zγu is over all u such that γu ∈ C1 , which is independent of M and can be factorized out. Therefore, in the rest of the proof we assume that γi ∈ G1 for all i ∈ [ℓ]. Using Claim D.4, we obtain zτ +

zγP · HeP

M ∈M(ℓ) uv∈M P ∈H(γu ,γv ) | {z M ̸=∅ Σ[γu ,γv ]

u∈M /

type(γu ) (zG1 ; Σ) −

ℓ Y

zγj ∈ span(zA1 \T1 ) .

(28)

j=1

}

By the recursive formula of the Wick products (Corollary 2.11), X

M ∈M(ℓ) uv∈M M ̸=∅

Σ[γu , γv ] · HeP

u∈M /

type(γj ) (zG1 ; Σ) + Hetype(τ ) (zG1 ; Σ) =

ℓ Y

zγj .

(29)

j=1

Combining Eqs. (28) and (29) concludes the proof. Finally, if we reduce Lemma 6.10 modulo the larger class of non-cactus diagrams (which are the negligible diagrams in expectation under the strong cactus property), we deduce that the joint moments of the diagrams in G1 have an asymptotically Gaussian structure.

Corollary 6.12. For all γ1 , . . . , γℓ ∈ G1 and σ1 , . . . , σk ∈ C1 , 

k Y

zσi 

i=1

ℓ Y

 X

zγj −

j=1

Σ[γx , γy ] ∈ span(zA1 \C1 ) ,

M ∈Mperf (ℓ) xy∈M

where Σ is defined in Eq. (27). Proof. Every non-treelike term in Eq. (25) is a fortiori not a cactus. Also, the only cactuses in the subtracted term occur when M is a perfect matching. In other words,  k Y i=1

ℓ Y

zσi  

 X

zγj −

j=1

M ∈Mperf (ℓ) Puv ∈H(γu ,γv ) ∀uv∈M

uv∈M

γPuv  ∈ span(zA1 \C1 ) .

 

Therefore, by Lemma D.3, we deduce zL uv∈M

γPuv

−

zγPuv ∈ span(zA1 \C1 ) ,

Y uv∈M

and the desired statement follows. 6.1.3

Proof of Theorem 6.2

Claim 6.13. Suppose that the traffic distribution of A exists. Then, for any α1 , . . . , αk ∈ A1 , the sequence E samp(zα1 · · · zαk ) converges as n → ∞. Proof. This is straightforward, as E samp(zα1 · · · zαk ) =

n 1 X zα [i] · · · zαk [i], E n i=1 1

and the inner polynomial is a scalar polynomial of A that can be expanded in the z-basis of scalar diagrams as a linear combination of various quotients of the scalar diagram formed by forgetting the identity of the root in α1 ⊕ · · · ⊕ αk . Claim 6.13 implies in particular that the sequence samp(zA1 ) is tight. In the rest of the proof, we show that the limit in distribution actually exists and characterize it. The following lemma is a direct consequence of the fundamental theorem of graph polynomials. Lemma 6.14. If ∥A∥ ≤ O(1), then for each α ∈ E1 , there exists Cα > 0 such that |samp(zα )| ≤ Cα . Proof. By Lemma 3.9 and Lemma 3.13, we can expand for some coefficients cβ = cβ (α) ∈ R, zα =

cβ wβ .

β∈E1

By Theorem 5.7, it holds for every β ∈ E1 that ∥wβ ∥∞ ≤ ∥A∥|E(β)| , which is at most Oα (1) by assumption. The lemma follows by the triangle inequality. 51

Lemma 6.15. Suppose that the traffic distribution of A exists and that Eq. (4) holds. Then, samp(zC1 ) converges in distribution to some stochastic process ZC∞1 . Proof. First, assume that supn≥1 ∥A(n) ∥ ≤ K holds almost surely, for some universal constant K > 0. All the moments of samp(zC1 ) converge by Claim 6.13. Since cactuses are 2-edge-connected, by Lemma 6.14, all random variables samp(zα ) for α ∈ C1 are uniformly bounded in n. Hence, the moments satisfy the growth condition Eq. (62), so that samp(zC1 ) converges in distribution by Theorem C.2. Finally, if we assume Eq. (4) rather than uniform boundedness, the result can be deduced from the latter case using Lemma C.3. Proof of Theorem 6.2. In the rest of the proof, we assume that the assumptions of Theorem 6.2 are satisfied. We start by analyzing convergence of the subtracted term from Corollary 6.12. By convergence in distribution of the cactuses (Lemma 6.15) and the continuous mapping theorem, we have for any γ1 , . . . , γℓ ∈ G1 and σ1 , . . . , σk ∈ C1 , 

samp 

k Y

 X

zσi

i=1

k Y

(d)

Σ[γx , γy ] −→

Zσ∞i

i=1

M ∈Mperf (ℓ) xy∈M

Σ∞ [γx , γy ] ,

(30)

M ∈Mperf (k) xy∈M

where we defined, for any γ1 , γ2 ∈ G1 , the “limiting” covariance matrix Σ∞ [γ1 , γ2 ] :=

Zγ∞P .

(31)

P ∈H(γ1 ,γ2 )

Since all joint moments converge by Claim 6.13, the sequence of random variables on the left-hand side of Eq. (30) is uniformly integrable. So we also get convergence of the mean, 

E samp 

k Y

i=1

 X

zσi



Σ[γx , γy ] −→ E  n→∞

M ∈Mperf (ℓ) xy∈M

k Y



Zσ∞i

i=1

Σ∞ [γx , γy ] .

M ∈Mperf (k) xy∈M

Combining with Corollary 6.12 and the strong cactus property, 

E samp 

k Y i=1

zσi

ℓ Y j=1





zγj  −→ E  n→∞

k Y

i=1



Zσ∞i

Σ∞ [γx , γy ] .

(32)

M ∈Mperf (ℓ) xy∈M

The right-hand side of Eq. (32) coincides with the moments of ZG∞1 ∪C1 . Recall that the law of ZG∞1 satisfies that after sampling ZC∞1 from its marginal (which is bounded almost surely by Lemma 6.14), then ZG∞1 conditioned on ZC∞1 is a Gaussian process with covariance kernel given by Eq. (31). This object satisfies the moment growth condition Eq. (62). So Theorem C.2 applies and we obtain convergence in distribution of samp(zG1 ∪C1 ) to ZG∞1 ∪C1 . By Lemma 6.6, the non-treelike diagrams converge in L2 to 0, so by Slutsky’s lemma, we obtain joint convergence in distribution, except for the remaining treelike, non-Gaussian trees. By Corollary 6.11, these are continuous images of cactuses and non-treelike diagrams, so by the ∞. continuous mapping theorem, all diagrams converge jointly in distribution to ZA 1

6.2

The treelike AMP algorithm

Now we turn to studying the dynamics of GFOM operations. Definition 6.16 (Asymptotic state). Let (xi )i∈I be a family of random vectors, xi ∈ Rn . We say that a stochastic process (Xi )i∈I is the asymptotic state of (xi )i∈I if, for any k ≥ 1, i1 , . . . , ik ∈ I, and any bounded continuous or polynomial function φ : Rk → R, n 1X E φ(xi1 [j], . . . , xik [j]) = E φ(Xi1 , . . . , Xik ) . n→∞ n j=1

lim

(33)

Definition 6.16 requires in particular for (xi )i∈I to converge in distribution to (Xi )i∈I . As with convergence in distribution in general, this suffers from the caveat that the law of the limit in distribution of (Xi )i∈I is unique, but the probability space on which the limit (Xi )i∈I is realized is not. Thus when we speak of “the asymptotic state” we refer to a specific law, not a specific collection of random variables. Nonetheless, the sampling procedure in Theorem 6.2 suggests a natural way to sample an asymptotic state of the iterates of a pGFOM, since, provided we know how to sample from ZC∞1 (which we must address on a case-by-case basis), the other Zα∞ are conditionally Gaussian or deterministic functions thereof. Translating the limiting variables Zα∞ from Theorem 6.2 to a construction of an asymptotic state, we find: Lemma 6.17. Assume that A = A(n) satisfies the assumptions of Theorem 6.2. Let x=

cα zα (A)

(34)

α∈A1

for some finitely supported coefficients (cα )α∈A1 . Then, X :=

cα Zα∞

(35)

α∈A1

is the asymptotic state of x. Moreover, if xt is of the form Eq. (34) for any t ≥ 1 and Xt is correspondingly defined as in Eq. (35), then (Xt )t≥1 is the asymptotic state of (xt )t≥1 . We emphasize that the index set t ≥ 1 is independent of n, and so our results hold for all fixed iterates t independent of n, in the limit n → ∞. Proof. The statement for bounded continuous test functions φ follows from Theorem 6.2 and the continuous mapping theorem. For polynomial φ, we proceed by a truncation argument. Let Sn := samp(x1 , . . . , xt ) and S := (X1 , . . . , Xt ). Fix a cutoff K > 0 and consider any bounded continuous function φK such that |φK | ≤ |φ|, φK (s) = φ(s) for all ∥s∥2 ≤ K and φK (s) = 0 for all ∥s∥2 > 2K (standard approximations show that such a function exists). First, |E φK (Sn ) − E φK (S)|

converges to 0 as n → ∞ by the bounded continuous case. Next, h

|E φ(Sn ) − E φK (Sn )| ≤ E |φ(Sn )|1∥Sn ∥2 >K

(Definition of the truncated function)

1 2

≤ (E φ(Sn )2 ) Pr(∥Sn ∥2 > K) 1

≤ (E φ(Sn )2 ) 2

(E ∥Sn ∥22 ) K

1 2

(Cauchy-Schwarz inequality)

1 2

(Markov inequality)

Note that these quantities are respectively equal to E φ(Sn )2 =

n 1X φ(x1 [i], . . . , xt [i])2 n i=1

and

E ∥Sn ∥22 =

t n X 1X

n i=1 s=1

xs [i]2 ,

which both converge as n → ∞ by existence of the traffic distribution, and in particular are bounded uniformly in n. Hence, there exists C > 0 independent of n and K such that lim sup |E φ(Sn ) − E φK (Sn )| ≤ n→∞

C , K

and the same holds for E φ(S) − E φK (S) by the same argument, using the fact that all moments ∞ . We obtain lim sup exist in the space generated by ZA n→∞ |E φ(Sn ) − E φK (Sn )| ≤ 2C/K, and 1 the claim follows from taking the limit K → ∞. By Claim 1.5, the iterates of a pGFOM are of the form Eq. (34), so they have an asymptotic state. By definition, these asymptotic states determine the limiting distribution of any (bounded continuous or polynomial) observable. Motivated by this, we introduce a family of approximate message passing algorithms whose asymptotic states are conditionally Gaussian. Theorem 6.18 (Treelike AMP). Assume that A = A(n) satisfies the assumptions of Theorem 6.2. Let ft : R → R be polynomial functions.12 Define: x0 = 1 , bs,t [i] :=

xt = Aft−1 − X

t−1 X

s=0 t−1 Y

bs,t · fs , !

A[ir−1 , ir ]fr′ [ir ]

is ,...,it−1 ∈[n] distinct r=s+1 is =i ft := ft (xt ) , ft′ := ft′ (xt ) ,

A[it−1 , is ] ,

(36)

f0 = 1 .

Then xt ∈ span(zG1 ∪(A1 \T1 ) (A)). Therefore, the asymptotic state of (xt )t≥1 defined in Eq. (35) is a centered Gaussian process conditionally on ZC∞1 . To prove Theorem 6.18, motivated by the results in Section 6.1, we introduce the following handy notations: ∞

Definition 6.19 (Equality modulo non-treelike diagrams). For x, y ∈ span(zA1 ), we write x = y if x − y ∈ span(zA1 \T1 ). We denote by cactus(x) the projection of x onto the span of the cactus diagrams C1 , and by gaussian(x) the projection of x onto the span of the Gaussian diagrams G1 . 12

For ease of exposition ft is assumed to be “memoryless”, meaning that it only takes the most recent xt as input.

The iterates of the treelike AMP algorithm Eq. (36) are engineered to asymptotically generate a self-avoiding walk. That is, whenever the algorithm performs a matrix multiplication operation, the Onsager correction terms in Eq. (36) (the subtracted terms involving bs,t ) are chosen to subtract off the terms in the resulting diagram expansion which (1) are treelike and (2) revisit an existing vertex in any diagram. Example 6.20 (Self-avoiding walk). For intuition, consider the case of Theorem 6.18 where ft (x) = x. Let πt be the t-path diagram and ρt the t-cycle diagram. We can expand exactly: Azπt = zπt+1 +

t X

zρs+1 ⊕πt−s .

s=0

For each term on the right-hand side, we have the approximate factorization (by Lemma 6.10) ∞ zρs+1 ⊕πt−s = zρs+1 · zπt−s , which holds up to non-treelike terms. Then, we define a self-avoiding version of power iteration by: x0 = 1,

xt+1 = Axt −

t X

zρs+1 · xt−s .

s=0 ∞

By construction, we have xt = zπt and therefore, assuming the conditions of Theorem 6.2, the asymptotic state Xt of xt is Gaussian. To analyze a general iteration in the proof of Theorem 6.18, we separate the diagram expansion of ft into its linear and nonlinear parts: ft =

cγ zγ (A) +

X γ∈G1

cτ zτ (A) +

X τ ∈T1 \G1

=:ft1

}

cα zα (A) .

X α∈A1 \T1

=:ft̸=1

}

We call ft1 = gaussian(ft ) the “linear part” since it should be thought of as the degree-1 part of the Hermite expansion of ft with respect to the Gaussian vectors zG1 (A), while ft̸=1 equals all of the other components of the Hermite expansion. More precisely, when ft is of the form ft (xt ) for some Gaussian vector xt , which is the situation for AMP, then ft1 has the following simple form. Lemma 6.21. Let x ∈ span(zG1 ) and let f : R → R be a polynomial. Then, ∞

gaussian(f (x)) = cactus(f ′ (x)) · x . Proof. Suppose that f (x) = xℓ for some integer ℓ ≥ 0 (the general case follows by linearity). By Lemma 6.10, the product of diagrams in G1 yields a diagram in G1 only when every diagram P except one is matched. Formally, write x = γ∈G1 cγ zγ , then by Lemma 6.10, ∞

f1 = ℓ

cγ1 . . . cγℓ

γ1 ,...,γℓ ∈G1

Viewing

M ∈Mperf (ℓ−1) Puv ∈H(γu ,γv ) ∀uv∈M

uv∈M

γPu,v ⊕γℓ .

(37)

uv∈M γPu,v as a fixed cactus, by Lemma 6.10, every term on the right-hand side satisfies

zL uv∈M

∞ L γPu,v ⊕γℓ = z uv∈M γPu,v · zγℓ .

(38)

Applying Lemma 6.10 one last time, it remains to observe that the cactus part of f ′ (x) = ℓxℓ−1 is ∞

cactus(f ′ (x)) = ℓ

cγ1 . . . cγℓ−1

γ1 ,...,γℓ−1 ∈G1

M ∈Mperf (ℓ−1) Puv ∈H(γu ,γv ) ∀uv∈M

uv∈M

γPu,v .

(39)

Combining Eqs. (37) to (39) yields the desired claim. The next key Lemma 6.23 derives an explicit asymptotic formula for the AMP iterates: xt is generated by taking a self-avoiding walk from each nonlinear term fs̸=1 . Definition 6.22. Let ct = cactus(ft′ (xt )). Define the self-avoiding walk matrix Bs,t generated by the iteration between time s and t to be: Bs,t [i, j] :=

A[is , is+1 ] · · · A[it−1 , it ] · cs+1 [is+1 ] · · · ct−1 [it−1 ] .

X is ,...,it ∈[n] distinct is =j, it =i

Recalling Definition 5.4, Bs,t is a linear combination of open cactus matrices in the z-basis (up to non-treelike terms which arise from intersections involving the cs+i ). For example, Bt−1,t equals A with the diagonal elements set to zero. We note the analogy between bs,t and Bs,t , which contain similar self-avoiding walks that return to the start and do not return to the start, respectively. Lemma 6.23. Define xt , ft by Eq. (36) and let ct = cactus(ft′ (xt )). Then for t ≥ 1: t−1 ∞ X

xt =

Bs,t fs̸=1

and

∞

ft = ct ·

s=0

t−1 X

Bs,t fs̸=1 + ft̸=1 .

s=0

Proof. First, note that for a fixed t, the second equation follows from the first: ∞

ft = ft1 + ft̸=1 ∞

= ct · xt + ft̸=1 ∞

= ct ·

t−1 X

(Lemma 6.21)

Bs,t fs̸=1 + ft̸=1

(first equation and Lemma D.3)

s=0

To establish the equations, we use induction on t. In the base case t = 1 we have x1 = Af0 −b0,1 ·f0 = B0,1 f0 as needed. Now, assume that the formulas hold for 0, . . . , t. Denote by Ct the diagonal matrix with entries ct . The equation for ft implies: t−1 ∞ X

Aft =

ACt Bs,t fs̸=1 + Aft̸=1

(40)

s=0

If we expand the matrix product ACt Bs,t we can partition the sum based on whether the matrix

A revisits a vertex already on the walk: (ACt Bs,t )[i, j] =

n X

A[i, k]ct [k]

is ,...,it ∈[n] distinct is =j, it =k

k=1

A[is , is+1 ] · · · A[it−1 , it ] · cs+1 [is+1 ] · · · ct−1 [it−1 ]

A[is , is+1 ] · · · A[it , it+1 ] · cs+1 [is+1 ] · · · ct [it ]

(41)

A[is , is+1 ] · · · A[it−1 , it ] · A[ir , it ] · cs+1 [is+1 ] · · · ct [it ] .

(42)

is ,...,it+1 ∈[n] distinct is =j, it+1 =i

t X

r=s is ,...,it ∈[n] distinct is =j, ir =i

The first term Eq. (41) is self-avoiding and equals Bs,t+1 [i, j]. In the second term Eq. (42), the term r = s is diagrammatically a cycle and is equal to bs,t+1 [i] when i = j, and 0 otherwise: Claim 6.24. We have:  ∞

bs,t =  

 X

A[is , is+1 ] · · · A[it−2 , it−1 ]A[it−1 , is ] · cs+1 [is+1 ] · · · ct−1 [it−1 ]  

is ,...,it−1 ∈[n] distinct is =i

i∈[n]

Proof. The only difference between this formula and the definition of bs,t is that the vectors ft′ at the internal vertices of the cycle have been replaced by ct . This holds up to non-treelike terms since placing a non-cactus diagram at any internal vertex of the cycle will create only non-treelike diagrams. The remaining terms in Eq. (42) are a cycle and a path joined together at vertex r: Claim 6.25. Let r ∈ {s + 1, . . . , t}. For i ∈ [n], let u[i] =

A[is , is+1 ] · · · A[it−1 , it ] · A[ir , it ] · cs+1 [is+1 ] · · · ct [it ]fs̸=1 [is ] .

X is ,...,it ∈[n] distinct ir =i

∞

Then u = br,t+1 · Cr Bs,r fs̸=1 . Proof. By expanding definitions, we can conveniently interpret br,t+1 · Cr Bs,r fs̸=1 [i] =

A[is , is+1 ] · · · A[it−1 , it ] · A[ir , it ] · cs+1 [is+1 ] · · · ct [it ]fs̸=1 [is ] .

is ,...,it ∈[n] is ,...,ir distinct ir ,...,it distinct ir =i

Since the diagram induced on {ir , . . . , it } is a cycle, any intersection between the vertices {is , . . . , ir } and {ir , . . . , it } would create a non-treelike diagram.

Plugging Claim 6.25 in to Eq. (40), we have: t−1 ∞ X

Aft =

Bs,t+1 fs̸=1 +

s=0

= ∞

t X s=0 t X

t−1 X

bs,t+1 · fs̸=1 +

s=0

̸ 1 Bs,t+1 fs= +

Bs,t+1 fs̸=1 +

s=0

t X r=0 t X

t−1 X t X s=0 r=s+1

br,t+1 · Cr

r−1 X

Aft̸=1

br,t+1 · (Cr Bs,r fs̸=1 ) +

| {z }

=Bt,t+1 ft̸=1 +bt,t+1 ·ft̸=1

! ̸ 1 Bs,r fs= + fr̸=1

s=0

br,t+1 · fr

r=0

The last equality is the inductive formula for fr . The Onsager correction subtracts off the second sum, leaving only the desired first sum for xt+1 . Proof of Theorem 6.18. We prove the following purely combinatorial claim about Eq. (36): xt is in the span of non-treelike diagrams and Gaussian treelike diagrams. By Lemma 6.17, this will imply that conditioned on ZC∞1 , the asymptotic state of (xt )t≥1 is a centered Gaussian process, as desired. To show the claim, we start from the conclusion of Lemma 6.23: t−1 ∞ X

xt =

Bs,t fs̸=1 .

s=0

The diagrams in Bs,t fs̸=1 are obtained by: (1) choose a diagram from fs̸=1 , (2) choose a cactus diagram from cr at each internal vertex of Bs,t (i.e. each internal vertex along a path of length t − s), (3) multiply these diagrams together. Since none of the diagrams in fs̸=1 have degree 1 at the root by definition, the only treelike terms in the product are formed by grafting the diagrams together without intersections. In particular, the root is the endpoint of the path in Bs,t and has degree 1. This concludes the proof. 6.2.1

Covariance structure of treelike AMP

While Theorem 6.18 shows that the treelike AMP iterates are asymptotically Gaussian, it does not identify their covariance. We calculate the covariance “combinatorially” by calculating the cactus diagrams appearing in the expansion of ⟨xs , xt ⟩. Proposition 6.26. Let xt follow the iteration Eq. (36). Then for any s, t ≥ 1, xs · xt −

s−1 t−1 X X

Bs′ st′ t (fs′ · ft′ ) ∈ span(zA1 \C1 ) ,

(43)

s′ =0 t′ =0

where for 0 ≤ s′ ≤ s, 0 ≤ t′ ≤ t, we define the matrix Bs′ st′ t ∈ Rn×n by Bs′ st′ t [i, j] :=

s−1 Y

A[ir , ir+1 ]

is′ ,...,is ,jt′ ,...,jt ∈[n] r=s′ distinct except is′ =jt′ =j, is =jt =i

t−1 Y r=t′

A[jr , jr+1 ]

s−1 Y r=s′ +1

cr [ir ]

t−1 Y

cr [jr ] .

r=t′ +1

When we apply this lemma, we will average Eq. (43) over the coordinates i ∈ [n] and over A. 58

Since the error terms are all in the span of non-cactus diagrams, all error terms converge to 0 by the strong cactus property. On the other hand, the average of the term xs · xt converges to the covariance E[Xs Xt ] which we want to calculate. The subtracted terms involving the Bs′ t′ st matrices converge to limits depending on the asymptotic values of the cactuses ZC∞1 . For some settings (such as Proposition 6.3), the values ZC∞1 are deterministic. For other settings, the values ZC∞1 are random, and we will condition on them in order to obtain the conditional covariance. Proof. Since xs and xt have degree exactly one at the root, in order to form a cactus in xs · xt , the paths from the root of xs and xt in the expansion from Lemma 6.23 must meet at some point. This intersection cannot happen at a vertex from fs̸=′ 1 or ft̸=′ 1 (that would create edges in two cycles). Let s′ ∈ {0, . . . , s − 1} and t′ ∈ {0, . . . , t − 1} denote the integers such that the first intersection corresponds to the indices is′ (for xs ) and it′ (for xt ) in Definition 6.22. Then, we can decompose xs · x t −

s−1 t−1 X X

Bs′ st′ t ((cs′ · xs′ + fs̸=′ 1 ) · (ct′ · xt′ + ft̸=′ 1 )) ∈ span(zA1 \C1 ) ,

s′ =0 t′ =0 ∞

and the conclusion follows from the equality (Lemma 6.21) fs = cs · xs + fs̸=1 . The cactus expansion of Bs′ st′ t (fs′ · ft′ ) can be obtained explicitly by combining a cycle of length s − s′ + t − t′ along the edges of Bs′ st′ t , a cactus from cr hanging at every vertex r in the cycle, and a homeomorphic matching of the tree components of fs′ and ft′ (Definition 6.8).

6.3

Examples of state evolution

In this section, we specialize Theorem 6.18 to obtain a more explicit description of the state evolution of the treelike AMP algorithm for several concrete matrix models. Notation 6.27. For a vector x ∈ Rn , we will use the following notation for empirical averages: n 1X xi . ⟨x⟩ := n i=1

Technically, most algorithms in this section are not pGFOM since they calculate empirical averages. Assuming that the traffic distribution concentrates and the vector x lies in the diagram basis, then the empirical average ⟨x⟩ concentrates, and we can replace ⟨x⟩ by its limit E X without changing the asymptotic state of the algorithm. This is formally proven in Lemma D.10. 6.3.1

Orthogonally invariant random matrices

In the special case that A is drawn from an orthogonally invariant random matrix ensemble, the treelike AMP algorithm recovers the orthogonal AMP algorithm of Fan [Fan22], giving a new proof of this result. Theorem 6.28 (State evolution for orthogonally invariant matrices). Let A = A(n) ∈ Rn×n sym be an orthogonally invariant random matrix converging in tracial moments in L2 to a probability measure with free cumulants (κq )q≥1 . Assume A satisfies Eq. (4). Let ft : R → R be polynomial functions

and define the iteration x0 = 1,

xt = Aft−1 (xt−1 ) −

t−1 X

t−1 Y

κt−s

s=0

⟨fr′ (xr )⟩

fs (xs )

∀t ≥ 1 .

(44)

r=s+1

Then, the asymptotic state of (xt )t≥1 is a centered Gaussian process (Xt )t≥1 with covariance E [Xs Xt ] =

s−1 t−1 X X s′ =0 t′ =0



κs−s′ +t−t′ 

s−1 Y



E fr′ (Xr ) 

r=s′ +1

t−1 Y



E fr′ (Xr ) E [fs′ (Xs′ )ft′ (Xt′ )]

∀s, t ≥ 1 ,

r=t′ +1

with X0 := 1. Proof. By Theorem 4.2, A satisfies the factorizing strong cactus property and its diagonal distribution exists, so the assumptions of Theorem 6.2 and Theorem 6.18 are satisfied. Therefore, the treelike AMP algorithm in Eq. (36) has Gaussian asymptotic state. We now specialize the Onsager correction term in Eq. (36) to this model. The term bs,t is represented by a cycle of length t − s, with fr′ (xr ) attached to the rth vertex of the cycle for each s < r < t. By Lemma D.3, we only need to look at treelike contributions in bs,t . Because of the base cycle, these are only cactuses, obtained by attaching cactuses from (fr′ (xr ))s<r<t along the Q ′ base cycle. By Proposition 6.3, bs,t has constant asymptotic state equal to κt−s t−1 r=s+1 E fr (Xr ). The cactuses in bs,t persist until the end of the algorithm, so that they will eventually contribute this value towards the asymptotic state. Hence it does not affect the asymptotic state to replace bs,t immediately by its limiting constant value. Moreover, by Lemma D.10 and Lemma B.7, we may replace E fr′ (Xr ) by the empirical average ′ ⟨fr (xr )⟩ to obtain Eq. (44) without affecting the asymptotic state. Now the asymptotic state Xt

of Eq. (44) matches that of Eq. (36), and we may apply Theorem 6.18 to deduce that Xt is Gaussian. To calculate the covariance E [Xs Xt ], we average Proposition 6.26 over the coordinates i ∈ [n] and take the limit n → ∞. On the right side of Eq. (43), the cycle of Bs′ st′ t contributes κs−s′ +t−t′ and the hanging diagrams fr′ (xr ) inside Bs′ st′ t contribute E fr′ (Xr ) by the factorizing cactus property. The cactuses in fs′ (xs′ ) · ft′ (xt′ ) contribute E [fs′ (Xs′ )ft′ (Xt′ )], which establishes the desired recurrence. Note that this proof only uses the strong factorizing cactus property and the concentration of the traffic distribution, which explains why Theorem 6.28 also holds for non-orthogonally invariant matrix models such as Wigner matrices (Section 4.1). 6.3.2

Punctured random and deterministic matrices

The punctured matrices studied in Section 5 do not satisfy the strong cactus property, so we cannot directly apply Theorem 6.18 to derive an AMP iteration for them. However, a reduction allows us to derive the state evolution of punctured orthogonally invariant random matrices from that of their unpunctured counterparts. These matrices are central because, by Theorem 5.3, they provide an intermediate step in deriving the state evolution of sequences of punctured deterministic matrices satisfying Assumption 5.2.

Note that a GFOM run on a punctured matrix must be initialized with a random vector x0 ∼ N (0, I), rather than x0 = 1, to avoid triviality. Theorem 6.29 (State evolution for punctured matrices). Let H = H (n) ∈ Rn×n sym be a sequence of orthogonally invariant random matrices satisfying Eq. (4) and converging in tracial moments in L2 to a probability measure with free cumulants (κq )q≥1 . Let A denote the puncturing of H (Definition 2.1). Let ft : R → R be polynomial functions with f0 (x) = x, and consider the pGFOM: x0 ∼ N (0, I) ,

xt = Aft−1 (xt−1 ) −

t−1 X

t−1 Y

κt−s

s=0

⟨fr′ (xr )⟩

(fs (xs ) − ⟨fs (xs )⟩1)

∀t ≥ 1 .

r=s+1

Then for any t ≥ 1 and any polynomial φ : Rt → R, we have lim

E ⟨φ(x1 , . . . , xt )⟩ = E φ(X1 , . . . , Xt ) ,

n→∞ H,x0

where (Xt )t≥1 is a centered Gaussian process with covariance given by E [Xs Xt ] =

s−1 t−1 X X



κs−s′ +t−t′ 

s′ =0 t′ =0

F0 := 1 ,

s−1 Y



E fr′ (Xr ) 

r=s′ +1

t−1 Y

 h

E fr′ (Xr ) E Fs′ Ft′

∀s, t ≥ 1 ,

r=t′ +1

Ft := ft (Xt ) − E ft (Xt )

∀t ≥ 1 .

By Theorem 5.1, the conclusion of Theorem 6.29 also holds for any sequence of deterministic matrices satisfying the delocalization assumption Assumption 5.2 and having a limiting diagonal distribution that factorizes over cycles (that is, matches the diagonal distribution of some orthogonally invariant random matrix ensemble). In particular, the conclusion holds for the Walsh-Hadamard matrices and the Discrete Cosine and Sine Transform matrices, for which the κq are the free cumulants of the ROM (Eq. (9)). The proof of Theorem 6.29 proceeds by reducing to the following iteration on the original, non-punctured matrix, initialized at the all-ones vector: u0 = 1 , where bs,t [i] :=

ut = Hft−1 −

f0 := u0 ,

bs,t · fs

∀t ≥ 1 ,

s=0 t−1 Y

H[ir−1 , ir ]fr′ [ir ]

X is ,...,it−1 ∈[n] distinct is =i

t−1 X

H[it−1 , is ]

∀t > s ≥ 0 ,

(45)

r=s+1

ft := Πft (ut )

∀t ≥ 1 .

Lemma 6.30. For any t ≥ 1 and any polynomial φ : Rt → R, lim

E ⟨φ(x1 , . . . , xt )⟩ = lim E ⟨φ(u1 , . . . , ut )⟩ .

n→∞ H,x0

n→∞ H

The proof of Lemma 6.30 is deferred to Section D.3. Proof of Theorem 6.29. We apply Theorem 6.28 to ut after replacing iteratively each occurrence of Πft (ut ) by ft (ut )−E [ft (Ut )]·1 (where Ut is the asymptotic state of ut as predicted by Theorem 6.28). 61

By Lemma D.10, this transformation does not change the asymptotic state of ut . The state evolution formula for polynomial test functions then transfers to xt by Lemma 6.30. 6.3.3

Block-structured random matrices

Our final example is the class of block-structured matrices whose blocks satisfy the factorizing strong cactus property, which we introduced in Section 4.3. As anticipated in Example 6.4, these matrices do not themselves satisfy the factorizing strong cactus property. Therefore, we start by describing the random limit ZC∞1 . (n)

×n

q q Lemma 6.31. Let q ∈ N. For r, c ∈ [q], let Ar,c = Ar,c ∈ Rsym be a sequence of symmetric random n×n matrices such that Ar,c = Ac,r . Let A ∈ Rsym be the block matrix with blocks (Ar,c )r,c∈[q] . Assume that each Ar,c satisfies the factorizing strong cactus property, and (Ar,c )1≤r≤c≤q are asymptotically {r,c} traffic independent. Let (κℓ )ℓ≥1 be the limiting free cumulants of Ar,c . Then,

(d)

(block(i), zC1 (A)[i]) −→ (R, ZC∞1 (R)) ,

i ∼ Unif([n]) ,

R ∼ Unif([q]) ,

where the (deterministic) sequence ZC∞1 (r) for r ∈ [q] is defined recursively by: ∞ (i) For the singleton cactus, Zsingleton (r) := 1.

(ii) Suppose σ ∈ C1 is rooted at a vertex u1 of degree 2. Let (u1 , . . . , uℓ ) be the cycle incident to the root. Let σ2 , . . . , σℓ ∈ C1 be the rooted cactuses attached to the vertices of the cycle. Then    ℓ ℓ  Y   X  {r,c} Y ∞    Zσk (r) Zσ∞k (c) κℓ  

if ℓ is even

 ℓ  Y   κ{r,r} Zσ∞ (r)  ℓ

if ℓ is odd

Zσ∞ (r) := c∈[q] 

k=2 k odd

k=2 k even

k=2

(iii) If σ ∈ C1 decomposes as σ =

Qℓ ∞ ∞ k=1 Zσk (r). k=1 σk , then Zσ (r) :=

Lℓ

In particular, the law of the limit ZC∞1 is Unif({ZC∞1 (r) : r ∈ [q]}). The proof is deferred to Section D.4. By unicity of the limit in distribution, Lemma 6.31, together ∞ . The joint convergence with block(i) clarifies the with Theorem 6.2, determines the law of ZA 1 source of the randomness of ZC∞1 : it arises because of the random choice of the block an entry belongs to in the samp(·) operation. Using Lemma 6.31, we can specialize the treelike AMP iteration and its state evolution to concrete block-structured models. We start with block GOE matrices (Definition 4.4). A family of AMP iterations for such matrices was derived in [Ran11, JM13]. As we will discuss below, these iterations have the same asymptotic state as treelike AMP. Theorem 6.32 (State evolution for the block GOE model). Let A ∼ BlockGOE(n, Σ), where Σ ∈ Rq×q ≥0 is a symmetric matrix. Given polynomial functions ft : R → R, let x0 = 1 ,

x1 = Af0 (x0 ) ,

′ xt = Aft−1 (xt−1 ) − (A⊙2 ft−1 (xt−1 )) · ft−2 (xt−2 )

∀t ≥ 2 .

(46)

Then the asymptotic state of (xt )t≥1 is a mixture 1q r∈[q] µr where µr denotes the law of a centered Gaussian process (Xt )t≥1 with covariance kernel Γr defined recursively by P

Γr [s, t] =

Σ[r, c]

c∈[q]

(XT )T ≥1 ∼µc

[fs−1 (Xs−1 )ft−1 (Xt−1 )]

∀r ∈ [q] ,

∀s, t ≥ 1 ,

with X0 := 1. Proof. By the discussion after Theorem 4.7, A has a traffic distribution and satisfies the strong cactus property13 , so it satisfies the assumption of Theorem 6.18. We consider the treelike AMP iteration xt in Eq. (36) applied to A. We show that this iteration has the same asymptotic state as Eq. (46) by simplifying the Onsager correction term. The free cumulants of the GOE are 0 except for κ2 , so by Lemma 6.31, the only asymptotically non-negligible cactuses are those such that every cycle is a 2-cycle. For any s < t − 2, bs,t contains an injective cycle of length larger than 2 that cannot be destroyed by later operations. For s = t − 2, we have: bt−2,t [i] =

n X

′ ′ ′ A[i, j]2 ft−1 [j] = (A⊙2 ft−1 )[i] − A[i, i]2 ft−1 [i] .

j=1 j̸=i ′ [i] and b Both A[i, i]2 ft−1 t−1,t contain a self-loop that also cannot be destroyed by later operations. In conclusion, the treelike AMP algorithm from Eq. (36) and the iteration in Eq. (46) are equal up to negligible diagrams. By Theorem 6.18, the asymptotic state (Xt )t≥1 of (xt )t≥1 in Eq. (46) exists and is Gaussian conditionally on ZC∞1 , and so, in the construction from Lemma 6.31, it is Gaussian conditionally on the random variable R.

Next, we specialize the covariance formula given by Proposition 6.26. Since only cactuses of 2-cycles are nonzero in the traffic distribution of A (this may be induced from Lemma 6.31), only the term for s′ = s − 1 and t′ = t − 1 is non-negligible in the expansion of n1 E xs · xt given by Proposition 6.26. The expansion into cactuses of that term is obtained by grafting together a 2-cycle at the root, and cactuses of 2-cycles from fs−1 and ft−1 at the child of the root. Applying the recursive formula for ZC∞1 (r) in Lemma 6.31, we obtain: E [Xs Xt | R = r] =

Σ[r, c] E [fs−1 (Xs−1 )ft−1 (Xt−1 ) | R = c]

∀r ∈ [q] .

c∈[q]

Thus, we have shown that, conditionally on R ∼ Unif([q]), (Xt )t≥1 is a Gaussian process with the required covariance. The result follows by taking µr to be the law of (Xt )t≥1 conditionally on R = r. To illustrate the modularity of our approach, we also study a different block-structured matrix model whose blocks are not all GOE. n

×n

q q Theorem 6.33 (State evolution for the community model). Let M ∈ Rsym be an orthogonally invariant random matrix converging in tracial moments to a probability measure with free cumulants (κq )q≥1 such that κ2 = 1q . Let A be the random symmetric n × n matrix with blocks (Ar,c )r,c∈[q] 13

One can also verify that it satisfies Eq. (4). But note that this was only needed in the proof of Theorem 6.2 to ensure the existence of the limit ZC∞1 , which we established directly in Lemma 6.31.

given by A1,1 = M and for all 1 ≤ r ≤ c ≤ q with (r, c) ̸= (1, 1), Ar,c are i.i.d. nq × nq GOE matrices with entries of variance n1 (and we set Ar,c = Ac,r ). Let xt be the treelike AMP iteration Eq. (36) run on A with arbitrary polynomial nonlinearities. Then the asymptotic state (Xt )t≥1 of (xt )t≥1 is the mixture (1 − 1q )µ0 + 1q µ1 , where µi is the law of a centered Gaussian process (Xt )t≥1 with covariance kernel Γi defined recursively by, for all s, t ≥ 1: Γ0 [s, t] = E [Fs−1 Ft−1 ] Γ1 [s, t] = E [Fs−1 Ft−1 ] +



s−1 t−1 X X

κs−s′ +t−t′ 

s′ =0 t′ =0



s−1 Y

E Fr′  

µ1

r=s′ +1



t−1 Y

E Fr′  E [Fs′ Ft′ ] ,

r=t′ +1

µ1

(s′ ,t′ )̸=(s−1,t−1)

Ft′ := ft′ (Xt ),

Ft := ft (Xt ),

X0 = 1 ,

where Eµ1 denotes expectation with respect to (Xt )t≥1 ∼ µ1 . Proof. The assumptions of Lemma 6.31 and Theorem 6.33 are satisfied. All blocks except the one in position (1, 1) have the same free cumulants (the GOE free cumulants, normalized so that κ2 = 1q ). Therefore, in the construction of Lemma 6.31, we have ZC∞1 (r) = ZC∞1 (s) for all r, s > 1. Let µ0 (resp. µ1 ) be the law of the asymptotic state (Xt )t≥1 of the treelike AMP iteration (xt )t≥1 conditioned on R > 1 (resp. R = 1). By Theorem 6.18, both µ0 and µ1 are the laws of centered Gaussian processes. It remains to specialize the formula of Proposition 6.26 for their covariance to the present setting. Conditionally on R > 1 (that is, outside the community), only 2-cycles at the root contribute to ZC∞1 . Thus, by combining the strong cactus property, Proposition 6.26, and Lemma 6.31, we obtain 1 (E [fs−1 (Xs )ft−1 (Xt ) | R = 1] + (q − 1) E [fs−1 (Xs )ft−1 (Xt ) | R > 1]) q = E [fs−1 (Xs )ft−1 (Xt )] .

E [Xs Xt | R > 1] =

Conditionally on R = 1 (that is, inside the community), we also obtain a contribution of E [fs−1 (Xs )ft−1 (Xt )] from the term s′ = s − 1 and t′ = t − 1 in Proposition 6.26 (again using the normalization κ2 = 1q inside the community). For all of the remaining terms s′ , t′ , when s − s′ + t − t′ is an even integer larger than 2, we obtain a contribution only from c = 1 in Lemma 6.31, namely κs−s′ +t−t′ E [Fs′ Ft′ | R = 1]

s−1 Y

Y t−1

E Fr′ | R = 1

E Fr′ | R = 1 .

r=t′ +1

r=s′ +1

When s − s′ + t − t′ is odd, Lemma 6.31 yields exactly the same expression as the even case. Altogether, we obtain the recursion E [Xs Xt | R = 1] = E [Fs−1 Ft−1 ] + s−1 t−1 X X

κs−s′ +t−t′ E [Fs′ Ft′ | R = 1]

s′ =0 t′ =0 (s′ ,t′ )̸=(s−1,t−1)

s−1 Y r=s′ +1

Y t−1

E Fr′ | R = 1

E Fr′ | R = 1 .

r=t′ +1

These are the desired covariance formulas for µ0 and µ1 , and the mixing weights of the events (R = 1) and (R > 1) are indeed 1q and 1 − 1q , respectively. 64

6.3.4

Further extensions

There are several possible technical extensions of the methods we have developed here, whose full development is left for future work. First, Lemma 6.31 applies to general orthogonally invariant distributions within the blocks, not just the GOE. In principle, one can then derive a corresponding state evolution formula mechanically for non-identically distributed orthogonally invariant blocks with arbitrary free cumulants, although the resulting expression is quite complicated. Second, for technical reasons, we assumed that the blocks are square and symmetric, so that we could work with undirected graphs. The results of [Mal20, CDM24] extend to general matrices, and our techniques should also extend to the setting of varying block sizes and asymmetric matrices, leading to non-uniform mixtures in the recursion for the covariance kernel. One caveat of the treelike AMP algorithm is that the Onsager correction term in Eq. (36) is not obviously efficient to compute in practice.14 On the other hand, the vectors bs,t have asymptotically constant entries in many settings, so that the Onsager correction can be replaced by a simpler asymptotically equivalent term, like in Theorems 6.28 and 6.29. This should also hold for blockstructured models, as in the generalized AMP algorithm of Javanmard and Montanari [JM13]. For example, Eq. (46) is expected to be asymptotically equivalent to: x0 = 1 ,

xt = Aft−1 (xt−1 ) − bt−2,t · ft−2 (xt−2 ) , bt−2,t [i] =

q X

Σ[block(i), c]⟨f ′ (xt−1 ) · 1block=c ⟩ ,

c=1

where 1block=c ∈ {0, 1}n indicates the entries in block c ∈ [q]. The treelike AMP algorithm for Theorem 6.33 is expected to be asymptotically equivalent to: x0 = 1 ,

′ xt = Aft−1 − ⟨ft−1 ⟩ft−2 −

t−1 X s=0 s̸=t−2

ft := ft (xt ) ,

κt−s

t−1 Y

⟨fr′ · 1block=1 ⟩ fs · 1block=1 ,

r=s+1

ft′ := ft′ (xt ) ,

where 1block=1 ∈ {0, 1}n indicates the entries in the first block. Because these expressions involve the blockwise indicators 1block=c , they could be represented and analyzed using an extended diagram basis in which certain indices are constrained to lie in a prescribed block. We leave the full development of this extension to future work. A final open question is to characterize traffic distributions satisfying the (not necessarily factorizing) strong cactus property. Sequences of block matrices with orthogonally invariant blocks provide one general construction of matrices with the strong cactus property. If a sequence of matrices has the strong cactus property, must its traffic distribution arise as the limit (in an appropriate sense) of traffic distributions of block matrices with orthogonally invariant blocks (allowing the number of blocks to tend to infinity)? 14 The bs,t can be approximated with high probability to negligible error for all 0 ≤ s < t ≤ T in time 2O(T ) poly(n) using the color coding technique [AYZ95, DSS20], but the exponential dependence on T makes this algorithm impractical to implement for large T .

References [ABKZ20]

Alia Abbara, Antoine Baker, Florent Krzakala, and Lenka Zdeborová. On the universality of noiseless linear estimation with respect to the measurement matrix. Journal of Physics A: Mathematical and Theoretical, 53(16):164001, 2020. 6, 12

[AGZ10]

Greg W. Anderson, Alice Guionnet, and Ofer Zeitouni. An introduction to random matrices. Cambridge University Press, 2010. 24

[AMS21]

Ahmed El Alaoui, Andrea Montanari, and Mark Sellke. Optimization of mean-field spin glasses. Annals of Probability, 49(6):2922–2960, 2021. 2

[AYZ95]

Noga Alon, Raphael Yuster, and Uri Zwick. Color-coding. Journal of the ACM, 42:844–856, 1995. 65

[Ban10]

Teodor Banica. The orthogonal Weingarten formula in compact form. Letters in Mathematical Physics, 91(2):105–118, 2010. 76

[BHX25]

Zhigang Bao, Qiyang Han, and Xiaocong Xu. A leave-one-out approach to approximate message passing. Annals of Applied Probability, 35(4):2716–2766, 2025. 12

[Bil95]

Patrick Billingsley. Probability and measure. John Wiley & Sons, third edition, 1995. 86

[BIPZ78]

Edouard Brézin, Claude Itzykson, Giorgio Parisi, and Jean-Bernard Zuber. Planar diagrams. Communications in Mathematical Physics, 59:35–51, 1978. 70

[BLM15]

Mohsen Bayati, Marc Lelarge, and Andrea Montanari. Universality in polytope phase transitions and message passing algorithms. Annals of Applied Probability, 25(2):753–822, 2015. 3, 11, 88

[BM08]

Adrian Bondy and U.S.R. Murty. Graph Theory. Springer, 2008. 33

[BM11]

Mohsen Bayati and Andrea Montanari. The dynamics of message passing on dense graphs, with applications to compressed sensing. IEEE Transactions on Information Theory, 57(2):764–785, 2011. 1, 11

[Bol14]

Erwin Bolthausen. An Iterative Construction of Solutions of the TAP Equations for the SherringtonKirkpatrick Model. Communications in Mathematical Physics, 325(1):333–366, 2014. 11

[BS10]

Zhidong Bai and Jack W. Silverstein. Spectral analysis of large dimensional random matrices, volume 20. Springer, 2010. 32

[BSK15]

Jean Barbier, Christophe Schulke, and Florent Krzakala. Approximate message-passing with spatially coupled structured operators, with applications to compressed sensing and sparse superposition codes. Journal of Statistical Mechanics: Theory and Experiment, 2015. 3, 6, 11

[CDM24]

Guillaume Cébron, Antoine Dahlqvist, and Camille Male. Traffic distributions and independence II: Universal constructions for traffic spaces. Documenta Mathematica, 29(1):39–114, 2024. 7, 25, 65, 75, 78, 84

[CGH+ 17] Jordan S. Cotler, Guy Gur-Ari, Masanori Hanada, Joseph Polchinski, Phil Saad, Stephen H. Shenker, Douglas Stanford, Alexandre Streicher, and Masaki Tezuka. Black holes and random matrices. Journal of High Energy Physics, 2017(5):1–54, 2017. 70 [CHS24]

Nicola Muca Cirone, Jad Hamdan, and Cristopher Salvi. Genus expansion for non-linear random matrix ensembles with applications to neural networks. arXiv preprint arXiv:2407.08459, 2024. 1

[CL21]

Wei-Kuo Chen and Wai-Kit Lam. Universality of approximate message passing algorithms. Electronic Journal of Probability, 26:1–44, 2021. 11

[CMP+ 23] Patrick Charbonneau, Enzo Marinari, Giorgio Parisi, Federico Ricci-Tersenghi, Gabriele Sicuro, Francesco Zamponi, and Marc Mézard. Spin Glass Theory and Far Beyond: Replica Symmetry Breaking after 40 Years. World Scientific, 2023. 3 [CMW20]

Michael Celentano, Andrea Montanari, and Yuchen Wu. The estimation error of general first order methods. In Conference on Learning Theory (COLT 2020), pages 1078–1141. PMLR, 2020. 1, 2, 11

[ÇO19]

Burak Çakmak and Manfred Opper. Memory-free dynamics for the Thouless-Anderson-Palmer equations of Ising models with arbitrary rotation-invariant ensembles of random coupling matrices. Physical Review E, 99(6):062140, 2019. 6, 12

[CŚ06]

Benoît Collins and Piotr Śniady. Integration with respect to the Haar measure on unitary, orthogonal and symplectic group. Communications in Mathematical Physics, 264(3):773–795, 2006. 76

[DFGZJ95] Philippe Di Francesco, Paul Ginsparg, and Jean Zinn-Justin. 2D gravity and random matrices. Physics Reports, 254(1-2):1–133, 1995. 70 [DJM13]

David L. Donoho, Adel Javanmard, and Andrea Montanari. Information-theoretically optimal compressed sensing via spatial coupling and approximate message passing. IEEE Transactions on Information Theory, 59(11):7434–7464, 2013. 11

[DLS23]

Rishabh Dudeja, Yue M. Lu, and Subhabrata Sen. Universality of approximate message passing with semirandom matrices. Annals of Probability, 51(5):1616–1683, 2023. 3, 6, 11, 12, 29

[DMM09]

David L. Donoho, Arian Maleki, and Andrea Montanari. Message-passing algorithms for compressed sensing. Proceedings of the National Academy of Sciences, 106(45):18914–18919, 2009. 1, 3, 11

[DMR14]

Yash Deshpande, Andrea Montanari, and Emile Richard. Cone-constrained principal component analysis. In Advances in Neural Information Processing Systems (NIPS 2014), pages 2717–2725, 2014. 2

[DSL24]

Rishabh Dudeja, Subhabrata Sen, and Yue M. Lu. Spectral universality in regularized linear regression with nearly deterministic sensing matrices. IEEE Transactions on Information Theory, 2024. 11

[DSS20]

Jingqiu Ding, Hopkins Samuel, and David Steurer. Estimating Rank-One Spikes from Heavy-Tailed Noise via Self-Avoiding Walks. In Advances in Neural Information Processing Systems (NeurIPS 2020), volume 33, pages 5576–5586, 2020. 65

[EM03]

Nicholas M. Ercolani and Kenneth McLaughlin. Asymptotics of the partition function for random matrices via Riemann-Hilbert techniques and applications to graphical enumeration. International Mathematics Research Notices, 2003:755–820, 2003. 75

[Eti24]

Pavel Etingof. Mathematical ideas and notions of quantum field theory. arXiv preprint arXiv:2409.03117, 2024. 72

[Fan22]

Zhou Fan. Approximate message passing algorithms for rotationally invariant matrices. Annals of Statistics, 50(1):197–224, 2022. 3, 10, 12, 59

[FVRS22]

Oliver Y. Feng, Ramji Venkataramanan, Cynthia Rush, and Richard J. Samworth. A Unifying Tutorial on Approximate Message Passing. Foundations and Trends in Machine Learning, 15(4):335–536, 2022. 1, 3, 11

[GHN26]

Mohammed-Younes Gueddari, Walid Hachem, and Jamal Najim. Approximate Message Passing for General Non-Symmetric Random Matrices. Journal of Theoretical Probability, 39, 2026. 12

[GP13]

Stavros Garoufalidis and Ionel Popescu. Analyticity of the planar limit of a matrix model. Annales Henri Poincaré, 14:499–565, 2013. 75

[GTM+ 24] Cedric Gerbelot, Emanuele Troiani, Francesca Mignacco, Florent Krzakala, and Lenka Zdeborová. Rigorous Dynamical Mean-Field Theory for Stochastic Gradient Descent Methods. SIAM Journal on Mathematics of Data Science, 6:400–427, 2024. 1 [Hol81]

Ian Holyer. The NP-completeness of some edge-partition problems. SIAM Journal on Computing, 10(4):713–717, 1981. 75

[IS24]

Misha Ivkov and Tselil Schramm. Semidefinite programs simulate approximate message passing robustly. In Proceedings of the 56th Annual ACM Symposium on Theory of Computing (STOC 2024), pages 348–357. ACM, 2024. 11

[Jan97]

Svante Janson. Gaussian Hilbert spaces, volume 129 of Cambridge Tracts in Mathematics. Cambridge University Press, Cambridge, 1997. 15

[JM12]

Adel Javanmard and Andrea Montanari. Subsampling at information theoretically optimal rates. In Proceedings of the International Symposium on Information Theory (ISIT 2012), pages 2431–2435. IEEE, 2012. 6

[JM13]

Adel Javanmard and Andrea Montanari. State evolution for general approximate message passing algorithms, with applications to spatial coupling. Information and Inference: A Journal of the IMA, 2(2):115–144, 2013. 3, 10, 11, 25, 62, 65

[JP25]

Chris Jones and Lucas Pesenti. Fourier analysis of iterative algorithms. In 52nd International Colloquium on Automata, Languages, and Programming (ICALP 2025), pages 102–1. Schloss Dagstuhl–LeibnizZentrum für Informatik, 2025. 11, 24, 26, 92

[KMW24]

Dmitriy Kunisky, Cristopher Moore, and Alexander S. Wein. Tensor cumulants for statistical inference on invariant distributions. In 65th Annual Symposium on Foundations of Computer Science (FOCS 2024), pages 1007–1026, 2024. 80

[LW22]

Gen Li and Yuting Wei. A non-asymptotic framework for approximate message passing in spiked models. arXiv preprint arXiv:2208.03313, 2022. 11

[LWF25]

Max Lovig, Tianhao Wang, and Zhou Fan. On universality of non-separable approximate message passing algorithms. arXiv preprint arXiv:2506.23010, 2025. 3

[Mal20]

Camille Male. Traffic distributions and independence: permutation invariant random matrices and the three notions of independence, volume 267:1300. American mathematical society, 2020. 5, 7, 23, 24, 26, 27, 65

[MFC+ 19] Antoine Maillard, Laura Foini, Alejandro Lage Castellanos, Florent Krzakala, Marc Mézard, and Lenka Zdeborová. High-temperature expansions and message passing algorithms. Journal of Statistical Mechanics: Theory and Experiment, 2019(11):113301, 2019. 7, 12, 19 [MM09]

Marc Mézard and Andrea Montanari. Information, Physics, and Computation. Oxford University Press, 2009. 3

[Mon12]

Andrea Montanari. Graphical models concepts in compressed sensing. In Yonina C. Eldar and Gitta Kutyniok, editors, Compressed Sensing: Theory and Applications, pages 394–438. Cambridge University Press, 2012. 1, 11

[Mon19]

Andrea Montanari. Optimization of the Sherrington-Kirkpatrick Hamiltonian. In 60th IEEE Annual Symposium on Foundations of Computer Science (FOCS 2019), pages 1417–1433. IEEE, 2019. 2

[MP17]

Junjie Ma and Li Ping. Orthogonal AMP. IEEE Access, 5:2020–2033, 2017. 12

[MPR94a]

Enzo Marinari, Giorgio Parisi, and Felix Ritort. Replica field theory for deterministic models: I. binary sequences with low autocorrelation. Journal of Physics A: Mathematical and General, 27(23):7615, 1994. 6

[MPR94b]

Enzo Marinari, Giorgio Parisi, and Felix Ritort. Replica field theory for deterministic models. II. A nonrandom spin glass with glassy behaviour. Journal of Physics A: Mathematical and General, 27(23):7647, 1994. 6

[MPV87]

Marc Mézard, Giorgio Parisi, and Miguel Angel Virasoro. Spin glass theory and beyond: An Introduction to the Replica Method and Its Applications, volume 9. World Scientific, 1987. 3

[MR15]

Andrea Montanari and Emile Richard. Non-negative principal component analysis: message passing algorithms and sharp asymptotics. IEEE Transactions on Information Theory, 62(3):1458–1484, 2015. 2

[MS12]

James A. Mingo and Roland Speicher. Sharp bounds for sums associated to graphs of matrices. Journal of Functional Analysis, 262(5):2272–2288, 2012. 32

[MW24]

Andrea Montanari and Yuchen Wu. Statistically optimal first-order algorithms: a proof via orthogonalization. Information and Inference: A Journal of the IMA, 13, 2024. 1

[MW25]

Andrea Montanari and Alexander S. Wein. Equivalence of approximate message passing and low-degree polynomials in rank-one matrix estimation. Probability Theory and Related Fields, 191:181–233, 2025. 2, 11

[NS06]

Alexandru Nica and Roland Speicher. Lectures on the combinatorics of free probability, volume 13. Cambridge University Press, 2006. 77

[OÇW16]

Manfred Opper, Burak Çakmak, and Ole Winther. A theory of solving TAP equations for Ising models with general invariant random matrices. Journal of Physics A: Mathematical and Theoretical, 49(11):114002, 2016. 12

[Pes26]

Lucas Pesenti. Algorithms beyond the union bound: polynomial optimization and discrepancy theory. PhD thesis, Bocconi University, 2026. 2

[Pet82]

L.C. Petersen. On the relation between the multidimensional moment problem and the one-dimensional moment problem. Mathematica Scandinavica, 51:361–366, 1982. 86

[PP95]

Giorgio Parisi and Marc Potters. Mean-field equations for spin models with orthogonal interaction matrices. Journal of Physics A: Mathematical and General, 28(18):5267, 1995. 6, 7

[Ran11]

Sundeep Rangan. Generalized approximate message passing for estimation with random linear mixing. In IEEE International Symposium on Information Theory (ISIT 2011), pages 2168–2172. IEEE, 2011. 1, 10, 62

[RFSK16]

Sundeep Rangan, Alyson K. Fletcher, Philip Schniter, and Ulugbek S. Kamilov. Inference for generalized linear models via alternating directions and Bethe free energy minimization. IEEE Transactions on Information Theory, 63(1):676–697, 2016. 1

[Rob39]

Herbert Ellis Robbins. A theorem on graphs, with an application to a problem of traffic control. The American Mathematical Monthly, 46(5):281–283, 1939. 33

[Rot64]

Gian-Carlo Rota. On the foundations of combinatorial theory: I. Theory of Möbius functions. Probability Theory and Related Fields, 2(4):340–368, 1964. 18

[RSF19]

Sundeep Rangan, Philip Schniter, and Alyson K. Fletcher. Vector approximate message passing. IEEE Transactions on Information Theory, 65(10):6664–6684, 2019. 10, 12

[RSFS19]

Sundeep Rangan, Philip Schniter, Alyson K. Fletcher, and Subrata Sarkar. On the convergence of approximate message passing with arbitrary matrices. IEEE Transactions on Information Theory, 65(9):5339–5351, 2019. 6

[RV18]

Cynthia Rush and Ramji Venkataramanan. Finite sample analysis of approximate message passing algorithms. IEEE Transactions on Information Theory, 64(11):7264–7286, 2018. 11

[Sch20]

Philip Schniter. A simple derivation of AMP and its state evolution via first-order cancellation. IEEE Transactions on Signal Processing, 68:4283–4292, 2020. 6, 12

[SSS19]

Phil Saad, Stephen H. Shenker, and Douglas Stanford. JT gravity as a matrix integral. arXiv preprint arXiv:1903.11115, 2019. 70

[Tak19]

Keigo Takeuchi. Rigorous dynamics of expectation-propagation-based signal recovery from unitarily invariant measurements. IEEE Transactions on Information Theory, 66(1):368–386, 2019. 12

[tH74]

Gerard ’t Hooft. A planar diagram theory for strong interactions. Nuclear Physics B, 72(3):461–473, 1974. 70, 73

[VSR+ 15]

Jeremy Vila, Philip Schniter, Sundeep Rangan, Florent Krzakala, and Lenka Zdeborová. Adaptive damping and mean removal for the generalized approximate message passing algorithm. In IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP 2015), pages 2021–2025. IEEE, 2015. 3

[WZF22]

Tianhao Wang, Xinyi Zhong, and Zhou Fan. Universality of approximate message passing algorithms and tensor networks. Annals of Applied Probability, 34:3943–3994, 2022. 3, 6, 7, 10, 11, 12, 28, 29

[ZK16]

Lenka Zdeborová and Florent Krzakala. Statistical physics of inference: Thresholds and algorithms. Advances in Physics, 65(5):453–552, 2016. 3

[ZWF24]

Xinyi Zhong, Tianhao Wang, and Zhou Fan. Approximate message passing for orthogonally invariant ensembles: Multivariate non-linearities and spectral initialization. Information and Inference: A Journal of the IMA, 13(3):iaae024, 2024. 3, 29

Traffic Distributions via Feynman Diagrams

One of our motivations is to connect graph polynomials with the celebrated Feynman diagram technique from physics. In quantum field theory, Feynman diagram expansion is used to reduce matrix integrals into graphical calculations. We show in this section that this method can (heuristically) derive the traffic distribution of orthogonally invariant distributions (Theorem 4.2). The matrix model that we consider in this section is specified by a potential function V : R → R, and has partition function Z Z :=

dA e− 2 Tr V (A) ,

(47)

where M := Rn×n sym is the space of symmetric n×n matrices. Equivalently, this is the partition function of the random matrix A ∈ M sampled from the probability measure µV (A) ∝ exp − n2 Tr V (A) , which is a special case of an orthogonally invariant distribution (Section 4.2). In physics, matrix integrals such as Eq. (47) are viewed as a 0-dimensional theory: the variable is a matrix, and the partition function is a finite-dimensional integral rather than a functional integral over fields on space-time. The large-n expansion of such integrals is organized into diagrammatic contributions indexed by Feynman diagrams. 1. In the limit n → ∞, only planar diagrams contribute at leading order, an observation going back to foundational work of ’t Hooft [tH74, BIPZ78]. Related planarity phenomena also appear in mathematics, for example in the connections between large random matrices and non-crossing pairings. 2. In special scaling limits of the potential with n, the Feynman diagram expansion can be interpreted in terms of physical theories such as 2D gravity and certain string-theoretic models [DFGZJ95, CGH+ 17, SSS19]. The combinatorial approach in this paper fits naturally into this perspective. First, our results are formulated in the large-n limit, and the dominant combinatorial objects in that limit are planar, as in the ’t Hooft limit. Second, we show that our w- and z-polynomials are planar dual to the Feynman diagrams traditionally used in physics. Third, while the Feynman diagram method is based on perturbative expansion around the GOE potential V (x) = x2 /2, our rigorous results Theorems 4.2 and 6.2 still remain valid beyond the radius of convergence for perturbative methods. We present in this section the traditional approach for computing Eq. (47) based on Feynman diagrams. The argument is “combinatorially rigorous” (true at the level of generating functions), but not sufficient to rigorously derive the probabilistic conclusions.

A.1

Calculation of the free energy

For now, we restrict to the case where the potential in Eq. (47) is V (A) = 12 A2 + g4 A4 , where the coupling constant g measures the strength of the quartic interaction in the model. Such potentials appear in string theory, statistical physics (the λϕ4 theory), and the theory of integrable systems. The quartic term g4 Tr A4 can be viewed as a correction term to the GOE model, for which R ZGOE = M dA exp −n Tr A2 /4 . The idea of the Feynman diagram technique is to perturbatively expand this correction term, reducing to a problem on Gaussian variables. We illustrate this by computing the free energy of the 70

quartic model, namely the quantity ln Z (this example can be found in physics textbooks). For an observable quantity O, we write ⟨O⟩ := EA∼µV [O], and ⟨O⟩GOE := EA∼GOE [O]. We have n gn 4 dA exp − Tr A2 − Tr A 4 8 M gn = ZGOE · exp − Tr A4 . 8 GOE

A simple calculation shows that ZGOE = 2 2 integrate term-by-term:

exp −

gn 4 Tr A 8

2π n

n(n+1)

= GOE

. We Taylor expand the remaining part and

∞ X 1 gn s D s=0

−

Tr A4

s E GOE

(48)

The quantities ⟨Tr A4 ⟩GOE on the right-hand side are expectations over Gaussian random variables, and can be computed by Wick’s lemma (Lemma 2.8) to be a sum over all Wick contractions between the variables (in graph-theoretic terms, a sum over all perfect matchings). The propagator for a single contraction with a GOE matrix is the covariance of the Gaussians, 1 1 ⟨A[i, j]A[k, ℓ]⟩GOE = δik δjℓ + δiℓ δjk , n n

where δij :=

 1

if i = j

0

otherwise

(49)

A Feynman diagram represents a combinatorial type of Wick contractions. In the graphical notations of this paper, we would visualize each Tr A4 as a square, with Wick contraction having the effect of gluing together edges of the squares. The ’t Hooft double line notation, which is more common in physics, represents each Tr A4 as a vertex with four incident double edges. These representations are dual to each other (in the sense of planar duality); see Fig. 5 for comparison.

(b) One of the Wick contractions appearing in 2 ⟨Tr A4 ⟩GOE . The edges of the squares are glued together according to the matching to make a “pillow”.

4 2

(a) Tr A represented as two squares in our notation, compared to the ’t Hooft double line notation.

Figure 5: Our Feynman diagram notation vs. the ’t Hooft double line notation. The delta functions in the propagator enforce that the vertices of the squares have a consistent index i when the edges of the squares are glued together. Note that the propagator in Eq. (49) for 71

the GOE model allows A[i, j], A[k, ℓ] to be glued in either orientation (in contrast to the Gaussian Unitary Ensemble which would only have one term). Therefore, we define a Feynman diagram for the GOE to be an oriented perfect matching between the edges of the squares. For each Feynman diagram γ, the contribution of γ to Eq. (48) is: (i) a factor n per vertex of γ, since each vertex holds an index from [n] which is summed over in Tr A4 . (ii) a factor n1 per paired edge of γ from the propagator, Eq. (49). 1 (iii) a factor − gn 8 per square face of γ from Eq. (48). There is also an overall factor of |F (γ)|! where |F (γ)| equals the number of square faces in γ.

For example, the s = 1 term in Eq. (48) is

gn − 8

g · ⟨Tr A ⟩GOE = − 8

· 2 · n2 + 5 · n + 5 .

The Feynman diagrams are enumerated in Fig. 6.

Figure 6: The 12 Feynman diagrams associated to ⟨Tr A4 ⟩GOE . The two red edges are matched in the orientation specified by the arrows, and similarly for the blue edges. Gluing either of the left two diagrams results in a “taco”. Gluing the remaining diagrams results in degenerate polyhedra. After gluing, the “tacos” have 3 vertices, the middle diagrams have 2 vertices, and the right diagrams have 1 vertex, respectively.

For a given Feynman diagram γ, the total factor of n is |V (γ)| − |E(γ)| + |F (γ)| =: χ(γ) which is the Euler characteristic of the polyhedron γ. In total, we obtain a Feynman diagram expansion for the partition function, 1 g Z = ZGOE − |F (γ)|! 8 γ∈Γ

|F (γ)|

nχ(γ) ,

(50)

where Γ is the set of Feynman diagrams, the set of polyhedra built from square faces. Formally, Γ = ⊔s≥0 Γs , where Γs is the set of oriented perfect matchings between the edges of s squares. Taking the logarithm has the effect of restricting the summation to connected Feynman diagrams; this is the linked cluster theorem in quantum field theory [Eti24, Section 3.5]. We obtain: ln

Z ZGOE

1 g = − |F (γ)|! 8 γ∈Γ

where Γc ⊆ Γ are connected Feynman diagrams. 72

|F (γ)|

nχ(γ) ,

(51)

A.1.1

Asymptotic limit n → ∞

As n → ∞, Eq. (51) significantly simplifies because only the planar diagrams survive, i.e. polyhedra γ with “no holes,” which have the maximum possible Euler characteristic among connected graphs (χ(γ) = 2). This foundational observation goes back to ’t Hooft [tH74].15 We obtain, at first order, Z 1 ln n2 ZGOE

g 1 − |F (γ)|! 8

X γ∈Γc planar

|F (γ)|

+ O(n−2 ) .

(52)

In summary, the Feynman diagram method shows that the non-Gaussian component of the matrix model can be replaced by a generating function for graphs/surfaces which, in the n → ∞ limit, restricts to a generating function for planar graphs/surfaces with genus 0. This restriction leads to significant simplifications in diagrammatic calculations, in the same way as our cactus property and treelike property in the rest of the paper.

A.2

Calculation of general observables: Argument for Theorem 4.2

We now assume that the potential V (A) has the general form V (A) = 12 A2 + j≥3 cj Aj (arbitrary coefficients on A and A2 can be handled by centering and rescaling, respectively). We compute the traffic distribution of A, which consists of all Sn -invariant observables of A. The z-polynomials are a basis for these observables where, for each multigraph α, P

X 1 1 ⟨zα (A)⟩ = n n i:V (α),→[n]

* Y

A[i(u), i(v)]

{u,v}∈E(α)

Separating out the Gaussian part of the action from the higher-order interactions: 1 1 ⟨zα (A)⟩ = · n n

zα (A) exp − D

exp −

j j≥3 cj n Tr A

j j≥3 cj n Tr(A )

E GOE

GOE

The dual Feynman diagrams are built from polygons with j ≥ 3 sides, each of which comes with a factor of −cj , generalizing the situation from the previous section. A small generalization of the argument shows that the denominator is *



exp−

X j≥3

+



cj n Tr Aj 

= GOE

Y (−cj )|Fj (γ)|

 γ∈Γ

j≥3

|Fj (γ)|!

  nχ(γ) ,

(53)

where |Fj (γ)| denotes the number of j-sided faces in γ. The numerator can also be calculated diagrammatically. The Wick contractions go between a collection of polygons as well as the additional edges A[i, j] in zα (A) . Let Γ(α) be the set of 15 ’t Hooft studies unitarily invariant matrix models instead of orthogonally invariant ones. He takes a further step √ by sending g → 0 at the rate Θ(1/ n), i.e., fixing λ = g 2 n to be constant. His claim is that λ is the only parameter characterizing the physical properties of observables in the large-n limit, and by taking λ → ∞ one gains some intuition on the physical phenomena of strongly interacting particles. The limit g → 0 is less interesting for us, since the traffic distribution (hence also the spectrum) is asymptotically the same as the GUE whenever g = o(1).

Feynman diagrams, visualized as polyhedra built on a set of “boundary” edges α. Then



zα (A) exp−

+

cj n Tr Aj 

j≥3



= GOE

Y (−cj )|Fj (γ)|

 γ∈Γ(α)

j≥3

|Fj (γ)|!

  nχ(γ) · (1 − O(n−1 )) .

(54)

Note α is considered a boundary and does not count towards the faces Fj (γ). To enforce that the labels of zα (A) are injective, we remove from Γ(α) any matching which causes two vertices of α to have the same label. The factor 1 − O(n−1 ) arises because each vertex is summed over n − O(1) indices to maintain injectivity, instead of precisely n which we had previously. We obtain the final result by dividing Eq. (54) by Eq. (53). This has the effect of restricting to the set of connected Feynman diagrams Γc (α) ⊆ Γ(α) by an alternate version of the linked cluster theorem. The final Feynman diagram formula is: 



X Y (−cj )|Fj (γ)| 1   · nχ(γ)−1 · (1 − O(n−1 )) . ⟨zα (A)⟩ = n |F (γ)|! j γ∈Γ (α) j≥3

(55)

Remark A.1. An alternative approach to the calculation would be to first symmetrize zα (A) over O(n) which is the symmetry group of the matrix model (and is larger than Sn ), then to plug in the values of the O(n)-invariant observables (the trace polynomials). We find it simpler to Taylor expand the action directly. A.2.1

Asymptotic limit n → ∞

In the asymptotic limit n → ∞, the only diagrams in Eq. (55) with constant-order magnitude are those such that α is a cactus graph, and γ consists of polyhedra with genus 0 attached to each cycle of the cactus, which has χ(γ) = 1. We prove this combinatorially in the forthcoming Lemma A.2. The large-n combinatorial summation factors over the cycles of the cactus, since the genus-0 polyhedra on each cycle can be chosen independently. We obtain Y 1 ⟨zσ (A)⟩ + O(n−1 ) if α is a cactus 1 n ⟨zα (A)⟩ = σ∈cycles(α)  n  O(n−1 ) otherwise    

The limiting value n1 ⟨zσ (A)⟩ of the q-cycle diagram σ is equal to κq + O(n−1 ) by the moment/free cumulant relation Eq. (8). Thus, Eq. (55) recovers Theorem 4.2. Lemma A.2. Let α ∈ A be a connected multigraph, and let γ ∈ Γc (α). Then χ(γ) = 1 if and only if α is a cactus and γ consists of genus-0 polyhedra attached to each cycle of α. Proof. The only α for which Γc (α) is nonzero are the Eulerian α, since a polyhedron γ ∈ Γc (α) with boundary α must have a boundary which is a union of cycles. Therefore, it remains to argue about Eulerian graphs α. For Eulerian α, the γ ∈ Γc (α) which maximize the quantity χ(γ) = |V (γ)| − |E(γ)| + |F (γ)| are given by decomposing α into the maximum number of simple cycles, then attaching a genus 74

0 polyhedron to each cycle. This achieves χ(γ) = |V (α)| − |E(α)| + C where C is the number of cycles.16 We argue that:

|V (α)| − |E(α)| + C ≤ 1

(56)

for all Eulerian graphs α and this is achieved if and only if α is a cactus. Fix a maximum cycle partition of α. The C cycles are edge-disjoint so we can remove one edge from each one while maintaining that the graph is connected. Let α′ be the resulting graph. Then |V (α′ )| − |E(α′ )| = |V (α)| − |E(α)| + C. Since α′ is still connected we have |V (α′ )| − |E(α′ )| ≤ 1. The final inequality is an equality if and only if α′ is a tree and hence α is a cactus. This proves Eq. (56) and completes the lemma.

A.3

Mathematical comments on the Feynman diagram method

The Feynman diagram method is not mathematically rigorous, with (in our opinion) the main obstruction being that intermediate summations such as Eqs. (50), (51) and (54) are divergent. The Euler characteristic grows with the number of disconnected polyhedra, but the method proceeds anyway to divide out the disconnected polyhedra, which ultimately yields a convergent summation in Eq. (52) (for sufficiently small values of the coupling constant g ≥ 0). The Feynman diagram method is a perturbative expansion because it holds for sufficiently small perturbations of the GOE density, up to the radius of convergence of the Feynman diagram summations [EM03, GP13]. On the other hand, Theorem 4.2 holds beyond the radius of convergence of the Feynman diagram expansion in Eq. (55), so it would be impossible to prove the theorem using a perturbative expansion alone.

Traffic Distributions via Weingarten Calculus

We now present different tools and calculations for the traffic distributions of orthogonally invariant matrices based on the Weingarten formula for the moments of entries of Haar-random orthogonal matrices. These essentially follow the ideas of similar calculations by [CDM24], but use the version of the Weingarten formula for the orthogonal group, which we review below.

B.1

Weingarten formula for orthogonal matrices

For i = (i1 , . . . , ik ) and a perfect matching α ∈ Mperf ([k]), define δα (i) =

 1

if iu = iv for all {u, v} ∈ α,

0

otherwise.

The Weingarten calculus expresses the moments of the Haar measure on O(n) in terms of a certain “Weingarten function” Wn (α, β) on pairs of matchings. 16

Note that the computational problem of, given an Eulerian graph α, compute a partition of E(α) into the maximum number of cycles, is NP-hard [Hol81].

Lemma B.1 (Weingarten formula). Let Q ∼ O(n) be a Haar-random orthogonal matrix. There exists a function Wn : Mperf ([k])2 → R such that E

Q∼O(n)

[Q[i1 , j1 ] · · · Q[ik , jk ]] =

Wn (α, β)δα (i)δβ (j).

X α,β∈Mperf ([k])

See [CŚ06, Ban10] for an explicit definition of Wn (α, β). We will only be interested in asymptotics for k constant and n → ∞, for which the approximations below will suffice. When k is odd, Mperf ([k]) = ∅, so the right-hand side above is zero, and indeed the left-hand side is easily seen to be zero without invoking the Weingarten formula, because Q has the same law as −Q. So, the only interesting case is k even. In that case, we give Mperf ([k]) the structure of a metric space, where ∆(α, β) is defined as the minimum number of swap operations needed to reach β from α (a swap replaces pairs {a, b}, {c, d} with pairs {a, c}, {b, d}). It is easy to check that ∆ is a metric (indeed, it is the distance on a certain graph structure defined on Mperf ([k])). Further, write cyc(α, β) for the set of even cycles formed by the disjoint union of α and β. Then, it is easy to show the alternative characterization ∆(α, β) =

k − |cyc(α, β)| . 2

As a sanity check, |cyc(α, β)| ≤ k2 with equality achieved if and only if α = β, which is precisely the case ∆(α, β) = 0. For α, β ∈ Mperf ([k]), let P(α, β) be the set of geodesic paths from p to q in Mperf ([k]), i.e., of P sequences α = γ0 , γ1 , . . . , γt = β with γi ̸= γi+1 for all i = 0, . . . , t − 1 and with t−1 i=0 ∆(γi , γi+1 ) = := ∆(α, β). For such a path P = (γ0 , . . . , γt ), write |P | t. Then, we define µ(α, β) =

(−1)|P |

P ∈P(α,β)

This may be viewed as a Möbius function of the partially ordered set whose chains are geodesics from a given “base” matching p to each other matching. An explicit formula from [CŚ06] is µ(α, β) =

(−1)

|C| −1 2

C∈cyc(α,β)

|C| −1 Cat 2

(57)

where Cat(·) are the Catalan numbers. The key asymptotic for the Weingarten function for our purposes is then the following: Proposition B.2 ([CŚ06]). For a fixed k and α, β ∈ Mperf ([k]), as n → ∞ we have

Wn (α, β) = n−k+cyc(α,β) µ(α, β) + O(n−1 ) . Note that the maximum possible scaling of this quantity is n−k/2 , which corresponds to the fact that with high probability the entries of Q are all roughly of order n−1/2 .

B.2

Möbius inversion on non-crossing partitions

Recall that NC(k) is the partially ordered set of non-crossing partitions, i.e., those whose parts do not cross when drawn as a partition of vertices of the k-cycle. We review some standard properties of this partially ordered set; see, e.g., [NS06] for a standard reference. Each non-crossing partition π ∈ NC(k) has a natural dual partition, called the Kreweras complement and denoted K(π). On the cycle graph Ck , this may be viewed as the maximal non-crossing partition of the midpoints of the edges of Ck that does not cross the boundaries of π. Alternatively, one may view both partitions as placed on a single cycle graph of twice the size, C2k , on alternating sets of vertices. We show this viewpoint with an example in Fig. 7. The map K : NC(k) → NC(k) is easily checked to be an involution. 8′

1′

7′

2′

6′

3′ 6

4 5′

4′

Figure 7: An illustration of the Kreweras complement operation on non-crossing partitions. The parts of a partition π ∈ NC(8) are drawn in blue, and the parts of the Kreweras complement K(π) ∈ NC(8) in red. We give NC(k) the usual partial ordering of refinement of partitions, written π ⪯ ρ, using that a refinement of a non-crossing partition remains non-crossing. This partial ordering has a minimal element 0 ∈ NC(k), the partition where every block is a singleton, and a maximal element 1 ∈ NC(k), the partition with just one block. The Kreweras complement is an anti-isomorphism of this ordering: it is a bijection that reverses the ordering, i.e. K(π) ⪯ K(ρ) if and only if π ⪰ ρ. In particular, K(0) = 1 and K(1) = 0. The Möbius function for the NC(k) poset gives values µ(π, ρ) for each pair π ⪯ ρ. The Kreweras complement interacts with the Möbius function in the following way that will be crucial for our purposes: µ(0, π) = µ(K(π), 1) . (58) Further, evaluations of the Möbius function as on the left-hand side may be expanded as products over the blocks of π, and the factors turn out to be the same as the combinatorial quantities appearing in Eq. (57); there is a combinatorial explanation for this coincidence but we will just need 77

to use that this indeed occurs: µ(0, π) =

(−1)|A|−1 Cat(|A| − 1) .

A∈π

Note that, applying Möbius inversion to Eq. (7), we obtain an explicit formula for the free cumulants in terms of the moments, as mentioned earlier in the main text: if mk are the moments of a probability measure, then the free cumulants κk are κk =

µ(π, 1k )

(59)

m|A| .

A∈π

π∈NC(k)

B.3

Tracial moments concentration

The result of [CDM24] also assumes the following formula for the joint moments of the various trace powers of a matrix, that we also use in our proof. We show that it follows from our assumptions. Lemma B.3 (Tracial moments concentration). Let A = A(n) ∈ Rn×n sym be random matrices that 2 converge in tracial moments in L to some µ. Then for any cycle diagrams ρ1 , . . . , ρk , 

lim E 

n j=1

n→∞

(n)

Proof. Let us write Tq = Tq 

E



k Y 1

wρj (A) =

q∈Q

Tq  −

1 lim E wρj (A) . n j=1 n→∞

:= n1 Tr Aq . For any finite multiset of integers Q, we can expand

 Y

k Y

 Y

E Tq =

E

∅̸=Q′ ⊆Q

q∈Q



(Tq − E Tq )

q∈Q′

E Tq .

q∈Q\Q′

Our goal is now to show that each term in the sum over Q′ converges to 0 as n → ∞. Fix Q′ ⊆ Q such that Q′ ̸= ∅, and select an arbitrary element q0 ∈ Q′ . By Cauchy-Schwarz, we have 2

  E 



(Tq − E Tq ) ≤ E(Tq0 − E Tq0 )2 · E 

q∈Q′

 Y

(Tq − E Tq )2  .

(60)

q∈Q′ \{q0 }

We know that E(Tq0 − E Tq0 )2 converges to 0 as n → ∞ by the L2 tracial moments convergence q/p assumption. For the remaining product of expectations from Eq. (60), we apply the bound Tq2 ≤ T2p ′′ ′ for all q ≤ p to get: for all Q ⊆ Q , Y

Tq2 ≤ T2 P

q∈Q′′

q q∈Q′′

Therefore, all terms in the expansion of Eq. (60) can be bounded by products of terms of the form E Tq for q ∈ N. These are all bounded as n → ∞, since convergence in L2 also implies convergence in expectation. Together, we deduce 

lim E 

n→∞

 Y

Tq  =

q∈Q

Y q∈Q

lim E Tq ,

n→∞

which is equivalent to the desired statement. Remark B.4. This property is a statement about concentration of the tracial moments. For an example where it does not hold, one can take A(n) = aIn for a ∼ Unif({±1}), in which case limn→∞ E n1 Tr(A) = limn→∞ E n1 Tr A3 = 0, while limn→∞ E[ n1 Tr(A) · n1 Tr A3 ] = 1. We will further show below that analogous formulas hold for joint moments of elements of the wand z-bases of polynomials, not just the cycle diagrams.

B.4

Traffic distribution of orthogonally invariant matrices

We now prove Theorem 4.2 by computing the traffic distribution of an orthogonally invariant matrix A, which we recall consists of the limits of expressions of the form n1 Ezα (A) for α ∈ A0 . First, for a graph α = (V (α), E(α)), define HE = HE(α) to be the set of half-edges in α, a set of size |HE| = 2|E| which may be identified with pairs (v, {v, w}) for each choice of v ∈ V and e ∈ Mperf (HE), {v, w} ∈ E(α). Then, to α itself is associated a distinguished perfect matching α which matches each pair (v, {v, w}) and (w, {v, w}) of half-edges that correspond to the same edge of α (this is the perfect matching that would realize α under the configuration model). We say that a matching β ∈ M(HE) is α-local if all of its matches are between half-edges of the form (v, e1 ), (v, e2 ), i.e., between pairs of half-edges associated to the same vertex (rather than e Let Loc(α) ⊆ Mperf (HE(α)) be the set of all α-local matchings. Note that the same edge) for α. Loc(α) ̸= ∅ if and only if α is Eulerian, i.e., if every vertex has even degree. e and the set Loc(α), which is minimized At the heart of the matter is the distance between α precisely by the cactus graphs α ∈ C:

Proposition B.5. For any graph α, not necessarily connected, all of whose connected components are Eulerian, we have e Loc(α)) = ∆(α,

min

β∈Loc(α)

e β) ≥ |V (α)| − |conn(α)| , ∆(α,

with equality if and only if every connected component of α is a cactus. Further, in that case, there is a unique β ∈ Loc(α) achieving equality, which is the (unique) such β that matches pairs of half-edges belonging to the same cycle in α. Proof. It suffices to consider α connected; the general case follows by considering each connected component separately. We may rewrite e Loc(α)) = ∆(α,

min

β∈Loc(α)

e β) = |E| − max |cyc(α, e β)| ∆(α, β∈Loc(α)

and therefore it suffices to show that, for all α-local matchings of half-edges β, we have (?)

e β)| ≤ |E| − |V | + 1 . |cyc(α, e and an α-local β is equivalently the number of cycles in The set of cycles in the disjoint union of α a cycle cover of α (i.e., a partition of its edges into cycles).

Figure 8: An illustration of the matchings involved in Proposition B.5 and the arguments afterwards. Gray regions represent vertices of a cactus graph α, three triangles joined at a vertex with one 4-cycle attached to one of those triangles at a different vertex. This graph has |V | = 10 and |E| = 13. e realizing the Black dots in those regions represent half-edges of α. In blue, we draw the matching α graph α; if the gray regions are each contracted to a point and only blue edges are retained, then the resulting graph is the cactus α. In red, we draw the unique α-local matching β that maximizes e β)| = |E| − |V | + 1 = 4; its being α-local corresponds to making matches only within the |cyc(α, gray regions. The bound is tight for cycles. Suppose C1 , . . . , Ck is a cycle cover of some connected multigraph α. Since α is connected, it is possible to order the Ci such that Ci+1 has a vertex in common with the union of C1 , . . . , Ci for each i = 1, . . . , k − 1. Adding each successive Ci then increases |E| − |V | + 1 by at least 1, so the bound follows. If the bound is tight, then in the above ordering Ci+1 must have exactly one vertex in common with the union of C1 , . . . , Ci , and thus is a cactus. In that case, there is only one cycle cover, and thus the minimizer β is unique and must be as specified in the statement. We now proceed to Theorem 4.2 by calculating the traffic distribution of a sequence of orthogonally invariant random matrices A = A(n) . We could view this as A = QDQ⊤ for a Haar-distributed orthogonal Q and some random diagonal D, but actually this will not be necessary. Instead, let us take a perspective similar to the calculations in, for instance, [KMW24], which we believe is useful in general. Our idea will be to average the zα (A) over a random rotation Q drawn independently of A. Regardless of the structure of A, this defines another family of polynomials: z̄α (A) := E zα (QAQ⊤ ) . Q

If Q is drawn from Haar measure, then the z̄α will be orthogonally invariant polynomials, a greater symmetry than permutation invariance of zα . In particular, since the invariants of matrices under the O(n) action are generated by traces of matrix powers, z̄α will be a polynomial in these. Proof of Theorem 4.2. Let Q ∈ O(n) be Haar-distributed and independent of A, and let α = (V, E) be a graph. As above, write HE = HE(α) for the set of half-edges. We start by directly expanding

the averaged polynomial z̄ introduced above: z̄α (A) = E zα (QAQ⊤ ) Q

=E Q

(QAQ⊤ )[i(v), i(w)]

i:V ,→[n] {v,w}∈E



=E Q



n X

Q[i(v), j1 ]A[j1 , j2 ]Q[i(w), j2 ]

 j1 ,j2 =1

i:V ,→[n] {v,w}∈E

Q[i(v), j(v, {v, w})]Q[i(w), j(w, {v, w})]A[j(v, {v, w}), j(w, {v, w})]

i:V ,→[n] {v,w}∈E j:HE→[n]



 Y

E

i:V ,→[n] j:HE→[n]

Q[i(v), j(v, {v, w})]Q[i(w), j(w, {v, w})]

{v,w}∈E

A[j(v, {v, w}), j(w, {v, w})]

{v,w}∈E

Here, we may use the Weingarten calculus, viewing the matchings involved as matchings of half-edges, provided that we view i : V → [n] as extended to i′ : HE → [n] by i′ (v, e) := i(v), i.e., labelling a half-edge by the vertex involved. This gives: =

W (β, γ)δβ (i′ )δγ (j)



 X

 j:HE→[n] β,γ∈Mperf (HE)

A[j(v, {v, w}), j(w, {v, w})]

{v,w}∈E

i:V ,→[n] β,γ∈Mperf (HE) j:HE→[n]

δβ (i′ ) W (β, γ)δγ (j)

A[j(v, {v, w}), j(w, {v, w})]

{v,w}∈E

i:V ,→[n]

The summation over i is zero unless β ∈ Loc(α), and in that case each choice of i contributes 1, for a total of n|V | (1 + O(n−1 )). So, we have

= 1+O

1 n



n|V |

 X

 γ∈Mperf (HE)

W (β, γ)

β∈Loc(α)

δγ (j)

A[j(v, {v, w}), j(w, {v, w})]

{v,w}∈E

j: HE→[n]

The remaining summation may be grouped into summations over the cycles in the disjoint union of e which gives γ and α, 1 = 1+O n



|V |

 X

 γ∈Mperf (HE)

W (β, γ)

β∈Loc(α)

Y C∈cyc(α e,γ)

Tr A

|C| 2

Now, we may use the asymptotic formula in Proposition B.2 and normalize the traces to get

= 1+O

1 n

 X

1 1+O n

β∈Loc(α)

= 1+O

1 n

 X

Tr A

µ(β, γ)

|C| 2

C∈cyc(α e,γ)

1 1+O n

1 = 1+O n

n|V |

−|HE|+cyc(β,γ)

γ∈Mperf (HE)

β∈Loc(α)







X γ∈Mperf (HE)



n|V |



n|V |

−|HE|+cyc(β,γ)+cyc(α e,γ)

|C| 1 Tr A 2 n

µ(β, γ)

C∈cyc(α e,γ)

X γ∈Mperf (HE)

1 1+O n

 β∈Loc(α)



n−∆(β,γ)−∆(αe,γ) µ(β, γ)

|C| 1 Tr A 2 . n

Y C∈cyc(α e,γ)

Let us pause to notice that we have achieved our initial goal, expressing the orthogonally invariant polynomial z̄α (A) as a polynomial in traces of powers of A. We now use that, if A was orthogonally invariant to begin with, then E zα (A) = E z̄α (A) A

and continue to determine the right-hand side as n → ∞. By the triangle inequality, we have e γ) ≥ ∆(β, α) e ≥ ∆(α, e Loc(α)) ≥ |V | − 1 . ∆(β, γ) + ∆(α,

Therefore, under our assumptions, all terms are negligible as n → ∞ except for those where equality is achieved throughout above. By Proposition B.5, we then find that if α is not a cactus then lim

E zα (A) = 0 .

n→∞ n A

So, suppose that α is a cactus. Then, using the factorization property (Lemma B.3), we have in the limit that 1 EA zα (A) = n→∞ n lim

µ(β, γ)

β∈Loc(α) γ∈Mperf (HE) ∆(β,α e)=|V |−1 ∆(β,γ)+∆(γ,α e)=∆(β,α e)

Y C∈cyc(α e,γ)

m|C|/2

where mk are the spectral moments. Letting η be the α-local matching of half-edges belonging to the same cycle around each vertex, by the uniqueness clause of Proposition B.5 we further have that only the term β = η contributes, giving =

µ(η, γ)

γ∈Mperf (HE) ∆(η,γ)+∆(η,α e)=∆(η,α e)

m|C|/2

C∈cyc(α e,γ)

e = k, and, rewriting the condition on γ in terms Suppose there are k cycles in α. Then, |cyc(η, α)| of cycle counts and using the explicit formula for the Möbius function from Eq. (57), we have

(−1)

|C| −1 2

γ∈Mperf (HE) C∈cyc(η,γ) |cyc(η,γ)|+|cyc(γ,α e)|=|E|+k

|C| Cat −1 2

m|C|/2

C∈cyc(α e,γ)

Now, we use that all γ appearing in the sum must only match half-edges belonging to the same e both have this property also, the various sets of cycles above all form partitions cycle. Since η and α of the cycles in α. Thus, the entire sum factorizes over the cycles of α. Further, those γ that are in e corresponding to non-crossing partitions the sum have both the partitions of cyc(η, γ) and cyc(γ, α) of each cycle of α, and these two non-crossing partitions are Kreweras complements of one another. Putting together all these combinatorial observations, we find: 

 X

Y  C∈cyc(α)

(−1)|A|−1 Cat(|A| − 1) ·

m|B|  .

B∈π

π∈NC(|C|) A∈K(π)

Now we use Eq. (58) and Eq. (59) to complete the proof: 

 X

 C∈cyc(α)

µ(0|C| , K(π)) ·

= =

 X

 C∈cyc(α)

m|B| 

B∈π

π∈NC(|C|)

 Y

µ(π, 1|C| ) ·

m|B| 

B∈π

π∈NC(|C|)

κ|C| ,

C∈cyc(α)

where we have at last identified the free cumulants, completing the calculation. We also note that, by exactly the same argument but using the disconnected case of Proposition B.5, we may equally well calculate suitably normalized limits of the values of disconnected diagrams in the z-basis, which factorize over their connected components: Proposition B.6. Let A = A(n) ∈ Rn×n sym be a sequence of orthogonally invariant random matrices that converge in tracial moments in L2 to a probability measure µ. Let D denote their limiting traffic distribution, which exists by Theorem 4.2 and is given by the explicit formula stated there. Then, for all k ≥ 1 and α1 , . . . , αk ∈ A, k Y 1 lim E zα1 ⊔···⊔αk (A) = D(αi ) . n→∞ nk i=1

B.5

Concentration of traffic observables

As a corollary, we may also conclude that the traffic distribution is concentrated in the sense of Definition 3.15. This also extends [CDM24, Theorem 4.7] to orthogonally invariant distributions. Lemma B.7. Let A = A(n) be orthogonally invariant random matrices that converge in tracial moments in L2 to a probability measure µ. Then the traffic distribution concentrates for A (in the sense of Definition 3.15). Proof. Let k ≥ 2 and α1 , . . . , αk ∈ A. Then, by Lemma 3.17, it suffices to show the concentration property in the z-basis, namely that: lim E

n→∞

" k Y 1

n i=1

(?)

zαi (A) = lim

n→∞

k Y

1 E zαi (A) . n i=1

Note that, upon expanding the summations in the z-basis polynomials, we have zα1 (A) · · · zαk (A) = zα1 ⊔···⊔αk (A) + zβ1 (A) + · · · + zβM (A) , where α1 ⊔ · · · ⊔ αk is the disjoint union, while the βi are various graphs formed by identifying subsets of the vertices of this disjoint union according to different non-trivial partitions of the vertices, provided that no two vertices of the same αj are identified. In particular, all βi have at most k − 1 connected components. Therefore, by Proposition B.6, we have lim

n→∞ nk

zβi (A) = 0

for all i ∈ [M ]. Thus, k Y 1 1 lim E zαi (A) = lim k E[zα1 ⊔···⊔αk (A)] , n→∞ nk n→∞ n i=1

and the result then follows by Proposition B.6.

B.6

Traffic distribution of punctured orthogonally invariant matrices

Since the r-ROM plays an important role in our main results, let us sketch how similar calculations can give an explicit combinatorial description of its traffic distribution, and indeed that of the puncturing of any orthogonally invariant random matrices. Recall that in the main text we relied entirely on the implicit description of this traffic distribution via Lemma 3.14. The closed form we give below is completely explicit, but, being in terms of a rather complicated summation over matchings, seems less useful than the implicit one. We follow the notation from the proof in the previous section. Additionally, for a graph α and a matching β of the half-edges of α, we write loc(β) for the set of edges of β that go between half-edges of the same vertex of α, and nonloc(β) for the set of edges of β that go between half-edges e is the matching of half-edges of α corresponding to the of different vertices of α. Recall also that α edges actually in the graph α.

Theorem B.8. Let A = A(n) ∈ Rn×n sym be a sequence of orthogonally invariant random matrices that converges in tracial moments in L2 to a probability measure µ. Write mk for the kth moment of µ and Π = Π(n) = I − n1 11⊤ . Then, for all α ∈ A, 1 EA zα (ΠAΠ) = n→∞ n lim

(−1)|nonloc(β)|

X β∈Mperf (HE(α)) α⊔nonloc(β) is a cactus

µ(β, γ)

γ∈Mperf (HE(α)) ∆(β,γ)+∆(γ,α e)=∆(β,α e)

m|C| .

C∈cyc(α e,γ)

Proof. Following the same calculations as in the proof of Theorem 4.2 above but now applied to ΠQAQ⊤ Π, we instead find: 1 EQ,A zα (ΠQAQ⊤ Π) n X 1 = EA Wn (β, γ) · n β,γ∈M (HE) perf

Tr A|C| · zG(β) (Π)

C∈cyc(α e,γ)

where G(β) denotes the graph formed by “wiring together” the matching of half-edges β (so e = α). Note that here if we replaced Π by I, we would get zG(β) (I) = that, for example, G(α) |V | −1 1β∈Loc(α) n (1 + O(n )), compatible with the previous calculation in the proof of Theorem 4.2, and indeed the above is true for an arbitrary symmetric matrix Π, not only the particular projection we are concerned with. But, in our particular case, since Π is constant on the diagonal and on the off-diagonal, we have X 1 = Wn (β, γ) · EA n β,γ∈M (HE) perf

Tr A

|C|

·n

C∈cyc(α e,γ)

|V |

1 − n

|nonloc(β)|

1 1− n

|loc(β)|

and now by the same asymptotics as before, =

X β,γ∈Mperf (HE)

1+O

1 n

n−∆(αe,γ)−∆(β,γ)−|nonloc(β)|+|V |−1 µ(β, γ)(−1)|nonloc(β)|

m|C|

C∈cyc(α e,γ)

e of its half-edges, and any other We claim that, for any connected α realized by the matching α matching β of the half-edges of α, we have e β) + |nonloc(β)| ≥ |V | − 1 . ∆(α,

As before, this is equivalent to having e β)| ≤ |E| + |nonloc(β)| − |V | + 1 . |cyc(α,

Consider an ancillary graph α′ constructed by adding edges to α for each non-local match in β. This graph is still connected, by parity considerations it must be Eulerian, and it has a total of e β)| is now the size of a cycle cover of α′ , and the claim then follows |E| + |nonloc(β)| edges. |cyc(α, by the bounds from the proof of Proposition B.5 applied to α′ . We also again have by the triangle inequality that e γ) + ∆(β, γ) ≥ ∆(α, e β) . ∆(α,

Thus, all terms in the sum above are of at most constant order. Further, those of constant order are those where the exponent of n is zero, which are those where the above bound is tight. By the characterization in Proposition B.5, this is precisely when α′ as formed above is a cactus, and the stated result follows after rearranging.

Convergence of Stochastic Processes

In Section 6, we deal with convergence in distribution of stochastic processes indexed by countably infinite set, intended as weak convergence in the product topology. Equivalently, this means that every finite-dimensional marginal converges in distribution. Definition C.1. Let A be a countable set. For random variables (x(n) )n≥1 and x∞ taking values in RA , we say that x(n) converges in distribution to x∞ and write (d)

x(n) −→ x∞ if, for every k ≥ 1 and α1 , . . . , αk ∈ A, we have (d)

(n) ∞ ∞ (x(n) α1 , . . . , xαk ) −→ (Xα1 , . . . Xαk ) .

To show convergence in distribution, we will use the method of moments [Bil95, Theorems 29.4, 30.1, 30.2]. The following theorem follows from Carleman’s conditions on moment-determinacy of a distribution on R, combined with [Pet82]. Theorem C.2 (Method of moments). Let (x(n) )n≥1 be a sequence of stochastic processes indexed by a countable set A. Assume that 1. All joint moments converge: for any k ≥ 1 and α1 , . . . , αk ∈ A, the limit of the joint moments lim E

n→∞

" k Y

x(n) αi

(61)

i=1

exists. 2. All marginals are subexponential: for every α ∈ A, there exists Cα > 0 such that for all p ≥ 1,

lim E x(n) α

n→∞

≤ (Cα p)2p .

(62)

Then x(n) converges in distribution to the unique law on RA with moments given by Eq. (61). Lemma C.3 (Truncation). Let (xn )n≥1 and (yn )n≥1 be sequences of random variables such that 1. For any K > 0, conditionally on |xn | ≤ K, (yn )n≥1 converges in distribution. 2. (xn )n≥1 is tight, i.e., supn≥1 Pr(|xn | > K) −→ 0. K→∞

Then, (yn )n≥1 converges in distribution. Proof. First, we prove: 86

Claim C.4. (yn )n≥1 is tight. Proof. For any K, L > 0, we have Pr(|yn | > L) ≤ Pr(|yn | > L | |xn | ≤ K) + Pr(|xn | > K). Pick K large enough so that the second term is bounded by ε uniformly in n. (yn )n≥1 is tight conditionally on |xn | ≤ K, so there exists L > 0 large enough so that the first term is also bounded by ε uniformly in n. By Claim C.4 and Prokhorov’s theorem, it remains to show that every subsequence of (yn )n≥1 that converges in distribution, converges to the same limit. Fix f : R → R to be a bounded continuous function and ε > 0. Then, by the law of total expectations, for any n ≥ 1, |E f (yn ) − E [f (yn ) | |xn | ≤ K]| = Pr(|xn | > K) (E [f (yn ) | |xn | > K] − E [f (yn ) | |xn | ≤ K]) ≤ 2∥f ∥∞ Pr(|xn | > K) ≤ε by setting K = K(ε) to be a large enough constant (with the second assumption). By the first assumption, there exists N ≥ 1 such that for any n, m ≥ N , |E [f (yn ) | |xn | ≤ K] − E [f (ym ) | |xm | ≤ K]| ≤ ε . In turn, this implies |E f (yn ) − E f (ym )| ≤ 3ε by the triangle inequality, so (E f (yn ))n≥1 is a Cauchy sequence, so it converges as n → ∞. This implies that every weak subsequential limit of (yn )n≥1 converges to the same limit, which concludes the proof.

C.1

Connection with convergence of the empirical distribution

Let us also remark on certain details concerning modes of convergence that are important to the use and interpretation of Theorem 6.2. Recall that we “stack” the zα (A) for α ∈ A1 into a single vector with more complicated entries, zA1 (A) ∈ (RA1 )n . Using our notation from Section 1, we then sample a random coordinate of this vector, forming a further random countably infinite vector samp(zA1 (A)) ∈ RA1 . This contains the ith entry of each zα (A), for a single shared randomly chosen i ∼ Unif([n]). Define the infinite ∞ similarly. Theorem 6.2 states that: random vector ZA 1 (d)

∞ samp(zA1 (A(n) )) −−−→ ZA , 1

(63)

n→∞

By the Cramér-Wold theorem, this is equivalent to: for any bounded continuous function φ and any finitely supported vector of coefficients cα , 







n X X 1X lim E φ cα zα (A)[i] = E φ  cα Zα∞  . n→∞ A n A i=1 α∈A α∈A 1

Alternatively, we may also make sense of this statement in terms of empirical distributions, which are just the laws of the random variables samp(x) discussed above. Definition C.5 (Empirical distribution). For x ∈ Rn , we write ed(x) := n1 empirical distribution of the entries of x. 87

i=1 δx[i] for the

Then, ed(zA1 (A)) is a random probability measure on the space RA1 , and the random variable samp(zA1 (A(n) )) is a single draw from this random probability measure. Its law is a deterministic probability measure on the space RA1 , which is the expectation of the random measure ed(zA1 (A)) (if µ is a random measure, then its expectation takes values (Eµ)(A) = E[µ(A)]). Thus, the above Eq. (63) is further equivalent to the weak convergence of probability measures (w)

∞ E ed(zA1 (A(n) )) −−−→ Law(ZA ). 1 n→∞

Again by the Cramér-Wold theorem, this is equivalent to, for any finitely supported coefficient vector of cα , having 

E ed 

 X

 (w)

cα zα (A(n) ) −−−→ Law 

α∈A1

n→∞

 X

cα Zα∞  .

α∈A1

In particular, since the output xt of a GFOM can be viewed in the above way, we see that the empirical distributions of xt = xt (A) are related to the asymptotic states Xt∞ by (w)

E ed(xt (A(n) )) −−−→ Law(Xt∞ ) . n→∞

Thus our results, interpreted in terms of convergence of the random empirical distributions of GFOM iterates, give convergence of the expectations of random measures. Often it is desirable to prove stronger modes of convergence in such situations, by proving that not only do we have n 1X φ(xt (A(n) )[i]) = Eφ(Xt∞ ) , n→∞ A n i=1

lim E

but also that the random variable inside the expectation concentrates over the randomness in A. We do not pursue this here, because it would require introducing additional assumptions on the matrices A involved, which may vary from application to application. As the example discussed in Remark B.4 shows, this kind of concentration does not follow automatically from the convergence in expectation that we show. An instructive example is the argument in [BLM15], which uses similar proof techniques to ours, but, to show that the above kind of convergence also happens in L2 uses a trick involving the entrywise independence of the Wigner matrices they work with (see their Proposition 5). In our much more general setting, it seems reasonable to ask instead for the convergence in the definition of the traffic distribution in Eq. (2) to happen in a stronger mode such as L2 . We leave the exploration of such conditions and the determination of which random matrix distributions they hold for to future work.

Omitted Proofs

D.1

Combinatorial lemmas

We gather here lemmas involving only graph combinatorics.

Lemma D.1. For all σ, σ ′ ∈ C1 and A ∈ Rn×n sym , zσ (A) · zσ′ (A) − zσ⊕σ′ (A) ∈ span(zA1 \C1 (A)) , where σ ⊕ σ ′ ∈ C1 is the grafting of σ and σ ′ at the root. Proof. In the z-basis expansion of zσ (A) · zσ′ (A), we sum over all possible partial matchings of the vertices of σ and σ ′ . The empty matching contributes exactly zσ⊕σ′ (A). Any other matching that merges some vertices u ∈ V (σ) and v ∈ V (σ ′ ) creates 4 edge-disjoint paths between the root and the merged vertex. Merging additional vertices of σ and σ ′ can only increase the number of edge-disjoint paths, so the resulting graphs cannot be cactuses. Lemma D.2. For all σ ∈ C1 , α ∈ A1 \ C1 and A ∈ Rn×n sym , zσ (A) · zα (A) ∈ span(zA1 \C1 (A)) . Proof. The proof is similar to Lemma D.1. In this case, the graph corresponding to the empty matching is not a cactus because α is not. All other matchings create at least 3 edge-disjoint paths between the root and the merged vertex. Lemma D.3. For each α ∈ A1 \ T1 and β ∈ A1 , zα (A) · zβ (A) ∈ span(zA1 \T1 ) . Proof. The non-treelike diagrams A1 \ T1 can be characterized as: Claim D.4. Let α ∈ A1 . Then α ∈ A1 \ T1 if and only if one of the following holds: (i) there exists a bridge edge which does not have a path to the root using only bridge edges, (ii) or there exist a pair of vertices with three edge-disjoint paths between them. Proof of Claim D.4. It is clear that either structure forbids α from being treelike. Conversely, if there are at most two edge-disjoint paths between all pairs of vertices, then the bridge edges of α go between cactuses. Then condition (i) characterizes whether all bridge edges are connected to the root. Using the claim, if α has a structure of type (ii) then this is preserved in the product terms with any β. Suppose then that α has a structure of type (i) and call the bridge edge e. Note that both α, β are connected by definition of A1 . If no descendants of e intersect with β, then the type (i) structure is preserved. Conversely, if any descendant of e intersects with β, then we obtain a new path from the descendant to the root through β which is disjoint from the other edges of α. Edge e has at least one ancestor which is not a bridge edge, hence there were already two edge-disjoint paths containing this ancestor. Together with the new path we obtain a structure of type (ii). In all cases, the product terms remain in A1 \ T1 . Proof of Lemma 6.9. First, if P matches an internal vertex of a hanging cactus, then it creates three edge-disjoint paths from the root to that vertex. These paths cannot be eliminated by merging

other vertices, so τP cannot be a cactus. Therefore, we may assume without loss of generality that τ1 and τ2 contain no hanging cactuses. It is straightforward to check that any homeomorphic matching yields a cactus. We focus on the converse. Specifically, suppose that we are given a matching P between the vertices of τ1 and τ2 such that τP ∈ C1 . We prove P ∈ H(τ1 , τ2 ) by induction on |V (τ1 )| + |V (τ2 )|. For the base case, suppose that τ1 or τ2 has only one vertex. Then τP can be a cactus only if both τ1 and τ2 consist of a single vertex. For the inductive step, let u11 , . . . , u1k be the children of the root of τ1 , and let u21 , . . . , u2ℓ be the children of the root of τ2 . A necessary condition for τP to be a cactus is that k = ℓ (and this is also necessary for P to be a homeomorphic matching). Moreover, after reordering u21 , . . . , u2k if necessary, we may assume that for all i ∈ [k], u1i and u2i lie on the same cycle in τP , and that these form exactly k distinct cycles in τP incident to the root. For each i ∈ [k] and j ∈ {1, 2}, let Sij denote the non-root vertices of τj that are mapped under P to the same cycle of τP as u1i , u2i . Claim D.5. For every i ∈ [k], there is exactly one vertex vi1 ∈ Si1 and exactly one vertex vi2 ∈ Si2 that are mapped to the same vertex of τP . Proof. Since τ1 and τ2 are acyclic, creating a cycle in τP requires identifying two other vertices than the root. Conversely, identifying more than one pair of vertices would create three edge-disjoint paths to the root in τP , contradicting the fact that the latter is a cactus. Claim D.6. For each i ∈ [k] and j ∈ {1, 2}, every pair in P incident to a vertex in the subtree rooted at vij has its other endpoint in the subtree rooted at vi3−j . Proof. Suppose for contradiction that there is a pair of P between a vertex w1 in the subtree rooted at vi1 and a vertex w2 in the subtree rooted at vi2′ for some i′ ̸= i. Then in τP there are three edge-disjoint paths from the image of vi1 to the root: two lie on the cycle formed by Si1 ∪ Si2 , and the third is obtained by concatenating the path from vi1 to w1 with the path from w2 to the root. This contradicts the fact that τP is a cactus. By Claim D.6, we may apply the induction hypothesis for each i ∈ [k] to the subtree of τ1 rooted at vi1 and the subtree of τ2 rooted at vi2 . Thus, the restriction of P to these subtrees is a homeomorphic matching. In particular, vi1 and vi2 have the same degree. Claim D.7. Let i ∈ [k]. Then vi1 and vi2 are either both in the core of their respective trees or both outside of it. Moreover, for each j ∈ {1, 2}, no vertex in Sij \ {vij } lies in the core of τj . Proof. For the first part, since vi1 and vi2 have the same degree, they are either both in the core or both outside the core. For the second part, suppose for contradiction that some w ∈ Sij \ {vij } lies in the core of τj . Since w has degree greater than 2, its image in the cactus τP is an articulation vertex. Let ρ be a cycle of τP incident to w that is distinct from the cycle induced by Si1 ∪ Si2 . Then the two neighbors of w in ρ are images of vertices of τj . Since τj is acyclic, the cycle ρ must contain a vertex w′ that is the image of a vertex of τ3−j . But then τP contains three edge-disjoint paths from w to the root: two through the cycle induced by Si1 ∪ Si2 , and a third obtained by following ρ from w to w′ and then the path from w′ to the root. This contradicts the fact that τP is a cactus. 90

Let i ∈ [k]. For j ∈ {1, 2}, let wij be the first descendant of uji that lies in the core of τj . By Claim D.7, there are only two cases: 1. Either vij = wij for both j ∈ {1, 2}. In this case, there are no non-core vertices to match on the path from uji to vij , so the induced matching is empty (and hence trivially order-preserving). 2. Or vij = ̸ wij for both j ∈ {1, 2}. In this case, by induction, the matching between vij and wij is order-preserving. Matching vi1 to vi2 and adding the matching from vij to wij yields an order-preserving matching from uji to wij . By induction, the restriction of P induces an isomorphism between the cores of τ1 and τ2 within each subtree rooted at vij . Since there is no core vertex on the path from uji to vij by Claim D.7, these local isomorphisms extend to an isomorphism between the cores of τ1 and τ2 globally. This concludes the proof. Proof of Lemma 6.10. Given γ1 , . . . , γℓ ∈ G1 , we can expand ℓ Y

zγj =

j=1

zγP ,

where P ranges over all partitions of V (γ1 ) ∪ . . . ∪ V (γℓ ) such that all roots are in the same block, but no two vertices of the same γi are in the same block. Suppose that γP is treelike. Claim D.8. Every internal vertex of a hanging cactus forms a singleton block. Proof. Suppose for contradiction that an internal vertex u of a hanging cactus in γ1 lies in the same block as some vertex v of γ2 . Let u′ be the attachment vertex of the cycle containing u. In γP , there are three edge-disjoint paths between the images of u and u′ : two are inherited from γ1 , while the third is obtained by following the path in τ2 from v to the root and then the path in τ1 from the root to u′ . This contradicts Lemma D.3, since γP is assumed to be treelike. By Claim D.8, we may temporarily delete the hanging cactuses from γ1 , . . . , γℓ and then reattach them in γP ; this does not affect whether γP is treelike. Hence, we may assume without loss of generality that none of γ1 , . . . , γℓ contains a hanging cactus. Claim D.9. Let M be the graph on [ℓ] with an edge between i, j ∈ [ℓ] if there exist u ∈ V (γi ) and v ∈ V (γj ) that lie in the same block of P . Then M is a matching. Proof. Suppose for contradiction that M is not a matching. Then there exist non-root vertices u ∈ V (γ1 ), v, v ′ ∈ V (γ2 ), and w ∈ V (γ3 ) such that u and v (resp. w and v ′ ) lie in the same block of P . Let v ′′ be the lowest common ancestor of v and v ′ in γ2 . Since γ2 ∈ G1 , v ′′ is not the root of γ2 . In γP , there are three edge-disjoint paths from the image of v ′′ to the root: one is the inherited path from v ′′ to the root inside γ2 ; the second follows the path in γ2 from v ′′ to v and then the path in γ1 from u to the root; and the third follows the path in γ2 from v ′′ to v ′ and then the path in γ3 from w to the root. This contradicts Lemma D.3, since γP is treelike by assumption.

By Claim D.9, it follows that ℓ Y

zγj −

j=1

zL uv∈M

M ∈M(ℓ) Pu,v ∀uv∈M

γPu,v ⊕

L u∈M /

γu ∈ span(zA1 \T1 ) ,

where for each edge uv ∈ M , the sum over Pu,v ranges over all partial matchings between V (γu ) and V (γv ) that fix the roots. Finally, note that unless Pu,v is empty, γPu,v cannot be a treelike diagram that is not a cactus. Indeed, no vertices in the hanging cactuses can be matched; otherwise it would create three edgedisjoint paths. Moreover, if we match two tree vertices, that would create two edge-disjoint paths to the root, and thus would force the diagram to be a cactus. Since the grafting of non-treelike diagrams is again non-treelike, the only treelike contributions arise when each factor γPu,v is a cactus. By Lemma 6.9, this forces Pu,v to be a homeomorphic matching. Hence, ℓ Y j=1

zγj −

zL uv∈M

M ∈M(ℓ) Puv ∈H(γu ,γv ) ∀uv∈M

γPu,v ⊕

u∈M /

γu ∈ span(zA1 \T1 ) ,

as desired.

D.2

Handling empirical averages

To represent expressions involving empirical averages, we allow the coefficients in a diagram representation to be formal polynomials in the quantities {⟨zα (A)⟩ : α ∈ A1 }. Another approach would be to use disconnected diagrams, as in [JP25]. Lemma D.10. Assume that A = A(n) satisfies the assumptions of Theorem 6.2, and furthermore, the traffic distribution concentrates for A (Definition 3.15). Let x=

cα zα (A)

(64)

α∈A1

for some finitely supported coefficients (cα )α∈A1 which are polynomials cα ∈ R[V] with V := {⟨zα (A)⟩ : α ∈ A1 }. Then, X ∞ X := cα (E ZA ) · Zα∞ (65) 1 α∈A1

is the asymptotic state of x. Moreover, if xt is of the form Eq. (64) for any t ≥ 1 and Xt is correspondingly defined as in Eq. (65), then (Xt )t≥1 is the asymptotic state of (xt )t≥1 . Proof. For polynomial test functions, the convergence in Eq. (33) follows directly from the concentration of the traffic distribution. Moreover, Lemma 3.16 implies that n1 zα (A) converges in L2 to a deterministic limit for any α ∈ A1 . So we can combine Lemma 6.17 with Slutsky’s lemma to obtain that Eq. (33) also holds for bounded continuous functions.

D.3

Proof of Lemma 6.30

In this section, we prove Lemma 6.30. We assume throughout that H satisfies the assumptions of Theorem 6.29. We will prove that xt and ut have the same state evolution by relating them to the following intermediate iteration: y0 ∼ N (0, I) ,

yt = HΠft−1 (yt−1 ) −

t−1 X

bs,t · (Πfs (ys ))

∀t ≥ 1 ,

(66)

s=0

where bs,t is defined in Eq. (45). Unless specified otherwise, all expectations in this section are taken with respect to both H and y0 . Theorem 6.18 does not apply to yt because of the Gaussian initialization y0 ∼ N (0, I) (instead of y0 = 1). To analyze this initialization, we extend the class of diagrams to generalized diagrams, that is, graphs α = (V (α), E(α)) together with an additional label p(v) ∈ N assigned to each vertex. The z-polynomial associated with a graph α is zα (A, y0 ) :=

A[i(u), i(v)]

i:V (α),→[n] {u,v}∈E(α)

y0 [i(v)]p(v) .

v∈V (α)

The collection of generalized vector diagrams A1 (y0 ) is defined analogously. Definitions such as T1 , ∞ G1 , and = extend to generalized diagrams by simply ignoring the labels p(v). As in the proof of Theorem 6.28, one caveat is that yt cannot be directly expanded as a linear combination of connected generalized vector diagrams, because the iteration involves the scalar quantity ⟨ft (yt )⟩. We therefore proceed as in Lemma D.10, viewing the coefficients in the diagram expansion as formal polynomials in these variables whenever necessary. Our first observation is that taking expectation over y0 in the z-basis turns (up to a scaling factor) a generalized diagram α into the same diagram α where the labels are ignored. Lemma D.11. For any generalized scalar diagram α (not necessarily connected) and any H ∈ Rn×n sym , E zα (H, y0 ) =

 Q 

v∈V (α) (p(v) − 1)!!

zα (H)

0

if p(v) is even for every v ∈ V (α) otherwise

Proof. In the z-basis, all vertices are assigned distinct labels. Therefore, we may take the expectation over y0 separately at each vertex, since the coordinates of y0 are independent. For each vertex v, we have EZ∼N (0,1) Z p(v) = (p(v) − 1)!! if p(v) is even or 0 if p(v) is odd. Next, we describe structural properties of the labels p(v) appearing in the diagram expansion of the iterates of the AMP iteration Eq. (66). ∞P

∞P

Lemma D.12. We have yt = τ cτ zτ (H, y0 ) and ft (yt ) = τ c′τ zτ (H, y0 ), where cτ and c′τ are supported on (generalized) treelike diagrams τ such that, for all v ∈ V (τ ):

p(v) =

   1   

if v is a leaf vertex of τ

0 or 2

if v is in a hanging cactus .

otherwise

Leaves of treelike diagrams are defined after removing hanging cactuses. Proof. First, the proof of Lemma 6.23 still goes through with the nonlinearities gt (y) = ft (y) − ⟨ft (yt )⟩, after extending the coefficient field from R to the ring of formal polynomials in {⟨zα (A)⟩ : α ∈ A1 }. Therefore, we obtain t−1 ∞ X

yt =

Bs,t (Πfs (ys ))̸=1 .

(67)

s=0

We now argue by induction on t. The base case is f0 (y0 ) = y0 which is the singleton with p(v) = 1. Now, suppose that the claim holds for yt . The treelike diagrams appearing in ft (yt ) are obtained by considering all possible products of treelike diagrams γ1 , . . . , γℓ ∈ G1 ∪ C1 appearing in yt . By Lemma 6.10, each such product can be written as a sum over matchings among the γi , where each γi is either paired into a cactus or does not intersect any other γj . In the second case, the values p(v) within γi are unchanged. In the first case, the values p(v) at the leaves are updated from 1 to 2, while all other values p(v) within γi remain unchanged. Moreover, no non-trivial intersection between Bs,t and (Πfs (ys ))̸=1 can produce a treelike diagram. Hence, the decomposition of yt+1 given by Eq. (67), together with the induction hypothesis, shows that in every treelike diagram appearing in yt+1 , the condition on p(v) is inherited directly from the corresponding property of fs (ys ). This completes the induction. Lemma D.13. For any t ≥ 1 and any polynomial φ : Rt → R, lim

E ⟨φ(x1 , . . . , xt )⟩ = lim

n→∞ H,x0

E ⟨φ(y1 , . . . , yt )⟩ .

n→∞ H,y0

Proof of Lemma D.13. An iteration involving A can be reduced to one involving H by expanding Aft (yt ) = Hft (yt ) − ⟨ft (yt )⟩H1 − ⟨HΠft (yt )⟩1 .

(68)

Set δt := ⟨HΠft (yt )⟩ and mt := ⟨ft (yt )⟩. We first compare yt with the following modified iteration, which differs from xt only in the formula for the Onsager correction term: ỹ0 = y0 ,

ỹt = Aft−1 (ỹt−1 ) −

t−1 X

bs,t · (Πfs (ỹs )) ,

s=0

where bs,t is defined in Eq. (45). Claim D.14. For any t ∈ N, we have ỹt − yt =

ct,α (δ0 , . . . , δt−1 , m0 , . . . , mt−1 )zα (H, y0 ) ,

(69)

where the sum runs over finitely many generalized vector diagrams, and each ct,α is a polynomial in δ0 , . . . , δt−1 , m0 , . . . , mt−1 that is divisible by δs for some s ∈ {0, . . . , t − 1}. Proof of Claim D.14. We argue by induction on t. For t = 0, y0 − ỹ0 = 0, establishing the base case. Let t ≥ 1 and suppose that Eq. (69) holds for all s < t. First, one easily verifies from the

induction hypothesis that the same property Eq. (69) holds for ∆s := fs (ỹs ) − fs (ys ) for every s < t. By Eq. (68), we can then write Aft−1 (ỹt−1 ) − HΠft−1 (yt−1 ) = HΠ∆t−1 − δt−1 1 − ⟨HΠ∆t−1 ⟩1 , and each of the three terms on the right-hand side satisfies a decomposition of the form Eq. (69) by the induction hypothesis. Finally, the correction terms differ by t−1 X

bs,t · (Πfs (ỹs )) −

s=0

t−1 X

bs,t · (Πfs (ys )) =

s=0

t−1 X

bs,t · Π∆s ,

s=0

which again satisfies the property Eq. (69) by the induction hypothesis. Combining these observations, we conclude that ỹt − yt satisfies Eq. (69), completing the induction. Next, fix any polynomial φ : Rt → R. By Claim D.14, we have ⟨φ(y1 , . . . , yt )⟩ − ⟨φ(ỹ1 , . . . , ỹt )⟩ =

cα (δ0 , . . . , δt−1 , m0 , . . . , mt−1 )⟨zα (H, y0 )⟩ ,

(70)

where the sum runs over finitely many generalized scalar diagrams, and each cα is a polynomial in δ0 , . . . , δt−1 , m0 , . . . , mt−1 that is divisible by some δs . In the remainder of the proof, we show that each term on the right-hand side of Eq. (70) converges to 0 in expectation. The reason is that each coefficient cα contains a factor δs , and these quantities converge to 0 in L2 : L2

Claim D.15. ⟨HΠft (yt )⟩ −→ 0 for any t ≥ 1. Proof of Claim D.15. The claim is equivalent to the statement that n12 E⟨H1, Πft (yt )⟩2 converges to 0. This quantity can be expanded as a linear combination of terms of the form 1 E [zα (H, y0 )zβ (H, y0 )] n2 where α, β ∈ A both belong to the support of the expansion of ⟨H1, Πft (yt )⟩. As in the proof of Lemma B.7, 1 1 E [zα (H, y0 )zβ (H, y0 )] = 2 E [zα⊔β (H, y0 )] + o(1) . n2 n Indeed, each identification of vertices across the two copies yields a connected diagram whose expectation, after normalization by 1/n2 , converges to 0 by the existence of the traffic distribution. This holds for every realization of y0 , and therefore also after taking expectation over y0 . Taking expectation over y0 and using Lemma D.11, each term either vanishes or becomes a constant multiple of EH [zα⊔β (H)], where α ⊔ β is viewed as an ordinary scalar diagram obtained by ignoring the labels p(v). By Lemma B.7 and the strong cactus property, the only terms that contribute to the limit are those for which both α and β are cactuses. Viewing H1 as a rooted tree with one edge, the cactuses in the z-basis expansion of ⟨H1, Πft (yt )⟩ arise when the child of H1 is merged with a leaf of a diagram from Πft (yt ). By Lemma D.12, such leaves satisfy p(v) = 1. Applying Lemma D.11 once again, we find that each of these cactus terms has expectation 0 over y0 , which concludes the proof. After taking expectation over H and y0 , any monomial appearing on the right-hand side 95

of Eq. (70) has the following form for some pi , qi ∈ N: "

E δs

t−1 Y

δipi mqi i ⟨zα (H, z0 )⟩

≤ E δs2

1 2

i=0

· E

"t−1 Y 2p i

δi

#! 21 2 i m2q i ⟨zα (H, z0 )⟩

(71)

i=0

where the inequality follows from Cauchy-Schwarz. The first factor on the right-hand side of Eq. (71) converges to 0 as n → ∞ by Claim D.15. The second factor can be expanded in the z-basis as a finite linear combination of products of generalized z-diagrams. Taking expectation over y0 and using Lemma D.11, each such term either vanishes or becomes a constant multiple of a product of ordinary scalar z-diagrams. By Lemma B.7, the normalized expectation of each of these terms has a finite limit as n → ∞ and in particular, is uniformly bounded in n. Therefore, the second factor on the right-hand side of Eq. (71) is bounded, and hence the right-hand side of Eq. (70) converges to 0 in expectation. In summary, we have shown: lim E⟨φ(y1 , . . . , yt )⟩ − E⟨φ(ỹ1 , . . . , ỹt )⟩ = 0 . (72) n→∞

Finally, as in the proof of Theorem 6.28, we may use the traffic concentration property (Lemma B.7) Q to replace bs,t by κt−s s<r<t ⟨fr′ (ỹr )⟩ in the iteration for ỹt without affecting the asymptotic state. This yields lim E⟨φ(x1 , . . . , xt )⟩ − E⟨φ(ỹ1 , . . . , ỹt )⟩ = 0 . n→∞

Combining this with Eq. (72) completes the proof. Proof of Lemma 6.30. First, we can replace every occurrence of Πf0 (y0 ) in yt by f0 (y0 ) using the traffic concentration property, since ⟨f0 (y0 )⟩ = ⟨y0 ⟩, which converges to 0 as n → ∞. After this update, the iterates yt and ut have the same generalized diagram expansion as functions of their initializations y0 and u0 . Note that this expansion is formal in the variables ⟨zα (A)⟩ for α ∈ A1 (y0 ), because the puncturing operation introduces terms of the form ⟨ft (yt )⟩. First, by the strong cactus property and Lemma D.11, all non-cactus terms in the generalized diagram expansions of ⟨φ(y1 , . . . , yt )⟩ and ⟨φ(u1 , . . . , ut )⟩ converge to 0 in expectation. Second, using Lemma D.12 and extending the same argument one further step to φ, all cactus diagrams in the generalized diagram expansions of yt and ut satisfy p(v) ∈ {0, 2}, since they have no non-root leaves (the iterates for t ≥ 1 have no singleton component). Therefore, by Lemma D.11, the expectations of the cactus terms remain unchanged as n → ∞ if we replace y0 by u0 = 1. Combining these facts with the traffic concentration property for H (Lemma B.7) shows that ⟨φ(y1 , . . . , yt )⟩ − ⟨φ(u1 , . . . , ut )⟩ converges to 0 in expectation, as desired.

D.4

Proof of Lemma 6.31

In this section, we prove the auxiliary lemmas for block matrices. Definition D.16. Let α ∈ C1 be a cactus diagram. For a coloring χ : V (α) → [q] of the vertices of α with q colors, we say that χ is valid if for every cycle ρ = (u1 , . . . , uk , u1 ) ∈ cyc(α), there exist r, c ∈ [q] such that χ(ui ) = r when i is even and χ(ui ) = c when i is odd, with r = c if k is odd. We write χ(ρ) = {r, c} in this case.

Our main diagrammatic calculation for block models is the following, which gives the traffic distribution on each block: Lemma D.17. Let A be as in the setting of Lemma 6.31. Then for all α ∈ A1 and r ∈ [q]:

q n

E zσ (A)[i] −→

 X Y χ(ρ)  κ|ρ|    χ:V (α)→[q] ρ∈cyc(α)

if α ∈ C1

  0

if α ∈ A1 \ C1

χ valid n→∞   χ(root)=r

i∈[n] block(i)=r

Proof. We partition the sum defining zα (A) according to the block of each vertex, as in the proof of Proposition 4.6: zα (A) =

zαχ ((Ar,c )r,c∈[q] )

X χ:V (α)\{root}→[q]

where αχ is a diagram whose edges are colored by the matrices Ar,c . For a fixed r′ ∈ [q], we get q n

X i∈[n] block(i)=r′

E zα (A)[i] =

q n χ:V (α)\{root}→[q] X

E zαχ ((Ar,c )r,c∈[q] )[i] .

i∈[n] block(i)=r′

By the definition of traffic independence (Definition 4.5), the limit as n → ∞ exists for each term indexed by χ on the right-hand side. Hence, the limit of the left-hand side also exists. Arguing as in the proof of Proposition 4.6, we find that the limit is zero for all α ∈ A1 \ C1 . For cactus diagrams α ∈ C1 , asymptotic traffic independence and the strong factorizing cactus property of the individual blocks imply that the only nonzero contributions arise when every cycle of αχ is monochromatic, in the sense that it involves only a single matrix Ar,c . This happens if and only if χ is a valid coloring, in which case the corresponding term contributes asymptotically χ(ρ)

κ|ρ| ,

ρ∈cyc(σ)

as desired. Proof of Lemma 6.31. Let LC1 (r) denote the values from Lemma D.17: Lσ (r) :=

χ(ρ)

κ|ρ| .

χ:V (σ)→[q] ρ∈cyc(σ) χ valid χ(root)=r

We first prove that all joint moments of zC1 (A)[i] conditioned on block(i) = r converge to the

moments of the deterministic sequence ZC∞1 (r). For any σ1 , . . . , σk ∈ C1 , we have q n =

q n

E [zσ1 (A)[i] · · · zσk (A)[i]]

X i∈[n] block(i)=r

E zσ1 ⊕...⊕σk (A)[i] + o(1)

(by Lemmas D.1 and D.17)

i∈[n] block(i)=r

= Lσ1 ⊕...⊕σk (r) + o(1) =

k Y

Lσj (r) + o(1) .

(by Lemma D.17)

j=1

So it remains to prove that LC1 (r) satisfies the same recursion as ZC∞1 (r). First, one readily checks that Lsingleton (r) = 1, as in (i). Next, suppose that σ is rooted at a vertex of degree 2, and let ℓ and σ1 , . . . , σℓ−1 be as in (ii). Then, by decomposing according to the value of cycle containing the root, we have

Lσ (r) =

   ℓ ℓ  X Y Y    {r,c}    Lσk (r) Lσk (c) κℓ   k=2

c∈[q]

if ℓ is even

k=2 k even

k odd   ℓ  Y  {r,r}   Lσk (r) κℓ

if ℓ is odd

k=2

just like the recursion in (ii). Similarly, (iii) follows from the fact that the definition of Lσ (r) factorizes over graftings at the root. Together, this shows that LC1 (r) = ZC∞1 (r). Since the limit is deterministic, we have shown that conditionally on block(i) = r, zC1 (A)[i] converges to ZC∞1 (r) in L2 . Since block(i) ∼ Unif([q]), it follows that (block(i), zC1 (A)[i]) converges in distribution to (R, ZC∞1 (R)), where R ∼ Unif([q]).

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