arXiv CSopen access

Free Decompression with Algebraic Spectral Curves

2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

machine learning, deep learning, neural networks

Free Decompression with Algebraic Spectral Curves Siavash Ameli1 , Chris van der Heide2 , Liam Hodgkinson3 , and Michael W. Mahoney4 1

arXiv:2605.03634v1 [stat.ML] 5 May 2026

ICSI and Department of Statistics, University of California, Berkeley [email protected]

Dept. of Electrical and Electronic Engineering, University of Melbourne [email protected] 3

School of Mathematics and Statistics, University of Melbourne [email protected]

ICSI, LBNL, and Department of Statistics, University of California, Berkeley [email protected]

Abstract Tools from random matrix theory have become central to deep learning theory, using spectral information to provide mechanisms for modeling generalization, robustness, scaling, and failure modes. While often capable of modeling empirical behavior, practical computations are limited by matrix size, often imposing a restriction to models that are too small to be realistic. This motivates the inference of properties of larger models from the behavior of smaller ones. Free decompression (FD) is a recently proposed method for extrapolating spectral information across matrix sizes, but its utility is currently limited by strong assumptions that preclude its implementation on more realistic machine learning (ML) models. We use algebraic spectral curve theory to provide a general FD methodology for spectral densities whose Stieltjes transform satisfies an algebraic relation, a modeling assumption that is more likely to hold in practice. This recasts FD as an evolution along spectral curves which can be readily integrated. Our framework enables the expansion of spectral densities that have multiple or multi-modal bulks, that exist at multiple scales, and that contain atoms, all characteristic of real-world data and popular ML models. We demonstrate the efficacy of our framework on models of interest in modern ML, including Hessian and activation matrices associated with neural networks and large-scale diffusion models.

Introduction

At the core of scientific computing and much of modern machine learning (ML) lies the challenge of estimating the eigenvalues of high-dimensional Hermitian matrices. Such matrices, including kernels, Hessians, and graph representations, encode the intrinsic geometry and connectivity of the data and models built on them, rendering the pursuit of efficient spectral techniques a primary concern for both theory and practice. Studying eigenspectra has become a prominent approach to understanding performance and guiding training in deep learning [10, 20, 36, 53]. In many cases, the spectra of such matrices have non-trivial structure, often containing spikes, multiple multi-modal bulks, and heavy-tails [14, 25]. Conventional algorithms to extract eigenvalue information from these matrices have required that the data are able to be stored in memory, scratch space, or can at least be accessed as an implicit operator (via matrix-vector products). More recently, a new class of algorithms has emerged that is able to provide highly-accurate estimates of the eigenvalues (or summary functionals thereof [2]) of matrices, even without implicit or explicit access to the full matrix, i.e., of so-called impalpable matrices [1]. One such method, termed Free Decompression (FD), shows great promise as a tool for gaining access to the spectral distributions of such impalpable matrices. The central premise is that by appropriately sampling a small sub-matrix from the large impalpable matrix of interest, one can evolve a partial differential equation (PDE) in the Stieltjes transform of a spectral density in the decompression ratio to the desired matrix dimension. Despite demonstrating impressive accuracy when 1

Mass Transfer Between Atoms and Bulks

Evolution of Spectral Edges

25 K

1.0

I1 (τ )

0.8

100

I2 (τ )

0.4

(x∗ , τ∗ ) 23 K

0.2

I(τ ),

τ < τ∗

10−2

Bulk mass Atom mass

0.0 22 K

23 K

24 K

Matrix size n

25 K

ρ(λ, τ )

10−1

n(τ )

Mass

24 K

0.6

Spectral edge Cusp point

22 K −1

Figure 1: Evolution for increasing matrix sizes τ of atomic mass (left), density and spectral edges (right) for the free compound Poisson density from Section 4. A complexity dealt with by our method is that the bulk with support I(τ ) splits at the cusp point (x∗ , τ∗ ) to reveal two evolving supports I1 (τ ) and I2 (τ ). considering synthetic examples that naturally arise in random matrix theory (RMT) and free probability, the methods in [1] are highly restrictive, assuming the distribution of eigenvalues can be described by a single smooth connected bulk. This assumption is inevitably violated when challenged with more realistic examples, preventing widespread adoption of FD-based methods. Instead, this work considers a general FD methodology based on the polynomial method [39]. The polynomial method is a tool from RMT that, in theory, greatly simplifies many computations involving spectral distributions of large random matrices and their transforms. In particular, it assumes that the Stieltjes transform of the limiting spectral distribution of a class of random matrices is algebraic, i.e., it can be represented as a root of a (typically low-order) bivariate polynomial relation. This recasts the FD problem as evolution on a spectral curve. All canonical random matrix ensembles satisfy this assumption (see Remark H.1), and by increasing the degree of the polynomial, an extraordinarily broad class of spectral densities can be captured in this way. For example, Pennington et al. [36, 38] highlight that the Stieltjes transforms of certain Hessian approximations and neural network gram matrices are, respectively, solutions to cubic and quartic polynomial equations. Conveniently, it turns out that a wide range of common operations and transforms preserve this structure, including FD. This fact can be exploited to circumvent many of the challenges faced by the approach taken in Ameli et al. [1], including the problem of analytic continuation, which is well understood for polynomials. Unfortunately, there has been no implementation of the polynomial method for matrix computations to date, due to major practical challenges when fitting the underlying polynomial, and identifying the roots of interest. These issues must first be overcome to accomplish FD at scale. Contributions. The central contribution of this work is to introduce a new general FD methodology and practical tools for performing spectral computations for large impalpable matrices. We develop the first robust implementation of the polynomial method and apply it to obtain a state-of-the-art implementation of FD. This enables us to directly trace a number of spectral quantities of interest; including the edges of the distribution, locations of atoms and their weights, as well as spectral moments (see Figure 1). This is highly valuable in many applications: for positive-definite matrices, the ratio of the largest to the smallest eigenvalues constitutes the condition number cond(A) of the matrix A, a critical component for understanding rates of convergence in optimization algorithms, for example. Similarly, it was demonstrated in Bonnaire et al. [10] that the edges of bulk regions of the spectrum play an important role in the training dynamics of diffusion models. In doing so, we • develop the polynomial method to introduce an algebraic FD framework (FDPM) that generalizes previous methods beyond quadratic Stieltjes transforms, extending previous work to accurately approximate multi-scale spectral distributions, which may be multi-modal and/or contain multiple bulk regions; • derive equations to accurately compute the dynamics of important spectral quantities (atoms, moments, 2

and edges), as well as track the mechanism of splitting and merging of bulks, without the need to compute the full spectral distribution; • demonstrate empirically that the method works on several classes of ML models; and • deploy an implementation into the open source software project freealg. The remainder of the paper is structured as follows. Section 2 outlines some background and related work, as well as tools from free probability. Section 3 describes free decompression via the polynomial method (FDPM), detailing FD in this context, and the numerical methods we use, including associated methods for atoms, spectral edges, and moments. Section 4 provides numerical examples demonstrating the utility of our method. We conclude and discuss limitations in Section 5. Background, proofs, additional information and experiments are provided in a detailed set of appendices.

Background on Free Decompression

Free decompression was recently introduced by Ameli et al. [1] as a method for gaining access to the spectral distribution of so-called impalpable matrices: cases where the user does not have complete access to the matrix, perhaps due to their sheer size, corruption, or data availability. The method relies upon tools from free probability to derive a PDE in the familiar Stieltjes transform, a standard tool in RMT that at this stage has been widely deployed in both ML theory and the development of practical algorithms. We briefly recall some relevant notation. Let A denote a matrix of size n × n, where n is assumed to be large, and let mA denote the corresponding empirical Stieltjes transform, given by mA (z) :=

n 1X 1 1 tr(A − zI)−1 = , n n i=1 λi − z

(1)

where z ∈ C \ {λi }ni=1 and λi is the i-th eigenvalue of A. It is generally assumed that A is an element of a n×n sequence (An )∞ , and as n → ∞, mAn (z) → m(z) n=1 where each Ai is contained in Aj for j > i, An ∈ R pointwise. Here, m isR the Stieltjes transform of a limiting spectral measure ν in the sense that m : C\supp(ν) → 1 C satisfies m(z) := R x−z ν(dx). The R-transform R(z) corresponding to m is the (possibly multi-valued) 1 function satisfying R(−m(z)) = z + m(z) for z ∈ C \ supp(ν). To model the process of matrix subsampling, let σn denote a uniformly random permutation of {1, . . . , n} k and for 1 ≤ k ≤ n, let A(k) n = (Aσn (i)σn (j) )i,j=1 . Then for any subsampling factor τ < 1, EmA([τ n]) (z) → mτ (z) n as n → ∞. Ameli et al. [1] asserts that the map m 7→ mτ is equivalent to the map that operates on the R-transform by R(z) 7→ R(τ z) via the Nica-Speicher Theorem [30].1 When τ < 1, as is typical in free probability theory, this is called free compression. However, this operation can be easily inverted by taking τ > 1, moving from the Stieltjes transform of the submatrix to that of the full matrix; this inverse operation is free decompression. While a number of obvious strategies exist to implement free decompression, they each have significant drawbacks. The multi-valued nature of the R-transform, and the need to eventually recover m, make it difficult to work with directly. The approach taken by Ameli et al. [1] relies on numerically solving a PDE in m, where the major challenges faced are accurate curve-fitting of the approximated spectral distribution of the initial matrix, and analytic continuation of the Stieltjes transform across the x-axis. Reparameterization of this PDE results in an inviscid complex Burgers equation, and while this presents another potential approach, it does not seem straight-forward to implement in practice. While the formation of singularities and shocks can be expected for equations of this type, we have identified a family of distributions for which the flow exists for all time in Appendix F.2. For further discussion on these difficulties, we refer the reader to Appendix F.3. 1

Note that this theorem does not hold at this level of generality. A rigorous evaluation suggests it may need unitary invariance (see Appendix F.1 for a proof of this in our notation). Fortunately, the assumption does appear to hold for many cases of interest in modern ML, including the ones considered here.

Implementing the Polynomial Method

An especially wide array of operations in random matrix theory, including free additive convolution, free multiplicative convolution, can be distilled into operations acting on various transforms of the Stieltjes transform, including the R-transform, the S-transform, and the moment transform. We have seen that free decompression is defined in terms of the R-transform. These transforms are all related to each other by algebraic relations. This is inconvenient for performing direct numerical conversions between the transforms themselves [12, 33]. However, it is very convenient if the Stieltjes transform is already defined through an implicit algebraic relation, as deriving relations for each of the transforms and their inverses offers no further complexity. This observation was first made by Rao & Edelman [39], who coined the polynomial method as a numerical framework that moves between these transforms for performing matrix computations. At the heart of this approach is the following ansatz. Assumption 1. The Stieltjes transform m is algebraic: there exists a nonzero bivariate polynomial P : C × C → C such that for any z ∈ C \ K where K is a finite set of singularities, P (z, m(z)) = 0.

(2)

We write dz for the degree of this polynomial in z, and s for its degree in m. In addition to the quadratic examples previously handled by Ameli et al. [1], this unlocks, for example, models of Hessian matrices and neural network gram matrices [36, 38], which we demonstrate in Section 4. Under Assumption 1, Cauchy, R-, S- and moment transforms, as well as common algebraic operations and stochastic transformations, can also be realized as the roots of corresponding polynomials. Despite its convenience, Rao & Edelman [39] highlight several fundamental obstacles to the polynomial method preventing its implementation. While m is easily estimated, finding P is nontrivial. Polynomial roots are sensitive to their coefficients, and because valid P for algebraic Stieltjes transforms are rare among bivariate polynomials, naïve polynomial fits often radically fail. Furthermore, identifying the correct numerical solution for m from P alone is nearly impossible. Consequently, to our knowledge, no general implementations of the polynomial method exist. In the following sections, we describe our approach to rectify these challenges and provide its first implementation. Fitting the Bivariate Polynomial. To find solutions to (1), we need to fix the degree of the polynomial in z and m, respectively denoted by dz and s. These degrees can be identified post-hoc using standard model criteria [21, Section 6.1.3]. It turns out that by appealing to further properties of the Stieltjes transform, we can demonstrate that its corresponding polynomial which solves (1) must have real-valued coefficients (Lemma D.1). Consequently, P can be written as P (z, m) :=

dz s X X

cij z i mj ,

cij ∈ R.

(3)

i=0 j=0

In order to extract the coefficients from (3) however, we require access to samples of a faithful approximation m̂ to the Stieltjes transform itself. One option would be to use the curve-fitting procedures of Ameli et al. [1] or Olver & Nadakuditi [33], treating each sub-interval of I separately, then combining them into the Stieltjes transform approximation via quadrature. However, we found that for reasonably sized matrices, we achieve better results by direct evaluation of the empirical Stieltjes transform (1). To obtain an accurate approximation for the coefficients cij , samples locations for m̂ need to be carefully chosen. After significant trial-and-error, we found success by sampling z from logarithmic-scaling high-radius Bernstein ellipses around individual bulks of the eigenspectrum. Further improvements result from incorporating constraints using empirical moments. Our approach is discussed in Appendix D, with the complete algorithm presented in Algorithm D.1. The Physical Branch. While algebraic manipulations of P are straightforward, a central numerical challenge is identifying, for each z ∈ C, which root of P (z, ·) is the physical Stieltjes value. The physical branch is characterized globally by the Herglotz property ℑm(z) > 0 for z ∈ C+ and by the normalization 4

m(z) = −z −1 + O(z −2 ) at infinity. However, these conditions do not by themselves provide a stable pointwise root-selection rule. In high-degree or multi-bulk examples, roots can come close, exchange ordering, or pass near branch points, and relying only on pointwise tests can easily lead to jumps between sheets. Branch Selection via Continuation. We therefore select the physical branch by geometric continuation on the algebraic spectral curve P (z, m) = 0. Starting from an anchor point where the physical root can be identified reliably, either from the normalization at infinity or by comparison with the empirical Stieltjes transform, we continue the solution along a path in the z-plane. Along a smooth portion of the curve, implicit differentiation gives dm ∂ P (z, m) =− z . dz ∂m P (z, m) This ordinary differential equation respects the local geometry of the sheet and provides a predictor for transporting the physical root from the anchor to the query point. The predictor is then corrected by Newton projection onto P (z, m) = 0, with step subdivision when the geometry becomes stiff. This continuation procedure is substantially more stable than independent pointwise root selection and allows the fitted polynomial to be used as a numerical Stieltjes transform even for the high-degree, multi-scale curves used in our experiments; details are given in Appendix E.1.

3.1

Free Decompression on Algebraic Spectral Curves

Since the R-transform that underpins free decompression satisfies an algebraic relation whenever Assumption 1 holds, the Stieltjes transform after free decompression also satisfies an algebraic relation. This is stated in our main Theorem, proved as Theorem F.4 in Appendix F.4. Theorem 1. Let m be a Stieltjes transform satisfying Assumption 1, τ ≥ 1 be a decompression ratio, and mτ be the free decompression of m by the factor τ . Then every point (z, mτ (z)) ∈ C2 satisfies τ −1 P z+ , τ mτ = 0. (4) τ mτ Especially noteworthy here is the lack of any further assumptions on the nature of m or on the explicit construction of the R-transform. Assumption 1 naturally allows free decompression to be viewed through a different lens, characterizing it as a family of rational reparameterizations of the spectral curve P (z, m) = 0. This characterization is convenient, as it circumvents the manual analytic continuation that plagued the PDE-based approach in Ameli et al. [1]. Indeed, for each fixed z away from the singular set, the Fundamental Theorem of Algebra guarantees that P (z, m) = 0 has s solutions in m, up to multiplicity. Viewing m(z) as a multivalued algebraic function, these roots represent the sheets of the analytic continuation of the Stieltjes transform. The physical Stieltjes transform is one such sheet, and free decompression transports it along the spectral curve. This geometric viewpoint is developed in Appendix F.6; Figure 2 illustrates the sheet structure and the gluing of branches in the Marchenko–Pastur case. Evolving the Physical Branch. The evolved algebraic relation in Theorem 1, and its polynomial form in Corollary F.5, gives all candidate roots of mτ (z), but it does not identify the physical branch. To select the physical root, we track the known physical sheet at τ = 1 along the free-decompression flow. For a fixed query point z, we work in the original spectral-curve coordinates (ζ, y), where P (ζ, y) = 0, and impose the characteristic relation ζ − (τ − 1)y −1 − z = 0. Equivalently, (ζ, y) is determined by the coupled system F1 (ζ, y) := P (ζ, y) = 0,

F2 (ζ, y; τ, z) := ζ − (τ − 1)y −1 − z = 0.

At τ = 1, this system reduces to ζ = z and y = m0 (z), where m0 (z) is already known from the physical-branch selection procedure above. Differentiating the system with respect to τ , while keeping z fixed, gives the tangent equation ∂ζ P (ζ, y) ∂y P (ζ, y) 0 ζ̇ = . 1 (τ − 1)y −2 ẏ y −1 5

(a) m1 on C

(b) m2 on C

(d) m1 on C− glued to m2 on C+

1.0

Im(z)

0.5

0.0

−0.5 −1.0 1.0

Im(z)

0.5

0.0

−0.5 −1.0 −0.5

0.0

0.5

1.0

Re(z)

1.5

2.0

2.5 −0.5

0.0

0.5

1.0

1.5

2.0

2.5

Re(z)

Figure 2: Visualization of multiple solutions to the algebraic relation (1) for the Marchenko–Pastur density 1 with λ = 50 (left), and the surface obtained by gluing the Stieltjes solution m1 to the secondary branch m2 (right); colors are interpreted as in Appendix A. Free decompression requires tracing families of continuous paths along this surface, from the top sheet down to the lower sheet. This ordinary differential equation on the spectral curve provides the predictor direction for transporting the physical sheet from one decompression ratio to the next. The predicted point is then corrected by Newton projection onto the nonlinear system (F1 , F2 ) = (0, 0), with step subdivision when the local geometry becomes stiff. After correction, the decompressed Stieltjes value is recovered as mτ (z) = τ −1 y. The complete geometric predictor–corrector procedure is given in Appendix F.6.

3.2

Evolution of Spectral Quantities

One of the advantages of viewing free decompression through the lens of the polynomial method is that the evolution of a number of quantities of interest becomes clear. In particular, the weight of atomic components of distributions, as well as their moments, are easily computable. Similarly, the edges of the distribution’s absolutely-continuous density components evolve in a deterministic way. Atoms. An atom of a distribution at a point x◦ corresponds to a simple pole of the Stieltjes transform m(z) at z = x◦ , whose weight is given in terms of its residue w = limz→x◦ (x◦ − z)m(z). Under Assumption 1, Ps writing P (z, m) = i=0 ai (z)mi , a necessary condition for the existence of a pole at finite real x◦ is that x◦ is a root of the leading coefficient, i.e., as (x◦ ) = 0. Computing the residue gives w(x◦ ) = as−1 (x◦ )/a′s (x◦ ). Under free decompression, we can track this weight as a function of the decompression ratio τ . This is depicted in Figure 1, which shows the proportion of probability mass occupied by the atomic component as the matrix decompresses. The precise dynamics are encoded in the following proposition, whose proof is given in Appendix G.4. Proposition 2. Let x◦ ∈ R be an atom of m(z) with initial mass w◦ ∈ [0, 1]. Under the free decompression flow, the atom location remains fixed at x◦ (τ ) = x◦ , while its weight evolves via w(τ ) = 1 − τ −1 (1 − w◦ ). Spectral Edges. Under FD, it is also possible to track the support of the spectral density, whose bulk regions may split or merge as the spectral density evolves. This turns out to be a simpler task than accurately estimating the entire spectral density, and can often be used to diagnose and correct poor density approximations. We assume that the spectral distribution can be decomposed into a finiteSset of atoms and k absolutely-continuous bulk components. The support of the latter can be written as I(τ ) = j=1 Ij (τ ), where each Ij (τ ) = [aj (τ ), bj (τ )]. The spectral edges are then these boundary points {aj (τ ), bj (τ )}kj=1 . Under Assumption 1, these points have a precise characterization in terms of the branch points of the spectral curve.

(a) Initial Density (n = 4K)

(b) Free Decompression

ρ(λ)

n = 22.0 K n = 22.5 K n = 23.0 K n = 23.5 K n = 24.0 K n = 24.5 K n = 25.0 K

Fitted Empirical Histogram

0.05

Prediction Empirical Histogram

0.04 0.03 0.02

0.01

0 0

0.00 0

Figure 3: Example of free decompression for recovering the ESD of a 32,000 × 32,000 matrix, generated from the compound free Poisson law described in Section 4. (a) ESD of the sampled 4000 × 4000 submatrix (histogram) and the fitted density (solid curve). (b) Free decompression flow over n. (c) Final predicted ESD at size n = 32,000, compared with the ESD of the full matrix. These are points at which different sheets of the spectral curve coalesce, satisfying P (ζ, y) = 0,

and

∂P ∂P (ζ, y) − (τ − 1) (ζ, y) = 0, ∂y ∂ζ

ζ ∈ R,

where y belongs to the physical branch, as derived in Proposition G.1 in Appendix G.1. To identify which candidate points are spectral edges requires branch selection via the continuity method from known points on the physical branch, similar to the method described in Appendix F.6. At τ = 1, we use nearby points to identify which roots come from the physical branch, and for τ > 1 we use points identified in previous time-steps. An algorithmic description of a procedure for finding the spectral edges is given in Algorithm G.1. For some positive-definite matrices, FD can predict that the left edge of the support converges to zero much faster than it does in practice. We can remedy this through a finite-size correction; see Appendix H for details. Merging and Splitting Bulks. As distributions evolve under free decompression, we may encounter situations where either multiple bulks merge into one continuous region, or a connected bulk region splits into separate bulks. Fortunately, we can identify points where this occurs. Structurally, these points are second-order branch points of the spectral curve, identified as satisfying y∗2 ∂y P − (τ∗ − 1) ∂ζ P = 0, and (∂ζζ P )(∂y P )2 − 2(∂ζy P )(∂ζ P )(∂y P ) + (∂yy P )(∂ζ P )2 + 2(∂ζ P )2 (∂y P ) = 0, P (ζ∗ , y∗ ) = 0,

y∗

as shown in Proposition G.2 in Appendix G.2. We derive similar equations for the dynamics of the gaps between bulk regions, as well as the support intervals in Appendix G.3. R Moments. Computing the moments µk = xk ν(dx) of a spectral distribution ν of interest is also of significant importance. Since moments are easy to estimate from eigenvalues, they comprise valuable criteria to ensure an effective polynomial fit to the Stieltjes transform P (see Appendix D.1). For example, a necessary condition for P to have a Stieltjes transform as a root is for i−j=emax cij (−1)j = 0 where emax := max{i−j : cij ̸= 0} (see Lemma C.1 in Appendix C.1). Proposition C.2 provides an efficient O(n2 ) recurrence relation to compute the first n moments µ1 , . . . , µn directly from P , without root selection. Furthermore, the moments can be evolved under FD without performing FD on the density itself. Specifically, the n-th moment of a Pn−1 ) k probability measure ν under FD with ratio τ satisfies µ(τ n = k=0 κn,k τ where the coefficients κn,k are computed in Proposition G.7 and depend only on the moments µ1 , . . . , µn of ν.

(a) Initial Density (n = 4K)

2.5

ρ(λ)

2.0

(b) Free Decompression

2.5

Fitted Empirical Histogram

2.0

1.5

1.0

0.5

0.0 0

Prediction Empirical Histogram

2.0 1.5 1.0 0.5

0.0 −1

2.5

n = 22.0 K n = 22.5 K n = 23.0 K n = 23.5 K n = 24.0 K n = 24.5 K n = 25.0 K

0.0 −1

−1

Figure 4: Example of free decompression for recovering the ESD of a 32,768 × 32,768 matrix, generated from the neural network Hessian described in Section 4. (a) ESD of the sampled 4000 × 4000 submatrix (histogram) and the fitted density (solid curve). (b) Free decompression flow from n0 to n. (c) Final predicted ESD at size n = 32,768, compared with the ESD of the original matrix.

Empirical Results

We now consider three examples that demonstrate the ability to decompress distributions of interest in ML. First, we consider an approximation of a Compound Free Poisson Law that contains both atoms and splitting bulk regions. Then, we consider the Pennington–Bahri model of the limiting spectral distributions of Hessian matrices of neural networks at critical points of the MSE loss [36]. Finally, we consider the activation matrix of a random Fourier features-based diffusion model. Compound Free Poisson. Our first example belongs to the family of compound free Poisson distributions (see Appendix I.1 for details). In our example we take ν(dx) = 34 δ2 (dx) + 14 δ5.5 (dx), and the distribution can be viewed as ( 14 δ2 (dx) ⊠ µM P (λ)) ⊞ ( 34 δ2 (dx) ⊠ µM P (λ)). The target density has aspect ratio λ = 0.1, and consists of a spike at 0, and two absolutely continuous bulk regions. A 6000 × 6000 matrix sampled with these properties, and a 1000 × 1000 submatrix was subsampled and used as initial data, and we remark that its corresponding law (and the spectral density approximated from its polynomial fit) consists of a single bulk region. Free decompression on this submatrix is depicted in Figure 3. This example clearly demonstrates both the atom dynamics and bulk-splitting behavior described in Section 3.2, as shown in Figure 1. Neural Network Hessian. Our second example is the Hessian matrix of the ℓ2 -loss of a two-layer neuralnetwork autoencoder model with CIFAR-10 data, at initialization. Pennington & Bahri [36] gave a free probabilistic approximation of the Hessian, which is treated separately in Appendix I.3. Noting that the Hessian matrix decomposes as H = H0 + H1 , with H0 being the squared-Jacobian component, and H1 the network Hessian term, the authors treated each of these components separately in a similar but smaller autoencoder model to the one we consider here. Full details on the architecture, training, and matrix formation can be found in Appendix I.4. For the purposes of decompression, the full 32,768×32,768 matrix was formed, and 4000×4000 submatrices were sampled and their eigenvalues combined to form an empirical spectral density. After fitting the bivariate polynomial, the Stieltjes transform was expanded under free decompression back to the full matrix, as shown in Figure 4. In Appendix I.4, we use the decompressed distributions to accurately track the index of the Hessian subsamples. We remark that while the full Hessian was computed as a baseline for our experiments, subsamples of Hessian matrices too large to practically compute can be readily obtained by masking using standard automatic differentiation tools. Diffusion Model. Our third example considers the expected activation matrix of a diffusion model studied in George et al. [16]; see Appendix I.5 for model and matrix details. Its large-scale limit, studied in Bonnaire et al. [10], has a highly multi-scale spectrum with several bulks, including a narrow atom-like bulk away from zero, making this a stringent log-scale test case for branch selection and density recovery. 8

(a) Initial Density (n = 4K)

Density evolution from n = 4K to n = 64K

Fitted Empirical Histogram

102

ρ(λ)

n = 22.0 K n = 22.5 K

10−1 10−4

n = 23.0 K

10−7

n = 23.5 K

(b) Final Density (n = 64K)

n = 24.0 K

ρ(λ)

n = 24.5 K n = 25.0 K n=2

5.5

n=2

6.0

10−5

Prediction Empirical Histogram

10−1 10−4

10−4

10−3

10−2

10−1 λ

100

101

102

103

10−7 10−5

10−4

10−3

10−2

10−1

100

101

102

103

Figure 5: Free decompression for recovering the ESD of an n × n activation matrix A of a random-features diffusion model with n = 64,000 random feature weights. Left: ridge plot of the FD density evolution on log scale (solid lines), compared with empirical spectral distributions of sampled submatrices (colored histograms). Right: (a) ESD of the sampled n0 × n0 matrix A0 , with n0 = 4000, together with the fitted density; and (b) final predicted ESD at n = 64,000, compared with the ESD of the original matrix A. The y-axis is logarithmic in all panels; the leftmost narrow bulk is atom-like and stands ten billion times higher than the right bulk.

n(τ )

ρ(λ, τ )

Evolution of Spectral Edges To perform free decompression, we 26 K 103 formed a 64,000 × 64,000 matrix and sampled a 4000 × 4000 submatrix. The eigenval101 25 K ues of the subsampled matrices were pooled to obtain a more accurate initial ESD. Fig10−1 ure 5 (left) shows ridge plots of the FD evo4 2 K 10−3 lution up to 64K, overlaid with empirical spectra from subsamples at each scale. The 23 K 10−5 right panels isolate the initial and final densities on the same vertical scale, showing 10−7 that the final FD prediction matches the 22 K −5 −4 −3 −2 −1 0 1 2 3 10 10 10 10 10 10 10 10 10 empirical spectrum across the full log-scale λ range. The bulk-edge behavior is also impor- Figure 6: Evolution of the ESD of the diffusion model. Solid tant: the leftmost spike-like bulk follows black curves show the predicted spectral edges, including a cusp atom-like mass dynamics, while its smallest point. The edge behavior closely matches predictions, except empirical eigenvalue is affected by finite-size for the leftmost edge. Our finite-size correction of this edge is extreme-value effects. Figure 6 shows the shown in the dashed black line. The blue region converges to edge evolution and bulk splitting. At these an atom. logarithmic scales, the leftmost empirical edge does not follow the macroscopic FD flow, but is captured by the finite-size correction of Appendix H. Overall, this example demonstrates that FDPM can evolve multi-scale, multi-bulk spectra even when branch selection, log-scale density recovery, and finite-size edge effects are all simultaneously active.

Conclusions

We have developed a new method of free decompression that generalizes previous work by representing the Stieltjes transform of the spectral density through a bivariate polynomial of arbitrary degree. Under this algebraic ansatz, free decompression becomes an evolution on spectral curves, allowing us to track the physical branch, evolve the density, and derive equations for spectral quantities such as atoms, moments, 9

Table 1: Comparison of direct eigendecomposition with free decompression (FD) for the diffusion-model random-feature matrix. The Direct column reports the process time for computing the eigenvalues at each matrix size. In the FD column, the first term is the time to compute eigenvalues of the initial 22 K-sized submatrix, and the second term is the FD evolution time. Distributional distances compare the empirical spectrum with the FD reconstruction, with W1 /L denoting the Wasserstein-1 distance on log-eigenvalues normalized by the log-width L of the plotted spectral range, and MMD denoting the Gaussian-kernel maximum mean discrepancy on log-eigenvalues. Moment errors are relative errors of the first moment and standard deviation. Log edge errors compare empirically detected spectral edges with FD-evolved edges on the log-λ axis, excluding the unstable left edge of the spike-like leftmost bulk. All reported errors are multiplied by 100 and shown as percentages. Size

Process Time (sec)

Distributional Distance

Moments Rel. Error

Log Edges Error

Direct

FD (ours)

W1 /L

MMD

∆µ1 /µ1

∆σ/σ

Mean

Max

22 K 23 K 24 K 25 K 26 K

32.2 431.43 3534.5 27056.2 192447.7

32.2 + 0.0 32.2 + 2.3 32.2 + 2.9 32.2 + 4.3 32.2 + 6.8

0.01% 0.22% 0.43% 0.70% 2.19%

0.01% 0.08% 0.12% 0.15% 0.53%

0.01% 0.02% 0.01% 0.05% 0.08%

0.01% 0.03% 0.00% 0.08% 0.09%

0.1% 0.2% 0.2% 0.4% 0.4%

0.2% 0.5% 0.2% 0.8% 1.2%

and bulk edges. This enables FD to handle substantially more challenging spectra than previous approaches, including multi-modal, multi-bulk, and multi-scale distributions arising in realistic ML matrices. Limitations. While the current method expands the range of matrices accessible to free decompression, some important regimes remain challenging. In particular, heavy-tailed spectral distributions are still elusive with the present tools. In addition, although the finite-size correction introduced here improves the behavior of near-singular edges at log scale, further refinement may be needed for highly accurate edge estimation in extreme cases. We hope these results provide a foundation for further development of algebraic spectral-curve methods for large-scale spectral inference.

References [1] Ameli, S., van der Heide, C., Hodgkinson, L., & Mahoney, M. W. (2025a). Spectral estimation with free decompression. In The Thirty-ninth Annual Conference on Neural Information Processing Systems. [2] Ameli, S., van der Heide, C., Hodgkinson, L., Roosta, F., & Mahoney, M. W. (2025b). Determinant estimation under memory constraints and neural scaling laws. In Forty-second International Conference on Machine Learning. [3] Anshelevich, M. (2008). Orthogonal polynomials with a resolvent-type generating function. Transactions of the American Mathematical Society, 360(8), 4125–4143. [4] Benaych-Georges, F. (2005). Classical and free infinitely divisible distributions and random matrices. The Annals of Probability, 33(3), 1134 – 1170. [5] Bercovici, H. & Voiculescu, D. (1993). Free convolution of measures with unbounded support. Indiana University Mathematics Journal, 42(3), 733–773. [6] Biroli, G., Bonnaire, T., de Bortoli, V., & Mézard, M. (2024). Dynamical regimes of diffusion models. Nature Communications, 15(1), 9957. [7] Blaizot, J.-P., Grela, J., Nowak, M., & Warchoł, P. (2015). Diffusion in the space of complex Hermitian matrices: microscopic properties of the averaged characteristic polynomial and the averaged inverse characteristic polynomial. Acta Physica Polonica. B, 46(9), 1801–1823. [8] Blaizot, J.-P. & Nowak, M. A. (2010). Universal shocks in random matrix theory. Phys. Rev. E, 82, 051115. [9] Blaizot, J.-P., Nowak, M. A., & Warchoł, P. (2013). Universal shocks in the Wishart random-matrix ensemble. Phys. Rev. E, 87, 052134. [10] Bonnaire, T., Urfin, R., Biroli, G., & Mezard, M. (2025). Why diffusion models don’t memorize: The role of implicit dynamical regularization in training. In The Thirty-ninth Annual Conference on Neural Information Processing Systems. [11] Burda, Z., Janik, R. A., Jurkiewicz, J., Nowak, M. A., Papp, G., & Zahed, I. (2002). Free random Lévy matrices. Phys. Rev. E, 65, 021106. [12] Chen, J. & Olver, S. (2026). Mathematics of Computation.

Computing inverses of Stieltjes transforms of probability measures.

[13] Couillet, R. & Liao, Z. (2022). Random Matrix Methods for Machine Learning. Cambridge University Press. [14] El Karoui, N. (2010). The spectrum of kernel random matrices. The Annals of Statistics, 38(1), 1 – 50. [15] Erdős, L., Krüger, T., & Schröder, D. (2020). Cusp universality for random matrices I: Local law and the complex Hermitian case. Communications in Mathematical Physics, 378(2), 1203–1278. [16] George, A. J., Veiga, R., & Macris, N. (2025). Denoising score matching with random features: Insights on diffusion models from precise learning curves. [17] Gibbs, A. L. & Su, F. E. (2002). statistical review, 70(3), 419–435.

On choosing and bounding probability metrics.

International

[18] Guionnet, A. (2009). Large random matrices: lectures on macroscopic asymptotics. [19] He, K., Zhang, X., Ren, S., & Sun, J. (2016). Deep residual learning for image recognition. Proceedings of the IEEE conference on computer vision and pattern recognition (pp. 770–778).

[20] Hodgkinson, L., Wang, Z., & Mahoney, M. W. (2025). Models of heavy-tailed mechanistic universality. In A. Singh, M. Fazel, D. Hsu, S. Lacoste-Julien, F. Berkenkamp, T. Maharaj, K. Wagstaff, & J. Zhu (Eds.), Proceedings of the 42nd International Conference on Machine Learning, volume 267 of Proceedings of Machine Learning Research (pp. 23290–23329).: PMLR. [21] James, G., Witten, D., Hastie, T., & Tibshirani, R. (2013). An introduction to statistical learning: with applications in R, volume 103. Springer. [22] Kesten, H. (1959). Symmetric random walks on groups. Transactions of the American Mathematical Society, 92(2), 336–354. [23] Krizhevsky, A. (2009). Learning Multiple Layers of Features from Tiny Images. Technical report, University of Toronto. [24] Liao, Z., Couillet, R., & Mahoney, M. W. (2021). A random matrix analysis of random Fourier features: beyond the Gaussian kernel, a precise phase transition, and the corresponding double descent. Journal of Statistical Mechanics: Theory and Experiment, 2021(12), 124006. [25] Martin, C. H. & Mahoney, M. W. (2021). Implicit self-regularization in deep neural networks: Evidence from random matrix theory and implications for learning. Journal of Machine Learning Research, 22(165), 1–73. [26] Marčenko, V. A. & Pastur, L. A. (1967). Distribution of eigenvalues for some sets of random matrices. Mathematics of the USSR-Sbornik, 1(4), 457. [27] McKay, B. D. (1981). The expected eigenvalue distribution of a large regular graph. Linear Algebra and its Applications, 40, 203–216. [28] Meixner, J. (1934). Orthogonale polynomsysteme mit einer besonderen gestalt der erzeugenden funktion. Journal of The London Mathematical Society-second Series, (pp. 6–13). [29] Muskhelishvili, N. (1992). Singular Integral Equations: Boundary Problems of Function Theory and Their Application to Mathematical Physics. Dover Books on Mathematics Series. Dover Publications. [30] Nica, A. & Speicher, R. (1996). On the multiplication of free N-tuples of noncommutative random variables. Amer. J. Math., 118(4), 799––837. [31] Nica, A. & Speicher, R. (2006). Lectures on the Combinatorics of Free Probability. London Mathematical Society Lecture Note Series. Cambridge University Press. [32] Olver, F. W. J., Olde Daalhuis, A. B., Lozier, D. W., Schneider, B. I., Boisvert, R. F., Clark, C. W., Miller, B. R., Saunders, B. V., Cohl, H. S., & McClain, M. A. (Accessed 3 May 2026). NIST Digital Library of Mathematical Functions. https://dlmf.nist.gov/. [33] Olver, S. & Nadakuditi, R. R. (2012). Numerical computation of convolutions in free probability theory. arXiv preprint arXiv:1203.1958. [34] Panaretos, V. M. & Zemel, Y. (2019). Statistical aspects of Wasserstein distances. Annual Review of Statistics and its Application, 6(1), 405–431. [35] Pastur, L. & Shcherbina, M. (2011). Eigenvalue Distribution of Large Random Matrices. Mathematical surveys and monographs. American Mathematical Society. [36] Pennington, J. & Bahri, Y. (2017). Geometry of neural network loss surfaces via random matrix theory. In D. Precup & Y. W. Teh (Eds.), Proceedings of the 34th International Conference on Machine Learning, volume 70 of Proceedings of Machine Learning Research (pp. 2798–2806).: PMLR. [37] Pennington, J., Schoenholz, S., & Ganguli, S. (2018). The emergence of spectral universality in deep networks. In A. Storkey & F. Perez-Cruz (Eds.), Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics, volume 84 of Proceedings of Machine Learning Research (pp. 1924–1932).: PMLR. 12

[38] Pennington, J. & Worah, P. (2017). Nonlinear random matrix theory for deep learning. In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, & R. Garnett (Eds.), Advances in Neural Information Processing Systems, volume 30: Curran Associates, Inc. [39] Rao, N. R. & Edelman, A. (2008). The polynomial method for random matrices. Foundations of Computational Mathematics, 8(6), 649–702. [40] Rudelson, M. & Vershynin, R. (2009). Smallest singular value of a random rectangular matrix. Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences, 62(12), 1707–1739. [41] Saitoh, N. & Yoshida, H. (2001). The infinite divisibility and orthogonal polynomials with a constant recursion formula in free probability theory. Probab. Math. Statist., 21, 159–170. [42] Senouf, D., Caflisch, R., & Ercolani, N. (1996). Pole dynamics and oscillations for the complex Burgers equation in the small-dispersion limit. Nonlinearity, 9(6), 1671. [43] Shlyakhtenko, D. & Tao, T. (2022). Fractional free convolution powers. Indiana University Mathematics Journal, 71(6). [44] Speicher, R. (2015). Free probability theory. In The Oxford Handbook of Random Matrix Theory. Oxford University Press. [45] Stein, E. M. & Shakarchi, R. (2003). Fourier Analysis: An Introduction, volume 1 of Princeton Lectures in Analysis. Princeton, NJ: Princeton University Press. [46] Sturmfels, B. (2002). Solving Systems of Polynomial Equations. Conference Board of the Mathematical Sciences Regional Confe. Conference Board of the Mathematical Sciences. [47] Trefethen, L. N. (2019). Approximation Theory and Approximation Practice. Philadelphia, PA: Society for Industrial and Applied Mathematics, extended edition edition. [48] Trefethen, L. N. (2023). Numerical analytic continuation. Japan Journal of Industrial and Applied Mathematics, 40(3), 1587–1636. [49] Ventura, E., Achilli, B., Silvestri, G., Lucibello, C., & Ambrogioni, L. (2025). Manifolds, random matrices and spectral gaps: The geometric phases of generative diffusion. In The Thirteenth International Conference on Learning Representations. [50] Voiculescu, D. (1991). Limit laws for random matrices and free products. Inventiones mathematicae, 104(1), 201–220. [51] Wachter, K. W. (1978). The strong limits of random matrix spectra for sample matrices of independent elements. The Annals of Probability, 6(1), 1–18. [52] Wegert, E. (2012). Visual Complex Functions: An Introduction with Phase Portraits. Mathematics and Statistics. Springer Basel. [53] Wei, A., Hu, W., & Steinhardt, J. (2022). More than a toy: Random matrix models predict how real-world neural representations generalize. In International Conference on Machine Learning (pp. 23549–23588).: PMLR. [54] Wigner, E. P. (1955). Characteristic vectors of bordered matrices with infinite dimensions. Annals of Mathematics, 62(3), 548–564.

Appendices Contents A Notation and Conventions

B Background 16 B.1 Stieltjes Transform and Boundary Values . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 B.2 Analytic Continuation Beyond the Physical Branch . . . . . . . . . . . . . . . . . . . . . . . . 17 C Algebraic Stieltjes Transforms 17 C.1 Moments from the Algebraic Relation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 D Fitting the Spectral Curve 20 D.1 Constraints at Infinity from Moments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 D.2 Constraints from Branch Points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 D.3 Numerical Reconstruction from Samples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 E Computing the Stieltjes Transform and Density 24 E.1 Physical Root Selection for the Stieltjes Transform . . . . . . . . . . . . . . . . . . . . . . . . 24 E.2 Density Recovery and δ-Extrapolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 F Free Decompression on Algebraic Spectral Curves 29 F.1 The Nica–Speicher Theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 F.2 Infinite Decompressibility . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 F.3 The PDE Approach and Challenges . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 F.4 Formal Free Decompression along Spectral Curves . . . . . . . . . . . . . . . . . . . . . . . . 34 F.5 Polynomial Evolution under Free Decompression . . . . . . . . . . . . . . . . . . . . . . . . . 35 F.6 Free Decompression on Spectral Curves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 G Evolution of Spectral Quantities 38 G.1 Evolving Spectral Edges . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 G.2 Spectral Cusps and The Mechanism of Merging Bulks . . . . . . . . . . . . . . . . . . . . . . 41 G.3 Dynamics of Gap and Support Intervals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 G.4 Dynamics of Atoms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 G.5 Evolution of Moments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 H Finite-Size Correction to the Spectral Curve I

Random Matrix Models 51 I.1 Compound Free Poisson Law . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 I.2 Free Lévy Law . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 I.3 Pennington–Bahri Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 I.4 Hessian of Two-Layer Network . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 I.5 Details of the Diffusion Model Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

J Implementation and Reproducibility Guide

Appendix A

Notation and Conventions

Nomenclature. We use boldface lowercase letters for vectors, boldface upper case letters for matrices, and normal face letters for scalars, including the components of vectors and matrices. Table A.1 summarizes the main symbols and notations used throughout the paper, organized by context.

Context

Table A.1: Common notations used throughout the manuscript. Symbol Description

Matrices & Models

A, An A(k) n I λi trace(·) trn (·) Qn H, J Σ W, X, Z U λ, c

General n × n matrix (often symmetric or Hermitian) k × k principal submatrix of An The identity matrix The i-th eigenvalue of a matrix Standard trace of a matrix Normalized trace, n−1 trace(·) Haar-distributed random orthogonal matrix Hessian matrix (and subcomponents) and Jacobian matrix Population covariance matrix Random matrices with i.i.d. entries (e.g., Wigner, Gaussian) Expected diffusion activation matrix (random Fourier features) Aspect ratios (p/n) or intensities/rates for random matrix models

Free Probability

ρ I m ω H[ρ] R S ν ν ⊞s Dc ν M (z) χ(z) rk (ν) µp , µp (t) ⊞ ⊠

Absolutely continuous empirical spectral density, R → R+ Compact support interval of density ρ, I = [λ− , λ+ ] ⊊ R R Stieltjes transform of density, m(z) = R ρ(y) y−z dy, C \ I → C Functional inverse of the Stieltjes transform, ω(m(z)) = z R ρ(y) Hilbert transform of density, H[ρ](x) = p. v. x−y dy, R → R 1 =z Voiculescu’s R-transform, R(z) = ω(−z) − z with R ω(m(z)) −1 zy S-transform, S(z) = 1+z ψ (z), where ψ(z) = ρ(y) dy z R 1−zy Probability measure associated with spectral density ρ s-fold free additive convolution power of ν Dilation of measure ν by c, i.e., push-forward under x 7→ cx Truncated moment-generating function Functional inverse of the truncated moment-generating function k-th free cumulant of the measure ν p-th moment of a probability measure (and its time evolution) Free additive convolution Free multiplicative convolution

Algebraic Curve

P (z, m) C

Algebraic equation for the Stieltjes transform, P (z, m(z)) = 0 Complex 1D spectral curve C = {(z, m) ∈ C2 | P (z, m) = 0}

Free Decompression

n ns τ, t ρ0 , m0

Size of the target (large) matrix Size of the sampled (small) sub-matrix, ns < n Decompression scale, τ = nn , t = log(τ ) s Density and Stieltjes transform at t = 0

Notation

C± ℜ, ℑ p. v.

Upper (lower) half complex plane Real and imaginary part of complex variable Cauchy principal value

Color code for plots. In our complex-plane plots, colors encode the values of a complex quantity using domain coloring: hue represents phase, while brightness represents magnitude; see the plain version in Figure A.1(b). For high-dynamic-range fields, we also use the tiled version shown in Figure A.1(a), in which phase and logarithmic magnitude are periodically reset, following the phase-portrait convention of Wegert [52, Section 2.5].

−1

−2

−2 −2

Re(z)

−1

−3

∞

−π

−2

1 2

|m(z)|

Im(z)

−3

(b) Plain Domain Coloring

arg(m(z))

(a) Tiled Domain Coloring

Re(z)

Figure A.1: Domain coloring of a Stieltjes-transform branch for the Free Lévy example in Figure I.3. Right: plain domain coloring, with hue encoding arg m(z) and brightness encoding |m(z)|. Left: tiled domain coloring, where phase and logarithmic magnitude are periodically reset to make high-dynamic-range structure more visible.

Appendix B B.1

Background

Stieltjes Transform and Boundary Values

Let ν be a probability measure on R with compact support supp(ν) ⊂ R. When ν is absolutely continuous with respect to Lebesgue measure dx, we denote its density by ρ(x) := ν(dx)/dx as the Radon–Nikodym derivative of ν. The Stieltjes transform of ν is the function m : C \ supp(ν) → C defined by Z 1 m(z) := ν(dx), z ∈ C \ supp(ν). (B.1) R x−z The Stieltjes transform is a Herglotz function: ℑm(z) > 0 for z ∈ C+ . Here, C+ := {z ∈ C | ℑz > 0}, and ℑw denotes the imaginary part for w ∈ C. Moreover, m satisfies the normalization at infinity: m(z) = −z −1 + O(z −2 ) as |z| → ∞. The function m admits non-tangential boundary values on the real line, denoted by m± (x) := limy→0± m(x+ iy) when the limits exist. If ν(dx) = ρ(x) dx in a neighborhood of x, the Sokhotski–Plemelj formula [29, p. 42] yields m± (x) = −πH{ρ}(x) ± iπρ(x), R ρ(t) where H{ρ}(x) := π1 p.v. R x−t dt is the Hilbert transform of ρ and p.v. denotes the Cauchy’s principal value of a singular integral. The density can be recovered from boundary values by ρ(x) = π −1 ℑm+ (x) = −π −1 ℑm− (x),

(B.2)

which is the inversion of the Stieltjes transform. The jump of m across the real axis is m+ (x) − m− (x) = 2iπρ(x). In particular, for x ∈ / supp(ρ), we have ρ(x) = 0, hence the boundary values coincide, and m is holomorphic across the real axis at x. In contrast, for x ∈ supp(ρ) the boundary values differ, which manifests as a 16

branch cut in the boundary values of m across supp(ρ). An example of a branch cut is shown in Figure 2(a), illustrated by the white line segment in the center of the plot. In the remainder, we denote the support of the absolutely continuous part of µ by I := supp(ρ), and we will typically assume that it consists of finitely many disjoint compact intervals, k [ I= Ij , Ij := [aj , bj ], j=1

where k ≥ 1, and aj < bj < aj+1 . We refer to k as the number of spectrum bulks, which is also the number of branch cuts of m. Also, any atomic component of µ produces isolated poles of m on the real axis rather than branch cuts, and will be treated separately when relevant.

B.2

Analytic Continuation Beyond the Physical Branch

Several applications, and in particular the method of characteristics used in free decompression [1], require evaluating the Stieltjes field m(z) along trajectories in the complex z-plane that cross the cuts I. Since the Stieltjes transform m(z) is defined as a single-valued holomorphic function only on C \ I, such evaluations necessarily invoke analytic continuation across I. It is therefore natural to regard m as belonging to a larger object: a multi-valued analytic function of z, whose single-valued realizations correspond to branches (or sheets) of a Riemann surface lying over the z-plane. We denote the Stieltjes transform defined by (B.1) by m1 and refer to it as the physical branch, characterized uniquely by the Herglotz property together with the normalization at infinity. Other branches, when they exist, are denoted by mi , i > 1, and are referred to as non-physical. Crossing a cut exchanges boundary traces of branches. Specifically, for x in the interior of an interval Ij , analytic continuation across Ij identifies the upper trace of one branch with the lower trace of another, and vice versa. Equivalently, two sheets are glued along Ij by matching boundary values with opposite orientations, so that a trajectory approaching Ij from C+ and crossing the cut may emerge on a different branch on the other side; see Figure 2 for an illustration. In our setting, the physical branch m1 is assumed known, either because µ is a model law or because ρ is estimated from an empirical spectrum and m1 is computed by (B.1). In contrast, the continuation of m1 beyond the physical domain is not determined uniquely by these data. In particular, the phrase “the other branch” is generally ambiguous: beyond matching boundary traces on I, there is no canonical choice of a global holomorphic continuation, and the ambiguity becomes more pronounced for multi-bulk spectra, where continuations around different endpoints {aj , bj }kj=1 can lead to distinct sheets. A convenient way to regularize this ill-posedness is to restrict attention to model classes in which the analytic continuation of m has only finitely many branches. This reduces continuation from an infinitedimensional ambiguity to a finite-sheeted branched covering of the z-plane. We adopt an algebraic model class, which we describe next.

Appendix C

Algebraic Stieltjes Transforms

Our structural assumption is that the Stieltjes transform m is algebraic: there exists a nonzero bivariate polynomial P such that P (z, m) = 0. (C.1) We write s := degm (P ) for the degree of P with respect to its second argument. For each z away from the branch locus, the algebraic equation P (z, m) = 0 admits s (not necessarily distinct) solutions in m, and the corresponding analytic branches are denoted by {mi (z)}si=1 . In this notation, m1 is the physical branch, singled out among the s branches by the Herglotz property together with the normalization at infinity, while the remaining branches are non-physical. Branch points occur at those z∗ for which two or more roots of P collide. Equivalently, there exists m∗ ∈ C such that P (z∗ , m∗ ) = 0, ∂m P (z∗ , m∗ ) = 0. (C.2) Eliminating m∗ yields the discriminant ∆(z) := Discm P (z, ·) ∈ R[z], whose zeros identify the branch locus of the multi-valued function defined implicitly by (C.1). In Stieltjes settings with conjugation symmetry, the 17

branch points z∗ on the real line include the spectral edges {aj , bj }kj=1 of the cut set I, though the algebraic model may admit additional branch points away from the real axis. A branch cut I is a choice of curves connecting branch points such that each branch mi becomes singlevalued and holomorphic on C \ I. Once I is fixed, each branch admits non-tangential boundary values on I, denoted by m± i . Across a cut component, analytic continuation can only permute branches. It is therefore convenient to collect branches into a vector, m(z) := [m1 (z), m2 (z), . . . , ms (z)]⊺ . For x ∈ Ij◦ , analytic continuation across Ij permutes the boundary traces, hence there exists a permutation matrix Πj ∈ {0, 1}s×s such that m+ (x) = Πj m− (x), x ∈ Ij◦ , where m± (x) denote non-tangential boundary values taken from C± . The collection {Πj }kj=1 encodes how the s sheets are glued along the real cuts and provides a finite description of the continuation rules needed to follow trajectories across I. In Appendix D, we describe the construction of a finite-sheet continuation from samples of the physical branch m1 by fitting (C.1). Before that, in Appendix C.1, we describe moment-based necessary conditions for the feasibility of such an algebraic representation.

C.1

Moments from the Algebraic Relation

One of the biggest challenges with the algebraic approach to performing free decompression is the identification of the root m(z) in the algebraic relation that corresponds to the principal branch m1 (z) of the Stieltjes transform we seek. This identification can be performed in one of two ways. The first is by homotopy, where one takes advantage of continuity and identifies the correct root at time t + δt using knowledge of the correct root at time t. This approach performs well when the roots are distinct, but can fail in regimes where the roots become close or “cross over”. The second approach involves the computation of moments; while this is less direct, it is arguably more concrete and provides a powerful check to ensure the correct solution is obtained. Let µp denote the p-th moment of the probability measure ν Z µp := xp ν(dx), p ∈ N0 , (C.3) R

where µ0 = 1 corresponds to unit mass normalization. If m1 denotes the Stieltjes transform of ν, then for sufficiently large |z|, m1 satisfies the Laurent expansion m1 (z) =

1 1 ν(dx) = − x − z z R

Z ∞ ∞ X µp 1 1X 1 p ν(dx) = − x ν(dx) = − p x p+1 , z p=0 z R R 1− z p=0 z

(C.4)

providing a concrete characterisation of the Stieltjes transform in terms of its moments. P Suppose that a polynomial parameterized by an index set A ⊂ N20 , P (z, m) = (i,j)∈A cij z i mj , is given. Our objective is to find a solution m(z) satisfying P (z, m(z)) = 0 and an expansion of the form m(z) = P∞ P∞ −k−1 for all sufficiently large |z|. Let w = 1/z and S(w) = k=0 ϑk wk so that m(z) = wS(w). If k=0 ϑk z m(z) corresponds to the principal branch m1 (z) with moments µk , then ϑ0 = −1 and µk = ϑk /ϑ0 = −ϑk for each k = 1, 2, . . . . P∞ −k−1 Lemma C.1. There exists a function m(z) satisfying P (z, m(z)) = 0 and m(z) = for k=0 ϑk z sufficiently large |z| if and only if ϑ0 is a root of the polynomial X L(ϑ) = cij ϑj , emax = max{i − j : (i, j) ∈ A, cij ̸= 0}. (i,j)∈A i−j=emax

Proof. Substituting m(z) = wS(w) into P (z, m(z)) and multiplying by wemax gives X cij wemax −(i−j) S(w)j = 0. (i,j)∈A

Taking w → 0+ gives L(ϑ0 ) = 0, as required. The existence of the root also guarantees the solution m(z) has the corresponding series expansion by the Newton–Puiseux Theorem. Example. As an example, consider the polynomial P (z, m) = λσ 2 zm2 + (z − σ 2 (1 − λ))m + 1 for the Marchenko–Pastur distribution with scale σ 2 and aspect ratio λ. In this case, emax = 0 and L(ϑ) = 1 + ϑ, which does indeed have a single root at ϑ = −1. In many practical cases, a root of P will only approximate a Stieltjes transform, and so it may only be the case that ϑ0 ≈ −1. Fortunately, ϑ0 can be estimated by finding the closest root of L(ϑ) to −1. As long as the remaining coefficients ϑk can be recovered, the “moments” of the corresponding Stieltjes transform can be estimated by µk ≈ ϑk /ϑ0 . Recovery of ϑk is trickier, but can be achieved using Proposition C.2 under the condition that L′ (ϑ0 ) ̸= 0. Proposition C.2. The coefficients ϑk can be computed using the recurrence relations  j    if j = k = 0, if j ≥ 0, k = 0,  ϑ0 1 ˜ dj,k = 0 if j = 0, k ≥ 1, if j = 0, k ≥ 1, dj,k = 0   Pk Pk−1 ˜ otherwise, otherwise. l=0 ϑl dj−1,k−l l=0 ϑl dj−1,k−l  ϑk = −



X X 1   cij d˜j,k +  L (ϑ0 ) p<e (i,j)∈A

′

max (i,j)∈A

i−j=emax

 ′ cij dj,k−emax +p   , L (ϑ) =

i−j=p

jcij ϑj−1 .

(i,j)∈A i−j=emax

Proof. Let dj,k denote the k-th power series coefficient of S(w)j so that S(w)j =

dj,k wk .

It is clear that dj,0 = ϑj0 and d0,k = 0 if k ≥ 1. Furthermore, X

dj,k =

(C.5)

ϑt1 · · · ϑtj .

t1 ,...,tj ≥0 t1 +···+tj =k

Next, we observe that since S(w)j = S(w) · S(w)j−1 , using the Cauchy product, ! ! ! k X X X X j k k S(w) = ϑk w dj−1,k w = ϑl dj−1,k−l wk , k

and so dj,k =

k X

l=0

(C.6)

ϑl dj−1,k−l .

l=0

To separate the effect of ϑk , let us consider the quantity dj,k in the case where ϑk = 0. Let d˜j,k be defined by Sk (w) =

k−1 X

ϑl w l ,

Sk (w)j =

l=0

X l

d˜j,l wl .

Similar reasoning implies that d˜0,0 = 1, d˜0,k = 0 for k ≥ 1 and d˜j,k =

k−1 X

ϑl d˜j−1,k−l .

(C.7)

dj,k − d˜j,k = jϑj−1 0 ϑk .

(C.8)

l=0

From (C.5), we find that

Now, recalling from the proof of Lemma C.1 that X cij wemax +j−i S(w)j = 0, (i,j)∈A

there is

X X

cij dj,k wemax +j−i+k = 0.

(i,j)∈A k

Rearranging gives X X

cij dj,k−emax +p wk = 0,

k p≤emax (i,j)∈A i−j=p

and therefore

cij dj,k−emax +p = 0.

(C.9)

p≤emax (i,j)∈A i−j=p

Substituting (C.8) into (C.9) gives X X cij d˜j,k + jϑj−1 ϑ k + 0

cij dj,k−emax +p = 0,

p<emax (i,j)∈A i−j=p

(i,j)∈A i−j=emax

and therefore

ϑk L′ (ϑ0 ) +

X (i,j)∈A i−j=emax

cij d˜j,k +

cij dj,k−emax +p = 0.

p<emax (i,j)∈A i−j=p

Rearranging for ϑk yields the result.

Appendix D

Fitting the Spectral Curve

In this section, we describe a general-purpose procedure for reconstructing non-physical branches of an algebraic Stieltjes transform from pointwise evaluations of the physical branch m1 . The basic idea is to fit an implicit algebraic relation P (z, m) = 0 as in (C.1) from samples {(zℓ , m1 (zℓ ))}nℓ=1 . Once such a relation is available, the candidate branch values at a given z ∈ C \ I are obtained as the s roots of the univariate polynomial m 7→ P (z, m), and non-physical branches can be accessed and tracked by analytic continuation away from the branch locus. Throughout, we assume that m is algebraic of degree s in the sense that there exists a nonzero bivariate polynomial P ∈ C[z, m], the ring of bivariate polynomials in (z, m) with complex coefficients, with degm P = s such that (C.1) holds for m = m1 on C \ I. Among all such relations, it is natural to seek a minimal one in the sense of having minimal degree in m, equivalently an irreducible P that does not contain spurious factors producing extraneous roots. Before turning to the numerical fitting procedure, we record structural constraints implied by the Stieltjes symmetry of m1 , starting with the fact that (C.1) admits a real-coefficient representative. Lemma D.1. Assume that P ∈ C[z, m] satisfies (C.1) for m = m1 and is irreducible in C[z, m]. Then there exists a nonzero scalar α ∈ C such that αP ∈ R[z, m]. In particular, after an irrelevant rescaling of P , one may assume P ∈ R[z, m]. 20

Proof. Since the distribution µ is real-valued, the Stieltjes transform m1 from (B.1) satisfies the Schwarz symmetry m1 (z̄) = m1 (z), z ∈ C \ I, where z̄ denotes the complex conjugate of z. Define the conjugate-reflected polynomial P ♯ ∈ C[z, m] by P ♯ (z, m) := P (z̄, m̄). Then, for every z ∈ C \ I, P ♯ (z, m1 (z)) = P (z̄, m1 (z)) = P (z̄, m1 (z̄)) = 0, so P ♯ also annihilates m1 identically on C \ I. By irreducibility, P and P ♯ must agree up to a nonzero scalar factor, hence there exists λ ∈ C \ {0} such that P ♯ = λP . Applying the ♯ operation twice yields P = (P ♯ )♯ = λ̄λP , so |λ| = 1. Choose α ∈ C \ {0} such that α/ᾱ = λ and set P̃ := αP . Then P̃ ♯ = ᾱP ♯ = ᾱλP = αP = P̃ , which implies that P̃ has real coefficients, that is, P̃ ∈ R[z, m]. Remark D.1. The rescaling in Lemma D.1 does not affect the roots of m 7→ P (z, m) for any fixed z, and has no impact on the reconstructed branches. A numerical normalization fixing the scalar ambiguity is described later in this section. △ Having established that (C.1) admits a real-coefficient representative, we henceforth assume P ∈ R[z, m]. To fit P from pointwise evaluations of m1 , we restrict P to a finite-dimensional model class by expanding it in a bivariate monomial basis. Concretely, we fix an index set A ⊂ N20 and parameterize X P (z, m) := cij z i mj , cij ∈ R. (D.1) (i,j)∈A

Let N := |A|, and collect the coefficients into a vector c ∈ RN according to a fixed ordering of A. A common choice for A is a rectangular support, A := {0, 1, . . . , dz } × {0, 1, . . . , s},

(D.2)

which corresponds to allowing degrees up to dz in z and up to s in m. More generally, the flexibility of A allows us to incorporate structural constraints on (C.1) such as moment constraints at infinity as in Appendix D.1 or reduced supports motivated by asymptotic cancellations.

D.1

Constraints at Infinity from Moments

Assuming that a finite number of moments of the target probability measure ν are available, the Laurent expansion of m1 at infinity induces linear constraints on the coefficients {cij }(i,j)∈A . Fix r ∈ N0 and define the moment vector µ := (µ0 , µ1 , . . . , µr ) ∈ Rr+1 . Definition D.1. For j ∈ N0 , define the j-fold discrete convolution µ∗j := µ ∗ · · · ∗ µ, | {z } j times

∗j with the convention µ∗0 = [1]. For q ∈ N0 , denote by µ∗j q the q-th entry of µ , equivalently

µ∗j q :=

j Y

µpℓ ,

p1 +···+pj =q ℓ=1 0≤pℓ ≤r ∗0 ∗j with the conventions µ∗0 0 = 1, µq = 0 for q ≥ 1, and µq = 0 whenever q < 0 or q > jr.

The following Proposition D.2 provides linear constraints on the coefficients using known moments. We note that Proposition D.2 is effectively a higher-order generalization of Lemma C.1. 21

Proposition D.2. Let P be parameterized as in (D.1) for a fixed choice of index set A ⊂ N20 , and assume that (C.1) holds for m = m1 . Suppose that the moments {µp }rp=0 are available, and let emax := max{i − j : (i, j) ∈ A, cij ̸= 0}.

(D.3)

Then the coefficients {cij }(i,j)∈A satisfy the linear constraints X

(−1)j cij µ∗j i−j−emax +l = 0,

l = 0, 1, . . . , r.

(D.4)

(i,j)∈A

Proof. Using the identity (x − z)−1 = −z −1 the truncated Laurent expansion m1 (z) = − Set ϕr (z) :=

p=0 µp z

−p

p=0 (x/z)

r X µp p=0 z

p+1 + O z

+ O z −r−2 in (B.1) and applying (C.3), we obtain

−r−2

|z| → ∞.

, so that m1 (z) = −z −1 ϕr (z) + O(z −r−2 ). For each j ∈ N0 , m1 (z)j = (−1)j z −j ϕr (z)j + O z −(j+r+1) ,

|z| → ∞.

Expanding the multinomial ϕr (z)j and collecting coefficients gives ϕr (z)j =

jr X

−q µ∗j , q z

q=0

where µ∗j q is the q-th entry of the j-fold discrete convolution of µ, as stated in Definition D.1. Consequently, for each (i, j) ∈ A, z i m1 (z)j = (−1)j

jr X

i−j−q µ∗j + O z i−j−(r+1) , q z

|z| → ∞.

q=0

Substituting this into (D.1) and canceling the coefficients of z emax −l for l = 0, 1, . . . , r in the truncated expansion of P (z, m1 (z)) yields (D.4). Remark D.2. In practice, while higher moments of an empirical distribution may be too noisy to estimate reliably, the leading constraint l = 0 in (D.4) is always available since µ0 = 1 holds unconditionally. This constraint enforces cancellation on the top (i − j)-edge of A, X (−1)j cij = 0. (i,j)∈A i−j=emax

For a full rectangular support (D.2), the top (i − j)-edge consists of the single monomial (i, j) = (dz , 0), and the above condition reduces to cdz ,0 = 0. It is therefore natural to exclude this monomial from A a priori. More refined choices of A, including slanted Newton polygon supports, are discussed later in connection with branch-point constraints. △

D.2

Constraints from Branch Points

Sk We now exploit the prior knowledge of the branch-cut geometry I = j=1 [aj , bj ]. As discussed in Appendix C, the branch locus of the algebraic relation (C.1) is identified by the zeros of the discriminant ∆(z), obtained by eliminating m∗ from (C.2). In the present setting, we assume that the spectral edges {aj , bj }kj=1 are among these branch points, which yields the nonlinear endpoint constraints ∆(aj ) = 0

and

∆(bj ) = 0, 22

j = 1, . . . , k.

(D.5)

Since ∆ ∈ R[z] is a univariate polynomial, each condition in (D.5) implies that (z − aj ) and (z − bj ) divide Qk ∆. If the endpoints {aj , bj }kj=1 are distinct, then with π(z) := j=1 (z − aj )(z − bj ) we obtain the divisibility condition π | ∆. This divisibility viewpoint suggests an alternative enforcement strategy, for instance by ˜ parameterizing ∆(z) = π(z)∆(z), although in our implementation we impose the endpoint conditions (D.5) directly. The endpoint constraints (D.5) also imply a simple Ps necessary richness condition on the polynomial model class. Recall s = degm (P ) and write P (z, m) = j=0 aj (z)mj with degz (aj ) ≤ dz . Then degz (∆) ≤ (2s − 1)dz . Indeed, up to a nonzero factor depending on the leading coefficient as (z), one has ∆(z) = Resm (P, ∂m P ), and the resultant is the determinant of the (2s − 1) × (2s − 1) Sylvester matrix in the variable m. Since all Sylvester entries are coefficients of P and ∂m P and satisfy degz (P ) = degz (∂m P ) ≤ dz , it follows that degz (∆) ≤ (2s − 1)dz . See, for example, [46, Section 4.1]. From the fact that a nonzero polynomial with 2k distinct real zeros must have degree at least 2k, this yields the necessary condition (2s − 1)dz ≥ 2k, (D.6) for P to accommodate 2k prescribed distinct endpoints. We emphasize that (D.6) is only a lower bound; the algebraic model may admit additional branch points away from {aj , bj }kj=1 , contributing further zeros of ∆.

D.3

Numerical Reconstruction from Samples

We now describe the numerical procedure used to fit an implicit relation P (z, m) = 0 from samples of the physical branch m1 , and to recover non-physical branches by root finding and continuation. Throughout, we fix an index set A ⊂ N20 and parameterize P by (D.1) with coefficient vector c ∈ RN , N := |A|. A convenient choice of A is the rectangular set (D.2), where the sheet count s is set empirically, while dz can be initialized using the lower bound (D.6) and increased as needed. In practice, we increase (s, dz ) until the fitting error stabilizes while avoiding over-parameterization. Sampling. We choose sampling points {zℓ }nℓ=1 ⊂ C \ I along multiple contours enclosing each Ij and the overall I, but away from the branch locus, in particular separated from the endpoints {aj , bj }kj=1 by a small buffer. At each sample, we evaluate m̂ℓ := m1 (zℓ ). The details of our approach for evaluating m1 from empirical measures are described in Appendix E. To preserve the Schwarz symmetry needed for real-coefficient structure as discussed in Lemma D.1, we sample in conjugate pairs: whenever zℓ is included, we also include z̄ℓ , which yields m1 (z̄ℓ ) = m1 (zℓ ). Linear system for the coefficients.

Define the design matrix A ∈ Cn×N with entries

Aℓ,(i,j) := zℓi m̂jℓ ,

(i, j) ∈ A,

ℓ = 1, . . . , n,

(D.7)

where the column index (i, j) follows a fixed ordering of A. The samples satisfy P (zℓ , m̂ℓ ) ≈ 0, hence Ac ≈ 0. We compute c as a nontrivial approximate null vector via the augmented homogeneous least-squares problem ĉ := argmin∥Aη c∥2 , c∈R

where

subject to

∥c∥2 = 1,

(D.8)

  ℜ(A) Aη := ℑ(A) , √ ηI

η ≥ 0,

(D.9)

and I ∈ RN ×N is the identity matrix. The case η = 0 recovers the unregularized fit, while η > 0 provides mild Tikhonov-type stabilization on noisy data. The solution is given by the right singular vector corresponding to the smallest singular value of Aη . The scalar ambiguity of P (see Remark D.1) is fixed by the normalization ∥c∥2 = 1 and a deterministic sign convention, for example by requiring the first nonzero entry of c to be positive. 23

Incorporating moment constraints. When moment constraints derived in Proposition D.2 are used, they can be written as a linear system Bc = 0, where B ∈ R(r+1)×N encodes the r + 1 enforced constraints {µp }rp=0 . To enforce these constraints, we restrict c to ker(B). Concretely, let Q ∈ RN ×Nc have orthonormal columns spanning ker(B), obtained for example from an SVD of B. Then, any feasible coefficient vector can be written as c = Qy. We then solve the reduced problem ŷ := argmin∥Ãη y∥2 y∈R

where

subject to

∥y∥2 = 1,

(D.10)

  ℜ(AQ) Ãη := ℑ(AQ) , √ ηI

followed by ĉ := Qŷ. In practice, enforcing only the leading constraint µ0 = 1 (that is, r = 0) is often beneficial even when higher-order moments are unreliable. Algorithm D.1 summarizes the coefficient-estimation step (with optional moment constraints), producing a fitted polynomial P̂ .

Appendix E

Computing the Stieltjes Transform and Density

This section describes the numerical procedure used to evaluate the Stieltjes transform and to recover the density from it. The first task is to evaluate the physical Stieltjes branch m(z). When empirical eigenvalues are available, this can be done directly by summing the Cauchy kernel. However, after fitting an algebraic spectral curve P (z, m) = 0, evaluating m(z) becomes a root-selection problem: for each query point z, one must identify the physical root among all algebraic roots of P (z, ·). This process is described in Appendix E.1. The second task, discussed in Appendix E.2, is the inverse Stieltjes problem: once m(z) has been evaluated near the real axis, we recover the density ρ from its boundary values.

E.1

Physical Root Selection for the Stieltjes Transform

Recall from (B.1) that the Stieltjes transform of a spectral density ρ is Z ρ(x) m(z) = dx, z ∈ C \ I. x −z R

(E.1)

For empirical eigenvalues {λ1 , . . . , λn }, this quantity is computed directly as n

mn (z) =

1X 1 . n i=1 λi − z

The more difficult problem arises when we evaluate the Stieltjes transform through a fitted algebraic relation P (z, m) = 0. For each fixed query point z ∈ C \ I, the polynomial m 7→ P (z, m) has s = degm (P ) roots. Only one of these roots is the physical Stieltjes value, while the remaining roots belong to non-physical sheets of the algebraic curve. Thus, evaluating m(z) from P is a nontrivial branch-selection problem, especially for high-degree curves and multi-bulk spectra, where roots may approach one another, exchange ordering, or pass near branch points. Figure E.1 illustrates this issue for the diffusion model example; see Appendix I.5. In the left panel, we evaluate all roots of m 7→ P (z, m) at z = λ + iδ, for a small fixed δ > 0 and a broad range of λ, and plot the candidate values π −1 ℑmi (λ + iδ), i = 1, . . . , s. Since the roots are unordered from one value of λ to the next, the result is an unordered scatter of algebraic candidates rather than a connected curve. At sufficiently fine resolution, the neighboring points visually trace several candidate branches. Only one of these traces corresponds to the physical density, 1 ρ(λ) ≈ ℑm(λ + iδ), π 24

Algorithm D.1: Implicit fitting of an algebraic relation P (z, m) = 0 from samples of m1 : Index set A ⊂ N20 such as in (D.2) with N := |A|; Oracle for the physical branch m1 (z) on C \ I (see Algorithm E.1); Number of sampling points n (even); Stabilization parameter η ≥ 0; Optional moments {µp }rp=0 (to enforce Proposition D.2). P Output : Coefficient vector ĉ ∈ RN defining P̂ (z, m) = (i,j)∈A ĉij z i mj ; Fitting residual error e := ∥Aĉ∥2 /∥A∥F for model selection.

Input

// Sample the physical branch away from the cut, in conjugate pairs 1 Initialize set of points Z ← {} and function values M ← {}. n/2

2 Generate n/2 points {zℓ }ℓ=1 ⊂ C

on arbitrary contours enclosing I, excluding ϵ-neighborhood from I. 3 Evaluate m̂ℓ ← m1 (zℓ ). // Total of n samples with conjugate points 4 Append (zℓ , m̂ℓ ) and (z̄ℓ , m̂ℓ ) to (Z, M). // Build the design matrix (D.7) from samples 5 Form A ∈ C

n×N

with columns indexed by (i, j) ∈ A: Aℓ,(i,j) ← zℓi m̂ℓj ,

(zℓ , m̂ℓ ) ∈ (Z, M).

// Optional: enforce linear moment constraints Bc = 0 from Proposition D.2 r

6 if moments {µp }p=0 are provided then 7 8 9

(r+1)×N

Build B ∈ R from (D.4) Compute an orthonormal basis Q ∈ RN ×Nc for ker(B) (e.g. using SVD) Replace A ← AQ and interpret unknowns as y ∈ RNc

// Nc = N − rank(B)

// Incorporate Tikhonov stabilization

  ℜ(A) 10 Form the augmented matrix Aη ← ℑ(A) in (D.9) √ ηI

// I is the identity of matching size

// Solve the homogeneous least-squares problem min∥Aη c∥2 in (D.8) or in (D.10) if constraints exit 11 UΣV ← Aη // Perform thin SVD where Σ = diag(σ1 , . . . , σN ) and σ1 ≥ · · · ≥ σN 12 v̂ ← V[:,N ] // The last column of V, i.e., the right singular vector for the smallest singular value σN 13 if moment constraints were enforced then

ĉ ← Qv̂ else ĉ ← v̂

14 ĉ ← ĉ/∥ĉ∥2

// Normalize, optionally, fix a deterministic sign convention

15 return ĉ and residual error e := ∥Aĉ∥2 /∥A∥F .

namely the trace that agrees with the empirical histogram shown in the background. The right panel shows the result after continuously selecting the physical branch. Thus, the purpose of the procedure below is to pass from the unordered scatter of algebraic roots in the left panel to the coherent physical curve in the right panel. The physical branch is characterized by the Herglotz property and by its normalization at infinity, ℑm(z) > 0

for z ∈ C+ ,

m(z) ∼ −

1 z

as |z| → ∞.

However, these conditions alone are not always sufficient for stable pointwise root selection. Instead, we treat the fitted relation as a spectral curve n o C = (z, m) ∈ C2 | P (z, m) = 0 , 25

(a) All Polynomial Roots

104

ρ(λ) = π1 =m1 (λ + δi)

1 =mi (λ + δi) π

104 10−1 10−6

10−11 10−16 10−5

10−1

10−1 10−6

10−11

Roots Empirical Histogram

10−3

(b) Density from Physical Root

101

103

Physical density Empirical histogram

10−16 10−5

10−3

10−1

101

103

Figure E.1: Physical branch selection for the fitted algebraic Stieltjes transform in the diffusion-model example at the initial size n = 22 K. Panel (a) shows all candidate densities π −1 ℑmi (λ + iδ) with δ = 10−6 , obtained from the seven roots mi of P (λ + iδ, m) = 0, where degm (P ) = 7. These unordered candidates include both physical and non-physical sheets and illustrate the difficulty of pointwise root selection in a multi-bulk, log-scale spectrum. Panel (b) shows that the continuation procedure successfully selects the physical branch, producing a coherent density curve that matches the empirical histogram on the bulks. The remaining sloped baseline in the spectral gaps is a finite-δ artifact and is addressed in Appendix E.2. and track the physical sheet geometrically on this curve from an anchor point to the desired target point. The anchor point za is chosen so that the physical root ma can be identified reliably. One option is to place za far from the spectrum, where the normalization m(z) ∼ −1/z selects the correct algebraic root. Another option is to place za closer to the real axis and select the algebraic root closest to the empirical value mn (za ) from (E.1). The latter approach is often more stable for complicated algebraic curves. Starting from (za , ma ) ∈ C, we choose a continuous path γ : [0, 1] → C+ with γ(0) = za and γ(1) = z, and lift this path to the spectral curve by seeking m(s) such that P (γ(s), m(s)) = 0,

m(0) = ma .

Root ordering alone is unreliable along such paths: near branch points, two sheets may come close, exchange order, or locally appear to cross. To avoid jumping between sheets, we follow the local geometry of the curve using a predictor–corrector procedure. Differentiating P (γ(s), m(s)) = 0 along the path gives ∂ P (γ(s), m(s)) ′ dm =− z γ (s). ds ∂m P (γ(s), m(s)) Equivalently, away from branch points, the local tangent satisfies dm ∂ P (z, m) =− z . dz ∂m P (z, m) This tangent relation provides a predictor for the next point znew along the path on the same sheet. The predicted value is then corrected by enforcing P (znew , mnew ) = 0, which amounts to projecting the predicted point back onto the spectral curve. This correction is performed by Newton iteration, viewing mnew ∈ C as a two-dimensional real variable. The selected root at znew is the algebraic root closest to this corrected value. If no root is sufficiently close, or if the local geometry is too stiff, the continuation step is subdivided and repeated with smaller increments until the root closeness criterion is satisfied. Repeating this predictor–corrector continuation along the path yields a continuous selection of the physical branch. This continuation-based, or homotopic, root selection is what allows the fitted polynomial to be used as a numerical Stieltjes transform. The resulting physical-branch selection procedure is summarized in Algorithm E.1. 26

Algorithm E.1: Physical branch selection for an algebraic Stieltjes transform Input

: Polynomial P ∈ R[z, m] (see Algorithm D.1); Nz ⊂ C+ ; Query points {zj }j=1 Anchor rule for selecting ma at za ; Newton tolerance ϵ > 0.

Output : Physical Stieltjes values m(zj ) satisfying P (zj , m(zj )) = 0. 1 for j = 1, . . . , Nz do 2 3 4

Choose an anchor point za and a path γ : [0, 1] → C+ with γ(0) = za , γ(1) = zj . Select the physical root ma at za using either m(z) ∼ −z −1 or the empirical value mn (za ). Initialize m(0) ← ma .

// Lift the path γ to the spectral curve P (z, m) = 0 5 6

for successive points zold , znew along γ do Compute the tangent predictor ṁ = −

7 8

∂z P (zold , m) ż. ∂m P (zold , m)

Predict mpred at znew . Correct mpred by Newton iteration on P (znew , m) = 0.

9 10 11 12 13

if Newton correction is accepted then Set m ← mnew . else Subdivide the path segment and repeat.

// Avoid sheet jumping

Store m(zj ) ← m. N

z 14 return {m(zj )}j=1 .

The continuation step presented here is not merely a numerical convenience. Without it, direct pointwise root selection is unstable even for relatively simple algebraic curves, because root orderings can change from one query point to the next. This instability becomes much more severe for high-degree, multi-bulk, or log-scale spectra, where several sheets may come close over many orders of magnitude. The diffusion-model example in Figure E.1 is deliberately chosen as a stringent stress test of this type: the fitted curve has degree degm (P ) = 7, the density spans many decades in both λ and ρ(λ), the geometry is highly stiff near the spectral edges, and the support contains several strongly unbalanced bulks, with the leftmost bulk orders of magnitude taller than the others. Nevertheless, the geometric continuation procedure reliably resolves the physical branch and produces a coherent density from the fitted algebraic curve. With the physical root m(z) selected consistently over a query grid z = λ + iδ, the density ρ(λ) can then be recovered from its boundary values using the inverse Stieltjes procedure described next.

E.2

Density Recovery and δ-Extrapolation

Once the physical branch m(z) has been selected, the density is recovered from the boundary values of the Stieltjes transform. Recall the Plemelj formula from (B.2): ρ(λ) =

1 lim ℑm(λ + iδ). π δ→0+

104

(a) Without Extrapolation

104 10−1

ρ(λ)

10−1 10−6 10−11 10−16 10−5

(b) With Extrapolation

Evaluated density Poisson kernel δ floor Empirical histogram

10−3

10−1

10−6 10−11

101

10−16 10−5

103

Evaluated density Empirical histogram

10−3

10−1

101

103

Figure E.2: Inverse Stieltjes recovery for the diffusion-model spectrum at the initial size n = 22 K. Panel (a) shows the direct finite-δ density ρδ (λ) = π −1 ℑm(λ + iδ), which agrees with the empirical histogram on the bulks but exhibits a Poisson δ-floor in spectral gaps. Panel (b) shows the density after extrapolating the finite-δ evaluations to δ = 0, removing the artificial gap baseline while preserving the bulk densities. In numerical computations, however, one cannot evaluate exactly on the real axis at δ = 0. Instead, one evaluates m at a small positive offset δ, giving the regularized density ρδ (λ) :=

1 ℑm(λ + iδ). π

The parameter δ plays the role of a smoothing scale: larger values of δ improve numerical stability but smear sharp features of the density, while smaller values give sharper reconstructions at the cost of increased sensitivity near branch points and spectral edges. This smoothing effect has a precise interpretation. Indeed, Z 1 1 δ ρ(ξ) dξ = (Pδ ∗ ρ)(λ), ρδ (λ) = ℑm(λ + iδ) = 2 π π (λ − ξ) + δ2 R where Pδ (λ) :=

1 δ 2 π λ + δ2

is the Poisson kernel in the upper half-plane [45]. Thus, evaluating the inverse Stieltjes formula at a single finite δ amounts to convolving ρ with Pδ , which introduces a finite-δ blur and a nonzero baseline outside the support. This effect is especially visible for spectra plotted on logarithmic scales. In spectral gaps, where the limiting density should vanish, the finite-δ Poisson tail creates a nonzero baseline. Away from the support, this baseline behaves like Pδ (λ) ∼ δ/(πλ2 ), which appears as a line of slope −2 on a log–log plot. We refer to this artifact as the Poisson δ-floor. Decreasing δ lowers this floor and makes the bulks stand out more sharply, but it also makes branch tracking and root selection more sensitive. The left panel of Figure E.2 illustrates this issue for the diffusion-model example. In the right panel, the baseline is eliminated by the extrapolation procedure described below, and the bulk edges are sharply resolved. To remove this artifact, we evaluate the physical branch on a ladder of offsets δ0 < δ1 < · · · < δr and extrapolate the sequence ρδi (λ) to δ = 0 at each fixed λ. A simple approach is to use a low-degree polynomial fit in δ. In practice, we use a Chebyshev basis for this extrapolation, which is more stable than a monomial basis when several δ-levels are used [47]. Concretely, after mapping the offset ladder {δi }ri=0 to points {ηi }ri=0 ⊂ [−1, 1], we approximate ρδi (λ) ≈

q X

ck (λ)Tk (ηi ),

k=0

where Tk are Chebyshev polynomials. The extrapolated density is then obtained by evaluating this fitted polynomial at the value of η corresponding to δ = 0. This approach is particularly useful for log-scale spectra, such as the diffusion-model example shown in Figure E.2. In this example, we used four δ-levels and a quadratic Chebyshev extrapolation. This was sufficient to reduce the artificial Poisson floor below 10−20 across the spectral gaps, without additional regularization. For densities represented on a uniform linear grid, one can alternatively exploit the convolution identity ρδ = Pδ ∗ ρ. In that case, the collection of finite-δ evaluations can be viewed as a regularized deconvolution problem, which can be solved efficiently in Fourier space. Indeed, on a uniform linear grid, the convolution relation gives F{ρδ }(ω) = F {Pδ }(ω) F{ρ}(ω) = e−δ|ω| F{ρ}(ω), where F is the Fourier transform. Thus, if ρδ were known exactly on a uniform grid, the formal inverse would be h i ρ = F −1 eδ|ω| F{ρδ }(ω) . The factor eδ|ω| inverts the Poisson smoothing, but it also amplifies high-frequency numerical noise. Hence, the Fourier recovery must be regularized, for instance by damping or truncating high-frequency modes. The δ-extrapolation is similarly important for densities with wide log–log scale spectra, such as the diffusion-model example shown in Figure E.2. Without it, the finite-δ Poisson floor obscures spectral gaps and makes narrow bulks appear artificially connected on a log–log plot. After extrapolation, the density retains the same bulk-scale agreement with the empirical histogram, while the gap baseline is pushed below the plotted range. This allows the multi-bulk structure and spectral edges to be resolved cleanly, even in this stringent log-scale test case.

Appendix F

Free Decompression on Algebraic Spectral Curves

In this section, we develop the mathematical and computational foundation of free decompression from an algebraic Stieltjes relation. We begin in Appendix F.1 with the Nica–Speicher theorem, which identifies the R-transform scaling that underlies free compression and motivates its inverse operation. We then define infinite decompressibility in Appendix F.2 and show its equivalence with free infinite divisibility. The next issue is computational. Although free decompression leads naturally to a nonlinear PDE for the Stieltjes transform, Appendix F.3 explains why solving this PDE directly is numerically unstable in realistic examples, especially when analytic continuation across branch cuts is required. The main purpose of the remainder of the section is to show that this difficulty can be avoided when the initial Stieltjes transform is represented by an algebraic relation P (z, m) = 0. In Appendix F.4, free decompression is lifted to a rational transformation of the associated spectral curve. In Appendix F.5, this transformation is converted into an evolved polynomial relation whose roots give all algebraic candidates at later times. The remaining challenge is to select the physical Stieltjes branch from these candidates. We address this in Appendix F.6 by tracking the physical sheet geometrically along the free-decompression flow. This leads to a predictor–corrector method on the spectral curve, formulated as a 2 × 2 Newton system in the natural curve coordinates. This viewpoint is the basis of the practical free-decompression algorithm used in the numerical experiments.

F.1

The Nica–Speicher Theorem

Our exposition begins with a complete proof of the R-transform relationship which underpins the free decompression procedure, also known as the Nica–Speicher Theorem. Although our conditions appear to be less general than required in practice, our emphasis here is on establishing a rigorous set of assumptions upon which free decompression can be guaranteed to hold. To do so, it is important to define several important objects pertaining to free probability theory. First, we begin with the definition of “limits” of matrices in the noncommutative algebraic sense. Definition F.1 (Free Algebra Limit). Let trn denote the normalized trace, given by trn (A) = n1 tr(A). A sequence of matrices An has a limit a if for every integer k, trn (Akn ) converges as n → ∞ to a value, denoted 29

∞ tr∞ (ak ). More generally, a family of sequences of matrices ({A1,n }∞ n=1 , . . . , {Ak,n }n=1 ) has limits (a1 , . . . , ak ) if for any (noncommutative) polynomial combination p(x1 , . . . , xk ), trn (p(A1,n , . . . , Ak,n )) converges to a value, denoted tr∞ (p(a1 , . . . , ak )).

Next, we consider the critical definition of asymptotic freeness, first introduced in Voiculescu [50]. ∞ Definition F.2 (Asymptotically Free). A family of sequences of matrices ({A1,n }∞ n=1 , . . . , {Ak,n }n=1 ) is ∞ ∞ asymptotically free from another family of sequences of matrices ({B1,n }n=1 , . . . , {Bl,n }n=1 ) if the two families have a joint limit (a1 , . . . , ak , b1 , . . . , bl ) and (a1 , . . . , ak ), (b1 , . . . , bl ) are freely independent in the following sense: for any integer m ≥ 1 and any “alternating” sequence of elements x1 , x2 , . . . , xm such that x2i+1 is a noncommutative polynomial combination of (a1 , . . . , ak ) and x2i is a noncommutative polynomial combination of (b1 , . . . , bl ), if tr∞ (xi ) = 0 for all i, then tr∞ (x1 x2 · · · xm ) = 0.

Finally, we will require another important transform, similar to the R-transform. R Definition F.3 (S-Transform). Let ν be a probability measure and let µn = xn dν(x) denote its moments for each integer n. The truncated moment-generating function of ν is given by M (z) =

∞ X

µn z n .

n=1

If χ(z) is the functional inverse of M so that M (χ(z)) = z, then the S-transform of µ is given by S(z) =

z+1 χ(z). z

The significance of the S-transform is its behavior for multiplicative free convolution. If a and b are freely independent, then the spectral distribution of b1/2 ab1/2 has S-transform Sab satisfying [13, Lemma 2.19] Sab (z) = Sa (z)Sb (z), where Sa and Sb are the S-transforms of a and b respectively. To relate the R- and S-transforms, we have the following lemma. Lemma F.1. Letting R(z) and S(z) denote the R- and S-transforms of a probability measure, if y = zS(z), then z = yR(y). Proof. From Couillet & Liao [13, Definition 5], if m(z) is the Stieltjes transform, then z+1 m = −zS(z). zS(z) Substituting y = zS(z) gives

z+1 y

= −y.

But since m(R(y) + y −1 ) = −y from Couillet & Liao [13, Definition 5], z/y = R(y) and so z = yR(y). Now we may prove the Nica–Speicher Theorem [30]. Theorem F.2 (Nica–Speicher). Let {An }∞ n=1 be a sequence of n × n deterministic real symmetric matrices whose empirical spectral distribution converges weakly as n → ∞ to a compactly supported probability measure kn with R-transform R(z) and S-transform S(z). Consider a sequence of integers {kn }∞ n=1 such that limn→∞ n = τ ∈ (0, 1). For each n, let Qn be a Haar-distributed random orthogonal matrix of size n × n and define kn Ãn = Q⊺n An Qn . Let σn be a random permutation of {1, . . . , n} and let Bn = (Ãσn (i)σn (j) )i,j=1 . Then as n → ∞, the empirical spectral distribution of Bn converges to a limiting probability measure whose R- and S-transforms Rτ (z) and Sτ (z) are given by Rτ (z) = R(τ z),

Sτ (z) = S(τ z). 30

Proof. Fix an integer n. For a permutation σ of {1, . . . , n} into itself, we may define the permutation matrix Πσ such that (Πσ )ij = 1 if i = σ(j) and 0 otherwise. In this way, Bn becomes the top-left kn × kn block of ˜ = Π⊺ Ã Π = Π⊺ Q⊺ A Q Π . Ã n σn n σn σn n n n σn Since Q̃n = Qn Πσn is also an orthogonal matrix and the Haar measure is invariant under permutations, D D ˜ = Ã . Consequently, we can assume without loss of generality that the permutations Q̃ = Q and so Ã n

σn are the identity and consider

B̃n = Pn Ãn Pn ,

where Pn is a diagonal matrix with (Pn )ii = 1 if i ≤ kn and 0 otherwise. Because Qn is orthogonal, the conjugation Ãn = Q⊺n An Qn leaves the eigenvalues invariant. Therefore, the empirical spectral distribution of Ãn is identical to that of An for all n, and converges weakly to the same compactly supported probability measure ν, which is uniquely characterised by its R-transform R(z). Furthermore, Pn clearly has an empirical spectral distribution converging weakly to a Bernoulli distribution: lim νPn = νp = (1 − τ )δ0 + τ δ1 .

n→∞

Let trn (A) = n1 tr(A) denote the normalize trace. We note that limn→∞ trn (Pn ) = limn→∞ knn = τ . Now, ∞ since An and Pn both have spectral distributions converging weakly, the family ({An }∞ n=1 , {Pn }n=1 ) has a ∞ ∗ ∞ joint limit (a, p). From Voiculescu’s Theorem [50, Theorem 3.8], ({Qn }n=1 , {Qn }n=1 ) has a limit (u, u∗ ) and ∞ ∞ ∗ ∞ ∞ ∞ ({An }∞ n=1 , {Pn }n=1 ) and {Qn }n=1 , {Qn }n=1 ) are asymptotically free. To show that {Ãn }n=1 and {Pn }n=1 ∗ are asymptotically free, let pi (u au) and qi (p) be centered polynomials such that tr∞ (pi ) = tr∞ (qi ) = 0. Note that pi (u∗ au) = u∗ pi (a)u since u∗ u is the identity. Now since tr∞ (xy) = tr∞ (yx), tr∞ (p1 (u∗ au)q1 (p)p2 (u∗ au)q2 (p) · · · ) = tr∞ (u∗ p1 (a)uq1 (p)u∗ p2 (a)uq2 (p) · · · ) = tr∞ (p1 (a)uq1 (p)u∗ p2 (a)uq2 (p) · · · u∗ ) = 0, since it is an alternating sequence of elements belonging to {a, p} and {u, u∗ }. Since pi and qi were arbitrary, ã = u∗ au and p are freely independent. To identify the distribution of b̃ = pãp, we note that p1/2 = p and so the S-transform of b̃, denoted Sb̃ satisfies Sb̃ = Sã Sp = Sa Sp . Therefore, we will need to find the S-transform of p. R To do P so, note that mn = xn νp (dx) = τ for any n. Therefore, the truncated moment-generating function ∞ τz z Mp (z) = n=1 τ z n = 1−z . The inverse is given by χp (z) = z+τ and so Sp (z) =

z+1 z+1 χ(z) = . z z+τ

Consequently, Sb̃ (z) =

z+1 S (z). z+τ a

(F.1)

Now observe that B̃n is a block matrix, with the top-left kn × kn block comprised of Bn , and every other elements is zero. Consequently, tr(Bkn ) = tr(B̃kn ) and so trn (B̃kn ) = n1 tr(Bkn ) = knn trn (Bkn ). The limit of B̃n is b̃ = pãp and so Bn has a limit b satisfying tr∞ (bk ) = Therefore, Mb̃ (z) =

∞ X

1 tr (b̃k ). τ ∞

tr∞ (b̃k )z k = τ

k=1

∞ X

tr∞ (bk )z k = τ Mb (z).

k=1

For the inverses, since z = Mb̃ (χb̃ (z)) = τ Mb (χb̃ (z)), we find that z χb̃ (z) = χb . τ 31

Now observe that z + 1 z Sb̃ (z) = χb z τ z z z z + τ z τ +1 Sb = z χb = , χb τ τ z τ τ and so by equating the two through χb , Sb̃ (z) =

z + 1 z . S z+τ b τ

Equating with (F.1) gives Sb (z) = Sa (τ z), completing the first part of the result. Now we need to express this relationship in terms of the R-transform. To do so, recall from Lemma F.1 that if y = zSb (z), then z = yRb (y). Since Sb (z) = Sa (τ z), y = zSb (z) = zSa (τ z) implies that τ y = (τ z)Sa (τ z). This immediately implies that

τ z = (τ y)Ra (τ y).

But since y = zSb (z) also implies that z = yRb (y), there is τ yRb (y) = τ yRa (τ y). Therefore, finally, Rb (y) = Ra (τ y) = R(τ y), completing the proof.

F.2

Infinite Decompressibility

The Nica–Speicher Theorem guarantees that the distribution obtained under free compression is well-defined through the R-transform by taking τ < 1. However, it does not imply that the inverse operation, taking τ > 1, is generally well-defined. Indeed, there are probability measures where free decompression need not be valid. In this section, we define the class of probability measures for which free decompression can be applied without bound on τ . This motivates the following definition. Definition F.4 (Infinitely Decompressible Measure). A probability measure ν with R-transform R is infinitely decompressible if for every τ ≥ 1, there exists a probability measure ντ with R-transform satisfying Rντ (z) = R(τ z). Since the R transform admits the series expansion X Rντ (w) = rk (ντ )wk−1 , k≥1

where rk (ντ ) is the k-th free cumulant of ντ , the Nica–Speicher Theorem implies rk (ντ ) = τ k−1 rk (ν),

k ≥ 1.

Let Dc ν denote the dilation of a measure ν, defined as the law of cX when X ∼ ν. This is also written as c∗ ν, as this operation is the push-forward of ν induced by the map x 7→ cx. Also, denote by ν ⊞s the s-fold free additive convolution power. The R-transform of dilation and convolution power operations satisfies the identities RDc ν (w) = c Rν (cw),

(F.2a)

Rν ⊞s (w) = s Rν (w),

(F.2b)

rn (Dc ν) = cn rn (ν), rn ν ⊞s = srn (ν).

(F.3a)

or, in terms of cumulants: (F.3b)

From the above, it is immediate that dilation and power operations commute. The fractional powers s in the above relations are due to Nica & Speicher [30, Proposition 1.2] (see also [43, Proposition 1.2]). 32

Definition F.5. The measure ν ∈ P(R) is freely infinitely divisible (FID, also called ⊞-infinitely divisible), if for any n ∈ N, there exists νn ∈ P(R) such that ν = νn ⊞ · · · ⊞ νn = νn⊞n . | {z } n times

A necessary and sufficient condition for a measure to be FID is that its R transform should have holomorphic extension defined on C+ with values in C− \ R [5, Theorem 5.10]. The condition that ν is FID is closely related to free decompression through the following proposition. Proposition F.3 (Infinitely Decompressible ⇐⇒ FID). A probability law ν is infinitely decompressible if and only if it is freely infinitely divisible. Proof. (Sufficiency) If ν is FID, then for each τ ≥ 1, the dilation Dτ ν is also FID, hence the free convolution ⊞τ

power ντ := (Dτ ν)

−1

exists. Using (F.2a) and (F.2b), we obtain Rντ (w) = R

⊞τ

(Dτ ν)

−1

(w) = τ −1 RDτ ν (w) = Rν (τ w),

so the free decompression flow exists. (Necessity) Assume ν is decompressible. Fix n ∈ N and set τ = n, yielding a measure νn such that Rνn (w) = Rν (nw). Let σn := Dn−1 νn . Then Rσn (w) = n−1 Rνn (n−1 w) = n−1 Rν (w), hence Rσ⊞n (w) = n

Rν (w) and thus ν = σn⊞n . Since this holds for all n, ν is FID. We also note that any FID distribution can be approximated by compound free Poisson distributions [44, p. 117] given in Appendix I.1.

F.3

The PDE Approach and Challenges

Let ν be a probability measure and let {ντ }τ ≥1 be a family of probability measures obtained by free decompression such that for each τ , denoting by Rτ the R-transform of ντ , Rτ (z) = R(τ z), where R is the R-transform of ν. For each τ ≥ 1, let mτ denote the Stieltjes transform corresponding to ντ . Ameli et al. [1] proved that m(t, z) := met (z) satisfies the PDE ∂m(t, z) ∂m(t, z) = −m(t, z) + m(t, z)−1 . ∂t ∂z

(F.4)

In the following exposition, the time change τ = et will be used repeatedly with reference to the flow (F.4). While free decompression is always achievable in principle via (F.4), in practice, it is extraordinarily difficult to perform for several reasons. The primary difficulty lies in the ill-posedness of (F.4), which can be seen in two ways. The first is through the relationship between (F.4) and the Burgers equation: if τ = et − 1 and u(τ, z) = −t −e /m(τ, z), then u satisfies the inviscid complex Burgers equation ∂τ u + u∂z u = 0. This equation is well-studied in the setting where m is real-valued, including in related problems in RMT [7, 8, 9], and a number of techniques have been developed to obtain accurate approximations of solutions. Arguably the most important consideration is the development of shocks, where after some time t, multiple solutions emerge. These shocks present fatal instabilities when solving the PDE directly. To deal with this issue, one often introduces regularization through the viscous complex Burgers equation ∂τ u + u∂z u = ξ∂zz u, taking ξ to be small [42]. This equation has a unique explicit solution, called the Cole-Hopf solution, and taking ξ → 0+ yields the entropy solution to the inviscid Burgers equation. In the real-valued case, the entropy solution is given by the Hopf-Lax formula. However, the current setting requires that m is complex valued, and methods developed for the real-valued setting [42] cannot be readily deployed. When performing free decompression, shocks inevitably emerge, and while the Cole-Hopf solution can be computed, it need not converge as ξ → 0+ to the solution where m(τ, z) remains a valid Stieltjes transform. To our knowledge, no examination of the complex inviscid Burgers equation has presented a method to reliably estimate this particular solution. The second is encountered when applying the method of characteristics to solve (F.4). This is the approach ultimately taken in Ameli et al. [1], but it is subject to several limitations. The characteristic curves for (F.4) 33

are straightforward to obtain, but at the same time where shocks occur in the Burgers equation formulation, the characteristic curves cross a branch cut in the complex plane. To extend beyond this cut requires analytic continuation of the Stieltjes transform m. Numerical analytic continuation is notoriously ill-posed, with worst-case error from naïve approaches growing doubly exponentially fast in t, or exponentially fast in the original matrix size [48]. In reality, free decompression exhibits worst-case error that grows no more than polynomially in the matrix size [1, Proposition 3], and Ameli et al. [1] presents an approach using polynomial approximation with error rates often growing only logarithmically in the size of the original matrix, with the run-time complexity being independent of t. At these rates, free decompression becomes impressively practical. Unfortunately, the approach is mathematically brittle, highly sensitive to the accuracy of the curve-fit used to obtain the initial data, and imposes assumptions on the form of the spectral density that rarely hold outside of the simplest synthetic cases. Our primary contribution is to identify an approach that maintains this high level of performance under less restrictive assumptions, rendering free decompression practical for more realistic examples.

F.4

Formal Free Decompression along Spectral Curves

Now, recall that the functional relation for the R-transform in terms of the Stieltjes transform is given by Couillet & Liao [13, Definition 5] as 1 R(−m(z)) = z + . m(z) According to the Nica–Speicher Theorem, the operation of free decompression (inverting the operation of free compression) coincides with the map ν 7→ ντ where the R-transforms R and Rτ of ν and ντ , respectively, satisfies Rτ (z) = R(τ z), τ > 1. If m is defined to satisfy the formal algebraic relation P (z, m) = 0 for some polynomial P , then, formally, we can track the behavior of the Stieltjes transform of ντ under free decompression. In doing so, we can track solutions to free decompression while avoiding any potential issues with m and/or R being multivalued. This is captured in Theorem F.4. Theorem F.4. For a polynomial P (z, m), define the following spectral curves: Cm = {(z, m) ∈ C2 | P (z, m) = 0},

(graph of m)

CR = {(−m, z + 1/m) | (z, m) ∈ Cm },

(R transform of m)

CRτ = {(w, r) ∈ C2 | (τ w, r) ∈ CR },

(R transform of mτ )

Cmτ = {(z, ω) ∈ C | ω ̸= 0 and (−ω, z + 1/ω) ∈ CRτ }. Then Cmτ ⊆ CPτ :=

(z, ω) ∈ C2 P

(graph of mτ )

τ −1 , τω τω

=0 .

Proof. We begin by unpacking Cmτ . Note that (z, ω) ∈ Cmτ if and only if (excluding the ω ̸= 0 case), (−ω, z + 1/ω) ∈ CRτ . This implies that (−τ ω, z + 1/ω) ∈ CR , which in turn implies that there exists (z0 , m) ∈ Cm such that 1 1 −m = −τ ω and z0 + =z+ . m ω Consequently, τ −1 m = τ ω and z0 = z + . τω Since (z0 , m) ∈ Cm implies that P (z0 , m) = 0, this implies that there exists z0 such that P (z0 , τ ω) = 0

and

z0 = z +

τ −1 . τω

Consequently, any (z, ω) ∈ Cmτ must also satisfy (z, ω) ∈ CPτ . Theorem F.4 immediately implies Theorem 1, as we identify mτ to be the free decompression of m by the factor τ . 34

F.5

Polynomial Evolution under Free Decompression

In this section, we use Theorem F.4 to derive a bivariate polynomial Qτ such that the Stieltjes transform mτ under free decompression by a decompression ratio τ satisfies Qτ (z, mτ (z)) = 0 for almost all z. Since this polynomial can be readily obtained from P and the roots obtained using standard methods, it provides a fast method of constructing candidate solutions for mτ . The polynomial is constructed in Corollary F.5. Corollary F.5. Suppose that m0 : C+ → C+ satisfies an algebraic relation of the form P (z, m0 (z)) = 0 for z ∈ C+ where dz s X X P (z, m) = cjk z j mk . j=0 k=0

Then for Cmτ defined in Theorem F.4, n o Cmτ ⊆ (z, m) ∈ C2 Qτ (z, m) = 0 , where Qτ (z, m) is the polynomial given by Qτ (z, m) :=

dz dz +s X X (τ ) k l ckl z m , k=0 l=0

and ) c(τ pq :=

dz X

cj,q−dz +j−p τ q−dz +j−p

j=p

j (1 − τ −1 )j−p , p

where cij = 0 if j ∈ / {0, . . . , s}. Proof. From Theorem F.4, we have the inclusion τ −1 2 Cmτ ⊆ (z, m) ∈ C P z + , τm = 0 . τm In particular, for (z, m) ∈ Cmτ , m ̸= 0 and mdz P

τ −1 , τ m = 0. τm

We will show that the left hand side of this expression is the polynomial Qτ (z, m). First note that P

X j dz s X τ −1 τ −1 cjk z + , τm = (τ m)k . τm τ m j=0 k=0

Expanding the binomial gives j X i X j j τ −1 j j−i τ − 1 j j−i z+ = z = z (τ − 1)i τ −i m−i . τm i τ m i i=0 i=0 Therefore, mdz P

X dz j s X X τ −1 j , τm = cjk (τ − 1)i τ k−i z j−i mk−i+dz . τm i j=0 i=0 k=0

Now we introduce new summation indices in place of i, k by letting p = j − i and q = k − i + dz . Then i = j − p and k =q + i − dz = q − dz + j − p and p ∈ {0, . . . , dz }, q ∈ {0, . . . , dz + s}. Using these indices, j since j ≥ p and ji = j−i , mdz P

X dz dz +s dz XX τ −1 j , τm = cj,q−dz +j−p (τ − 1)j−p τ q−dz z p mq , τm p p=0 q=0 j=p 35

where we have adopted the convention that cj,k = 0 if k ∈ / {0, . . . , s}. Since τ q−dz (τ − 1)j−p = τ q−dz +j−p (1 − −1 j−p τ ) . Hence X dz dz +s X (τ ) p q τ −1 mdz P z + , τm = cpq z m , τm p=0 q=0 as required. As stated, Qτ need not have the same degree in m as the initial polynomial P . However, we recall that polynomial relations of this form are not unique, and they can be turned into a unique minimal form by P ) dividing by common factors of z and m, and normalizing so that p c(τ pq = 1. We refer to the endpoint of this reduction pipeline as the polynomial Pτ . For common examples (e.g., Marchenko-Pastur or semicircle laws), Pτ often has the same degree as P , due to Lemma F.6 below. Lemma F.6. If the index set A for P is upper-triangular, then degm (Pτ ) ≤ degm (P ). Proof. Note that it will suffice to show that P (z + (1 − τ −1 )m−1 , τ m) is a bivariate polynomial in z and m, as this will imply that mdz can be factorized from Qτ . If A is upper-triangular, then every surviving monomial z j mk in the expansion for P satisfies k ≥ j. But this implies that P

X j dz s−j X τ −1 τ −1 z+ , τm = cj,j+k z + (τ m)j+k τm τ m j=0 k=0

dz s−j X X

cj,j+k (τ zm + τ − 1)j (τ m)k ,

j=0 k=0

which has only positive powers of z and m. For a bivariate polynomial of the form P (z, m) = q(z)m2 −p(z)m+1, the condition in Lemma F.6 suggests that degz (q) ≤ 2 and degz (p) ≤ 1 implies that Pτ does not have a higher degree than P . Many ensemble models exhibit bivariate polynomial relations of this form, including [54] semicircle law (free Gaussian); [26] (free Poisson); [22] and [27] (free binomial); [51] (free Jacobi); and the general free-Meixner family [3, 41] (free analogs of the classical [28] laws). The Q and P for these ensemble models can be found in [1, Table B.3]. However, P need not be upper-triangular in general; for example, the deformed Wigner ensemble [35, Sections 2.2 and 18.3] is not upper-triangular. This increase in degree exacerbates the most significant difficulty with performing free decompression using Corollary F.5: identifying the correct root of Pτ is extraordinarily challenging. For this reason, we will require a more refined approach by tracing along spectral curves to identify the correct root.

F.6

Free Decompression on Spectral Curves

The previous subsection shows that, once the initial algebraic curve P (z, m) = 0 is given, free decompression produces an evolved algebraic relation whose roots provide all candidate values of the transformed Stieltjes function at time τ . However, this polynomial evolution alone does not identify the physical branch. This is illustrated in Figure F.1(a), where the roots of the evolved polynomial form a cloud of candidate densities, only one of which agrees with the empirical spectrum. This root-selection problem is fundamentally different from the one in Appendix E.1. There, for the initial curve, one may choose an anchor point directly in the z-plane and continue along that curve. After free decompression, by contrast, the evolved polynomial at a fixed τ > 1 provides no new anchor identifying the physical sheet. Instead, the physical branch must be continued from the known physical branch at τ = 1 along the free-decompression flow itself. We now formulate this continuation geometrically. We write mτ (z) for the Stieltjes transform after decompression ratio τ , equivalently at time t = log τ . In particular, m0 denotes the initial Stieltjes transform

1 =mi (λ + δi) π

104

(a) All Polynomial Roots

(b) Density from Physical Root

10−1 10−6

10−11 10

Physical density Empirical histogram Poisson kernel δ floor

Roots Empirical Histogram

Physical density Empirical histogram

−16

10−5

10−3

10−1

101

103

10−5

10−3

10−1

101

103

10−5

10−3

10−1

101

103

Figure F.1: Free decompression on the fitted algebraic spectral curve for the diffusion-model example at the final size n = 64,000, evolved from the initial n = 4000 submatrix. Panel (a) shows all candidate densities π −1 ℑmi (λ + iδ) with δ = 10−6 , obtained from the seven roots mi of the evolved polynomial relation with degm (P ) = 7, at the final decompression time. These unordered candidates form a cloud of physical and nonphysical sheets. Panel (b) shows the result after tracking the physical branch along the free-decompression flow, yielding a coherent density curve that matches the empirical histogram on the bulks, but still exhibits the finite-δ Poisson floor in the spectral gaps. Panel (c) shows the final density after δ-extrapolation to δ = 0 using four δ-levels and a quadratic Chebyshev fit, which removes the artificial gap baseline while preserving the bulk densities. at t = 0, or τ = 1, and should not be confused with the indexed algebraic roots mi used earlier. Let the initial fitted spectral curve be n o C = (ζ, y) ∈ C2 | P (ζ, y) = 0 . Here, y = m0 (ζ) denotes a branch of the initial algebraic Stieltjes transform. Under free decompression by ratio τ ≥ 1, the method of characteristics gives the change of variables z=ζ−

τ −1 , y

y , τ

where w = mτ (z) is the decompressed Stieltjes value. Equivalently, ζ =z+

τ −1 , τw

y = τ w.

Substituting this inverse map into P (ζ, y) = 0 gives the evolved algebraic relation τ −1 P z+ , τ w = 0. τw This is precisely the relation obtained in Theorem F.4. After clearing denominators, it gives the evolved polynomial whose roots provide all algebraic candidates for w = mτ (z), but not the physical one. To track the physical root, we work instead in the original curve coordinates (ζ, y) ∈ C. For a fixed target point z, define F1 (ζ, y) := P (ζ, y), τ −1 F2 (ζ, y; τ, z) := ζ − − z. y Thus, a point (ζ, y) ∈ C maps to the target z at decompression ratio τ precisely when F1 (ζ, y) = 0,

F2 (ζ, y; τ, z) = 0.

At τ = 1, this system reduces to ζ = z and y = m0 (z), so the initial physical point is known from the branch-selection procedure in Appendix E.1. The task is then to continue this solution in τ , while keeping z fixed. 37

In practice, evolution along the spectral curve can pass through regions where sheets come close, roots change order, or the geometry becomes stiff. To resolve this, we evolve the above system using a predictor– corrector approach, where the predictor uses the local tangent along the constrained curve. Namely, differentiating F1 = 0 with respect to τ , while keeping z fixed, gives ∂ζ P (ζ, y)ζ̇ + ∂y P (ζ, y)ẏ = 0. Similarly, differentiating F2 = 0, namely ζ− gives

τ −1 − z = 0, y

ζ̇ − y −1 + (τ − 1)y −2 ẏ = 0,

or equivalently,

ζ̇ + (τ − 1)y −2 ẏ = y −1 .

Combining these two relations gives the 2 × 2 tangent system ∂ζ P (ζ, y) ∂y P (ζ, y) 0 ζ̇ = . 1 (τ − 1)y −2 ẏ y −1 This tangent system provides the predictor direction for advancing the physical point from one value of τ to the next. The predicted point is then corrected by Newton iteration on the nonlinear system (F1 , F2 ) = (0, 0). After correction, the decompressed Stieltjes value is recovered as mτ (z) = w =

y . τ

The resulting geometric continuation procedure is summarized in Algorithm F.1. In the implementation, the Newton correction is damped and the τ -steps are adaptively subdivided when the local geometry becomes stiff; these safeguards prevent jumps between nearby sheets without changing the underlying geometric continuation described above. The effect of this procedure is shown in Figure F.1. The evolved polynomial alone produces the unordered candidate cloud in panel (a), while the geometric continuation in τ selects the physical sheet shown in panel (b). This is a stringent test case: the density is log-scale, multi-bulk, and highly unbalanced, with nearby sheets interacting over several orders of magnitude. Nevertheless, the continuation remains on the physical branch and produces a coherent density matching the empirical bulks; the remaining finite-δ floor is then removed by the extrapolation step in panel (c).

Appendix G

Evolution of Spectral Quantities

Free decompression primarily evolves the full density ρ(t, x), or equivalently its Stieltjes transform m(t, z). In many applications, however, one only needs derived spectral quantities such as edge locations, cusps (merging/splitting points), gap widths, atoms’ weights and locations, or spectral moments. This section shows that these quantities can often be evolved directly from the spectral curve, without reconstructing the density pointwise.

G.1

Evolving Spectral Edges

The spectral edges of a density are of particular interest; for instance, [10] shows that the locations of bulk edges play an important role in the dynamics of diffusion models during training. Informally, if the support of the density ρ(t, x) is a union of k disjoint intervals I(t) :=

k [

Ij (t),

Ij (t) := [aj (t), bj (t)],

j=1

Algorithm F.1: Free decompression by geometric continuation on a spectral curve : Polynomial P ∈ R[z, m] such that P (z, m0 (z)) = 0 (see Algorithm D.1); Nz ⊂ C+ ; Query grid {zj }j=1 Decompression ratios 1 = τ0 < τ1 < · · · < τT ; Oracle for the initial physical branch m0 (z) on P (z, m) = 0 (see Algorithm E.1); Newton tolerance ϵ > 0.

Input

Output : Physical Stieltjes values mτℓ (zj ) for all requested τℓ and zj . // Initialize the physical sheet at τ0 = 1

1 for j = 1, . . . , Nz do (0)

← m0 (zj )

(0) ζj ← zj

mτ0 (zj ) ← yj

// Select physical root using Appendix E.1 // At τ = 1, the map is identity in z

(0)

// Continue the physical sheet in decompression ratio τ 5 for ℓ = 1, . . . , T do 6 for j = 1, . . . , Nz do 7

(ℓ−1)

Set (ζ, y) ← (ζj

(ℓ−1)

, yj

// Predict along the tangent of the constrained system 8

9 10 11

Solve

0 ∂ζ P (ζ, y) ∂y P (ζ, y) ζ̇ = −1 . y 1 (τℓ−1 − 1)y −2 ẏ

∆τ ← τℓ − τℓ−1 . ζpred ← ζ + ∆τ ζ̇. ypred ← y + ∆τ ẏ. // Correct by projecting back to the decompression system

Starting from (ζpred , ypred ), apply damped Newton to solve P (ζ, y) = 0,

13 14 15 16 17

ζ − (τℓ − 1)y −1 − zj = 0.

if Newton converges to (ζ∗ , y∗ ) within tolerance ϵ then (ℓ) (ℓ) ζj ← ζ∗ , yj ← y∗ . mτℓ (zj ) ← τℓ−1 y∗ . // Recover evolved Stieltjes value else Subdivide [τℓ−1 , τℓ ] and repeat the predictor–corrector step. // Avoid sheet jumping in stiff regions

18 return {mτℓ (zj )}ℓ,j .

then the edges are the endpoints {aj (t), bj (t)}kj=1 , where ρ(t, x) vanishes outside I(t) and is positive inside the interior of each Ij (t). When the Stieltjes transform m(t, z) is algebraic, the notion of edge admits a precise characterization in terms of the branch points of the spectral curve. Concretely, suppose m is a branch m = m1 (z) of the algebraic relation P (z, m) = 0, and define the spectral curve n o C := (ζ, y) ∈ C2 | P (ζ, y) = 0 . The curve C defines a branched covering of the ζ-plane, with local sheets given by the roots y = yi (ζ) of 39

P (ζ, y) = 0. A branch point is a point (ζ∗ , y∗ ) ∈ C where two (or more) local sheets coalesce; equivalently, P (ζ∗ , y∗ ) = 0,

and

∂y P (ζ∗ , y∗ ) = 0.

Branch points determine how sheets glue across branch cuts. Among branch points, spectral edges are distinguished by the physical sheet. A real branch point x∗ ∈ R is a spectral edge if the associated branch cut connects the physical sheet to a non-physical sheet, so that crossing the cut switches between m1 and a different branch. In contrast, other branch points whose cuts glue two non-physical sheets do not bound the support of the density (even if they are real branch points). In multi-bulk situations, this distinction is essential: not every real branch point corresponds to an endpoint of the support. With this notion in mind, we first characterize the evolution of branch points under free decompression, and then interpret which of the resulting real branch points correspond to spectral edges via their sheet connectivity to the physical branch. Proposition G.1 (Evolving Branch Point under Free Decompression). Let (z∗ , m∗ ) ∈ C2 be a branch point of m(t, z) with m∗ := m(t, z∗ ). Define the pullback (ζ∗ , y∗ ) = ϕ−1 t (z∗ , m∗ ). Then, (ζ∗ , y∗ ) ∈ C, when exists, satisfies P (ζ∗ , y∗ ) = 0,

(G.1a)

y∗2 ∂y P (ζ∗ , y∗ ) − (τ − 1) ∂ζ P (ζ∗ , y∗ ) = 0,

(G.1b)

followed by the pushforward z∗ = ζ∗ − (τ − 1)y∗−1 . In particular, real solutions z∗ ∈ R of (G.1) include the evolving spectral edges of the real density. Proof. The spectral edges of the real density at time t correspond to real branch points of the map (ζ, y) 7→ z = ζ − (τ − 1)y −1 restricted to the spectral curve C. Along the spectral curve C, implicit differentiation of P (ζ, y) = 0 gives dP (ζ, y) =

∂P (ζ, y) ∂P (ζ, y) dζ + dy. ∂ζ ∂y

On the spectral curve, P (ζ, y) = 0, hence dP (ζ, y) = 0. Whenever ∂y P (ζ, y) ̸= 0, it holds ∂ζ P (ζ, y) dy =− . dζ ∂y P (ζ, y) On the other hand, differentiating the projection π ◦ ϕt (ζ, y) = z(ζ, y) = ζ − (τ − 1)y −1 along C yields ∂ζ P (ζ, y) dz dy = 1 + (τ − 1)y −2 = 1 − (τ − 1)y −2 . dζ dζ ∂y P (ζ, y) The critical points of the projection z = π ◦ ϕt are characterized by dz/dζ = 0 together with P (ζ, y) = 0, which is equivalent to (G.1). Since spectral edges are a subset of branch points, Proposition G.1 applies in particular to the evolution of the edges {aj (t), bj (t)}kj=1 . For each fixed t ≥ 0, the algebraic system (G.1) may admit a set of solutions n op B(t) := z∗(j) (t)

j=1

where typically p ≥ k. The real subset B(t) ∩ R provides candidates for spectral edges; one then selects those candidates whose cuts connect the physical sheet to a non-physical sheet. Operationally, this can be checked by whether the physical branch exhibits a nontrivial jump across the real axis at that point (equivalently, by the Stieltjes inversion formula): a candidate x∗ ∈ R corresponds to an edge of the density precisely when the boundary values m1 (t, x∗ + i0+ ) and m1 (t, x∗ + i0− ) switch sheets across the associated cut. Algorithm G.1 summarizes the resulting procedure. 40

Algorithm G.1: Evolution of spectral edges under algebraic free decompression : Polynomial P ∈ R[z, m] with degm (P ) = s such that P (z, m0 (z)) = 0 (see Algorithm D.1); Time grid 0 = t0 < t1 < Sk·0· · < tT ; Initial support I(0) = j=1 [aj (0), bj (0)] (optional, for initialization); Oracle for the physical branch m0 (z) on C \ I (see Algorithm E.1); // Only needed if I(0) is

Input

provided.

Tolerance ϵ > 0. Output : Edge trajectories {aj (tℓ ), bj (tℓ )} stored as T + 1 by 2k0 with NaNs after merging; Active bulk count k(tℓ ) (number of connected components). // Define the edge system from Proposition G.1 t

1 Set τ ← e

and define the functions G1 (ζ, y; t) := P (ζ, y) G1 (ζ, y; t) := y 2 ∂y P (ζ, y) − (τ − 1)∂ζ P (ζ, y).

// Initialization at t0 = 0 2 if initial edges I(0) are provided then (j)

(j)

0 Initialize guesses {(ζ0 , y0 )}j=1 as follows: ▷ Sample 2k0 points ζ near each endpoint

3 4

// Any reasonable edge-local seeding works

(j)

▷ Choose Herglotz-consistent root y0 ← m0 (ζ0 ).

5 6 else

(j)

0 Obtain initial guesses {(ζ0 , y0 )}j=1 by solving (G1 , G2 ) = (0, 0) at t0 using multiple seeds.

// Track edge preimages (ζ, y) by warm-started Newton over time 8 for ℓ = 0, 1, . . . , T do 9 10 11 12 13 14

τ ← etℓ foreach j = 1, · · · , 2k0 not yet discarded do (j) (j) Use previous solution (ζℓ−1 , yℓ−1 ) as initial guess. (j)

(j)

Apply damped Newton to solve (G1 (ζ, y; tℓ ), G2 (ζ, y; tℓ )) = (0, 0) for (ζℓ , yℓ ). (j) (j) zℓ ← ζℓ − (τ − 1)/yℓ (j) (j) (j) if |ℑzℓ | < ϵ then Store xℓ ← ℜzℓ

// Pushforward to the z-plane

else mark failure.

// Handle merged edges of disjoint bulks (inner edges annihilate 16 17

P ← {a0 , b1 , . . . , ak0 , bk0 } foreach adjacent pair (bj , aj+1 ) in P do if bj (tℓ ) > aj+1 (tℓ ) then (bj (tℓ ), aj+1 (tℓ )) ← (NaN, NaN)

k(tℓ ) ← number of remaining connected components.

19 return Spectral edges

n oT k0 {aj (tℓ ), bj (tℓ )}j=1

ℓ=1

G.2

and active bulk count {k(tℓ )}Tℓ=1 ;

Spectral Cusps and The Mechanism of Merging Bulks

In multi-bulk densities (i.e., k > 1), adjacent bulks may merge under free decompression. We describe this merger through the geometry of the spectral curve, emphasizing that branch points persist while their sheet connectivity (monodromy) changes. Let Ij (t) := [aj (t), bj (t)] denote the j-th bulk. Assume that at t = 0 the bulks are disjoint, a1 (0) < b1 (0) < a2 (0) < b2 (0) < · · · < ak (0) < bk (0), so that Ij (0) ∩ Ij+1 (0) = ∅. Suppose that at some time t = t∗ two consecutive bulks collide, meaning 41

bj (t∗ ) = aj+1 (t∗ ). For t > t∗ the union Ij (t) ∪ Ij+1 (t) becomes a single interval, and the two points that collided cease to be spectral edges of the density because the density no longer vanishes between them. Importantly, however, these points do not disappear from the spectral curve: Proposition G.1 continues to evolve them as bona fide branch points of C, with bj (t) > aj+1 (t) at t > t∗ . The correct interpretation is that the merger is a phase transition in sheet gluing. For t < t∗ , the two inner real branch points connect the physical sheet to non-physical sheets, and hence they bound two bulks of the physical density. After t > t∗ , these same branch points typically glue two non-physical sheets to each other; the physical sheet is then connected only through the remaining outer branch points. In this sense, post-merger “ghost edges” are still branch points, but they are no longer edges of the physical density because they no longer involve the physical sheet. At the collision time t∗ , the density develops a cusp (a cubic singularity at a sharp local minimum), associated with the Pearcey universality class [15]. Geometrically, the colliding branch points correspond to a higher-order critical point of the projection of the spectral curve. The following proposition characterizes this cusp point by algebraic conditions on C. Proposition G.2 (Cusp point at Merging Bulks). Suppose two consecutive real edges collide at time t∗ , namely, there exists x∗ ∈ R such that x∗ = bj (t∗ ) = aj+1 (t∗ ). Let m∗ := m(t∗ , x∗ ) and define the pullback (ζ∗ , y∗ ) := ϕ−1 t (x∗ , m∗ ). Then, (ζ∗ , y∗ ) ∈ C is a higher-order critical point of the projection z(ζ, y) = ζ − (τ∗ − 1)y −1 restricted to C, and hence satisfies P (ζ∗ , y∗ ) = 0

(G.2a)

y∗2 ∂y P − (τ∗ − 1) ∂ζ P = 0,

(G.2b)

y∗ (∂ζζ P )(∂y P ) − 2(∂ζy P )(∂ζ P )(∂y P ) + (∂yy P )(∂ζ P )

+ 2(∂ζ P ) (∂y P ) = 0,

(G.2c)

where all derivatives are evaluated at (ζ∗ , y∗ ), and τ∗ := et∗ . Moreover, the cusp location is the pushforward x∗ = ζ∗ − (τ∗ − 1)y∗−1 . Proof. Consider the restriction of the map z(ζ, y) = π ◦ ϕt (ζ, y) = ζ − (τ − 1)y −1 to the spectral curve C = {(ζ, y) ∈ C2 | P (ζ, y) = 0}. Near (ζ∗ , y∗ ) ∈ C, let y = y(ζ) be a local parameterization of C. At cusp point occurs at time t∗ corresponds to a critical point of multiplicity two of the projection ζ 7→ z(ζ, y(ζ)) = z(ζ) along C, i.e., dz(ζ∗ ) d2 z(ζ∗ ) = 0, = 0, dζ dζ 2 together with P (ζ∗ , y∗ ) = 0. The constraint P (ζ∗ , y∗ ) = 0 implies dP/dζ = 0 and d2 P/dζ 2 = 0 at (ζ∗ , y∗ ). Implicit differentiation of P (ζ, y(ζ)) = 0 along C gives dP ∂P ∂P = + y′ , dζ ∂ζ ∂y d2 P

∂2P

dζ

∂ζ

2 =

′ 2 + 2y

∂2P ∂2P ∂P + (y ′ )2 2 + y ′′ , ∂ζ∂y ∂y ∂y

where y ′ := dy/dζ and y ′′ := d2 y/dζ 2 . Solving dP/dζ = 0 and d2 P/dζ 2 = 0 for y ′ and y ′′ yields y′ = −

∂ζ P , ∂y P

(G.3a)

y ′′ = −

∂ζζ P + 2y ′ ∂ζy P + (y ′ )2 ∂yy P , ∂y P

(G.3b)

whenever ∂y P (ζ, y) ̸= 0. 42

Next, differentiating z(ζ) = ζ − (τ − 1)y(ζ)−1 along C yields dz = 1 + (τ − 1)y −2 y ′ , dζ d2 z −3 ′ 2 (y ) + y −2 y ′′ . 2 = (τ − 1) −2y dζ

(G.4a) (G.4b)

Imposing dz/dζ = 0 and using y ′ from (G.3) gives the second equation in (G.2), y 2 ∂y P − (τ − 1)∂ζ P = 0. Next, imposing d2 z/dζ 2 = 0, substitute y ′ and y ′′ from (G.3), and clear denominator by multiplying by y(∂y P )3 . After simplification, this yields the third equation in (G.2), y ∂ζζ P (∂y P )2 − 2(∂ζy P )(∂ζ P )(∂y P ) + (∂yy P )(∂ζ P )2 + 2(∂ζ P )2 (∂y P ) = 0, where all partial derivatives are evaluated at (ζ∗ , y∗ ). Finally, the cusp location in the physical coordinate is the pushforward x∗ = z(ζ∗ , y∗ ) = ζ∗ − (τ∗ − 1)y∗−1 at time t∗ = log(τ∗ ). Practically, (G.2) can be used to solve for the triple (ζ∗ , y∗ , τ∗ ) and hence recover both the merging time t∗ = log(τ∗ ) and the merge location x∗ . This provides a principled and more accurate alternative compared to monitoring edge crossing on a discrete time grid. Once t∗ is identified, one can treat the post-merge configuration as having one fewer bulk, and interpret subsequent real branch points through the updated sheet gluing relative to the physical branch.

G.3

Dynamics of Gap and Support Intervals

Beyond tracking the density itself, it is often useful to characterize whether a bulk interval widens or narrows, and whether a gap between two bulks shrinks toward a merger or expands. Define the bulk widths sj and gap widths gj by sj (t) := |Ij (t)| = bj (t) − aj (t), gj (t) := aj+1 (t) − bj (t). The dynamics of these intervals can be obtained from the relative velocity of edges. We derive an expression for the velocity of an edge, x′∗ (t), in terms of the geometry of the spectral curve, and then specialize at t = 0, where it admits a particularly simple interpretation in terms of the Hilbert transform. Proposition G.3 (Edge Velocity under Free Decompression). Let x∗ (t) ∈ R be a regular spectral edge of the density, and let (ζ∗ (t), y∗ (t)) denote a corresponding preimage on the spectral curve C = {(ζ, y) ∈ C2 | P (ζ, y) = 0} under the free decompression map z = ζ − (τ − 1)y −1 with τ = et . Suppose (ζ∗ , y∗ ) is a simple critical point of the projection (so it is not a cusp point), hence it satisfies P (ζ∗ , y∗ ) = 0,

(G.5a)

F (ζ∗ , y∗ , τ ) = y∗2 ∂y P (ζ∗ , y∗ ) − (τ − 1)∂ζ P (ζ∗ , y∗ ) = 0,

(G.5b)

and the edge location is x∗ (t) = ζ∗ (t) − (τ − 1)y∗ (t)−1 . Then, x∗ (t) is differentiable and its velocity is x′∗ (t) = ζ∗′ (t) − τ y∗ (t)−1 + (τ − 1) y∗′ (t) y∗ (t)−2 , where (ζ∗′ (t), y∗′ (t)) is obtained by solving the linear system ′ ∂ζ P ∂y P ζ∗ 0 , = −τ ∂ζ F ∂y F (ζ ,y ,τ ) y∗′ ∂τ F (ζ ,y ,τ ) ∗

∗

with ∂τ F = −∂ζ P (ζ, y), ∂ζ F = y 2 ∂ζy P − (τ − 1)∂ζζ P, ∂y F = 2y∂y P + y 2 ∂yy P − (τ − 1)∂ζy P. 43

(G.6)

(G.7)

Proof. We differentiate the defining relations (G.5) along the solution curve t 7→ (ζ∗ (t), y∗ (t)). From P (ζ∗ , y∗ ) = 0, we obtain ∂ζ P (ζ∗ , y∗ ) ζ∗′ + ∂y P (ζ∗ , y∗ ) y∗′ = 0. From F (ζ∗ , y∗ , τ ) = 0 and τ ′ (t) = τ (t), we obtain ∂ζ F (ζ∗ , y∗ , τ ) ζ∗′ + ∂y F (ζ∗ , y∗ , τ ) y∗′ + ∂τ F (ζ∗ , y∗ , τ ) τ ′ (t) = 0, which yields the linear system (G.7). Finally, differentiating x∗ (t) = ζ∗ (t) − (τ (t) − 1)y∗ (t)−1 gives (G.6). Corollary G.4 (Initial edge velocity). Let x∗ (0) be a regular spectral edge at t = 0, and let Z 1 ρ(t, y) H(t, x) := H{ρ(t, ·)}(x) = p.v. dy, π R x−y be the Hilbert transform of ρ(t, ·). Then, x′∗ (0) =

1 . πH(0, x∗ )

(G.8)

Proof. At t = 0, we have τ = 1 and x∗ (0) = ζ∗ (0). Also y∗ (0) = m(0, x∗ ). Using the Plemelj formula, m(0, x∗ ) = −πH(0, x∗ ) + iπρ(0, x∗ ), but on the spectral edge ρ(0, x∗ ) = 0, so y∗ (0) = −πH(0, x∗ ). Also, from Proposition G.3, we have x′∗ (0) = ζ∗′ (0) − τ (0) y∗ (0)−1 . Moreover, at t = 0 the criticality condition in (G.5) reduces to ∂y P (ζ∗ (0), y∗ (0)) = 0. If x∗ (0) is a regular edge (not a cusp), then ∂ζ P (ζ∗ (0), y∗ (0)) ̸= 0, and differentiating P (ζ∗ (t), y∗ (t)) = 0 at t = 0 implies ζ∗′ (0) = 0. Hence x′∗ (0) = −y∗ (0)−1 = 1/(πH(0, x∗ )). The dynamics of bulk widths sj (t) and gap widths gj (t) follows immediately. At t = 0, Corollary G.4 yields the explicit identities 1 1 − , πH(0, bj (0)) πH(0, aj (0)) 1 1 gj′ (0) = − . πH(0, aj+1 (0)) πH(0, bj (0)) s′j (0) =

(G.9a) (G.9b)

The sign of the Hilbert transform at inner edges depends on the global arrangement of outer mass across bulks; therefore, inner gaps may either shrink (leading to a merger at a finite cusp time) or expand. In contrast, for the two outermost edges of the total support I(t) = [a1 (t), bk (t)], the signs are fixed: for |x| large, m(0, x) ∼ −1/x, hence H(0, x) ∼ 1/(πx). As a result, H(0, a1 (0)) < 0 and H(0, bk (0)) > 0, implying a′1 (0) < 0 and b′k (0) > 0. Thus, the total support length |I(t)| = bk (t) − a1 (t) is increasing at t = 0.

G.4

Dynamics of Atoms

An atom of µ corresponds to a simple pole of the Stieltjes transform m(z) at z = x◦ , with weight w = − Res m(z) = lim (x◦ − z)m(z). z=x◦

z→x◦

Ps When m is algebraic and satisfies P (z, m) = 0, where P (z, m) := j=0 aj (z)mj , a necessary condition for a pole at a finite real point x◦ is that x◦ is a root of the leading coefficient, i.e., as (x◦ ) = 0. Conversely, a real root of as corresponds to an atom on the physical branch provided that m has a simple pole at x◦ (i.e., the singularity is not canceled and is not of higher order). In that case, if x◦ is a simple root of as , coefficient matching yields a (x◦ ) w = s−1 . a′s (x◦ )

Proposition G.5. Let x◦ ∈ R be an atom of m(z) with initial mass w◦ ∈ [0, 1]. Under the free decompression flow, the atom location remains fixed x◦ (t) = x◦ , while its weight evolves as w(t) = 1 − (1 − w◦ )e−t .

(G.10)

Proof. Let x◦ (t) denote the (a priori time-dependent) location of the atom with weight w(t). The Laurent expansion of m(t, z) near z = x◦ (t) is m(t, z) = − Differentiating gives ∂t m(t, z) = −

w(t) + O(1). z − x◦ (t)

w′ (t) w(t) x′◦ (t) + O(1), + z − x◦ (t) (z − x◦ (t))2

and ∂z m(t, z) =

w(t) (z − x◦ (t))2

+ O(1).

Moreover, the leading term of m−1 near z = x◦ (t) is m(t, z)−1 = −

z − x◦ (t) + O (z − x◦ (t))2 . w(t)

Substituting these expansions into the PDE ∂t m = −m + m−1 ∂z m yields −

w(t) w′ (t) w(t) x′◦ (t) 1 + O(1) = + − + O(1). z − x◦ (t) (z − x◦ (t))2 z − x◦ (t) z − x◦ (t)

(G.11)

The right-hand side has no second-order pole (z − x◦ (t))−2 , hence the coefficient of (z − x◦ (t))−2 on the left-hand side must vanish, implying w(t) x′◦ (t) = 0. If w(t) > 0, then x′◦ (t) = 0, and therefore x◦ (t) = x◦ (0) is constant. Finally, equating the coefficients of the first-order poles in (G.11) gives the ODE w′ (t) = 1 − w(t), which is solved by (G.10) with the initial condition w(0) = w◦ . Remark G.1. The conclusion of Proposition G.5 is local: it does not depend on whether x◦ lies inside, outside, or at the boundary of the absolutely-continuous support of µ0 . Moreover, if µ0 has finitely many (i) na (i) (i) atoms {(x(i) ◦ , w◦ )}i=1 at distinct locations x◦ , then each location x◦ remains fixed and each weight evolves independently according to (G.10). △ Remark G.2. In backward flow with t < 0 (free compression), there is a finite time t◦ := log(1 − w◦ ) < 0 at which the atom weight vanishes, w(t◦ ) = 0. For t < t◦ , the pole (atom) ansatz is no longer valid, reflecting that the atom has disappeared. In particular, the atom dynamics are forward deterministic while the atom is present (w◦ > 0); however, once compression reaches the extinction time w◦ = 0, the compressed measure has lost the memory of the atom location. Consequently, decompression from a post-extinction measure is non-unique unless additional structural information is imposed (e.g., a known atom location such as x◦ = 0 in rank-deficient PSD settings). A concrete example is provided by rank-deficient matrices. Suppose A is an n × n matrix with nullity q < n, so its empirical spectral measure has an atom of weight w◦ = q/n at the origin. Then a principal submatrix A◦ of size n◦ becomes full rank once n◦ ≤ n − q, and hence its empirical spectral measure has no atom at zero. △

G.5

Evolution of Moments

We now turn to computing the moments of the evolving measures ντ under free decompression starting from ν ≡ ν1 . For n ≥ 0, define the time-dependent raw moment Z (τ ) µn := xn dντ (x). 45

These moments can be evolved without reconstructing ντ : the free decompression PDE from Ameli et al. [1, Proposition 1] induces a closed triangular system of ODEs for µn (t) = µexp(t) , allowing each moment to be n computed directly from lower-order moments. Throughout this section, we let m(t, z) denote the Stieltjes transform of νet and assume that m(t, z) ̸= 0 along the free decompression flow. Lemma G.6 (ODEs for Moments). The evolution of the moment µn (t) is governed by the system of the triangularly coupled ODEs, X n 1 d +1 µi (t)µn−i (t) = (n + 1)µn (t), n ∈ N, (G.12) 2 dt i=0 and, in particular, µ0 (t) = 1 and µ1 (t) = µ1 (0) for all t > 0. Proof. Since by the hypothesis m(t, z) ̸= 0, we can rearrange the PDE from Ameli et al. [1, Proposition 1] to 1 ∂ ∂ + 1 m(t, z)2 = m(t, z). (G.13) 2 ∂t ∂z At |z| > max(λ) and under mild regularity conditions, we can expand a convergent Laurent series of m(t, z) by ! Z Z X ∞ ∞ X dσ(λ) λn µn (t) m(t, z) = =− dσ(λ) = − n+1 n+1 . λ−z n=0 z n=0 z The Laurent series of m(t, z)2 becomes 2

m(t, z) =

∞ n X X n=0

Also,

! µi (t)µn−i (t)

i=0

n+2 .

∞ X ∂ µ (t) m(t, z) = (n + 1) nn+2 . ∂z z n=0

(G.14)

Substituting the Laurent series into (G.13) and rearranging for the powers of z −(n+1) , we obtain (G.12). For the first few terms, we find that µ′0 = 1 − µ0 , and since µ0 (0) = 1, µ0 (t) = 1 for all t > 0. Similarly, we find µ′1 (t) = 0 and so µ1 (t) is also constant. The triangular structure of Lemma G.6 suggests solving the moment equations recursively. In particular, the lower-order moments determine the forcing term in the ODE for µn (t), and this forcing is a finite linear combination of powers of et . We therefore seek µn (t) in the form µn (t) =

n−1 X

κn,k ekt ,

or equivalently

k=0

) µ(τ n =

n−1 X

κn,k τ k ,

k=0

where the coefficients κn,k are determined recursively from the initial moments. The following proposition gives the resulting recurrence. Proposition G.7 (Moment Recurrence). For any t ≥ 0 and n ≥ 1, the n-th moment µn (t) satisfies Pn−1 µn (t) = k=0 κn,k ekt where κ1,0 = µ1 (0),

κn,k =

 1 n−1 k XX   2k + 1  κi,j κn−i,k−j ,  n − k − 1

0 ≤ k ≤ n − 2,

    µn (0) −

k = n − 1,

i=1 j=0 n−2 X

κn,j

j=0

and where κn,k = 0 if k ≥ n. 46

(G.15)

Proof. Starting from (G.12),

Substituting µn (t) =

1 d +1 2 dt

Pn−1

k=0 κn,k e

"n−1 X

# µi (t)µn−i (t) + 2µn (t) = (n + 1)µn (t).

i=1

gives

(n + 1)µn (t) = =

1 d +1 2 dt

"n−1 X

µi (t)µn−i (t) + 2

i=1

n−1 X

# kt

κnk e

k=0

n−1 n−1 X X 1 d +1 µi (t)µn−i (t) + κnk (k + 2) ekt , 2 dt i=1 k=0

and hence

n−1 n−2 X X 1 d +1 µi (t)µn−i (t) = (n − k − 1) κnk ekt 2 dt i=1

(G.16)

k=0

Rearranging terms gives n−1 X

µi (t)µn−i (t) =

i=1

n−1 X i=1

   i−1 n−i−1 X X  κij ejt   κn−i,j ejt  j=0

n−1 k X n−2 XX

j=0

κij κn−i,k−j ekt .

(G.17)

i=1 k=0 j=0

Substituting (G.17) into (G.16) gives n−2 X k=0

n−1 n−2 k X XX 1 (n − k − 1) κnk ekt , κij κn−i,k−j ekt = k+1 2 i=1 j=0 k=0

and (G.15) follows by equating terms for 0 ≤ k ≤ n − 2. Only the expression for κn,n−1 remains, which Pn−1 follows from j=0 κn,j = µn (0).

Appendix H

Finite-Size Correction to the Spectral Curve

It is critical to note that the discussion so far has focused on cases where the spectrum of a finite size matrix is approximated by that of a matrix with “equivalent” spectral distribution, but is of infinite size, and therefore has a continuous spectrum. This is typically fine for most applications, but it can fail in cases where a positive-definite matrix is considered with a nonzero minimal eigenvalue, yet the limiting spectrum possesses a left edge at zero. In these cases, especially to estimate the condition number of matrix, accounting for finite size effects becomes essential. This is intractable in the general cases considered here, so we shall instead develop an effective heuristic based on a special class of random matrix distributions. We note that the discussion below is strictly informal, as a robust proof is well beyond the scope of this paper. Let A be an n × n real symmetric matrix and assume that the probability measure over the elements of A satisfies n n Y dAij , p(A) ∝ exp − trV (A) 2 i,j=1 where V is an arbitrary polynomial. Our objective is to estimate m̄n (z) = Emn (z) where mn (z) = n1 trG(z) and G(z) = (A − zI)−1 is the resolvent for z ∈ C \ R. Consider Stein’s Lemma: for any matrix-valued smooth function f = (fij ) that grow subexponentially, Z X n ∂ (fij (A)e− 2 trV (A) )dA = 0, ∂Aij i≤j

and so

    X X ∂fij (A) ∂ n  = nE  fij (A) trV (A) = E[tr(f (A)V ′ (A))]. E ∂Aij 2 ∂Aij 2 i≤j

i≤j

Here, we want to take f (A) = G(z). Note that since V is a polynomial, Q(x, z) = (V ′ (z) − V ′ (x))/(z − x) is a bivariate polynomial. Then G(z)V ′ (A) = (A − zI)−1 V ′ (A) = V ′ (z)(A − zI)−1 + Q(A, z) = V ′ (z)G(z) + Q(A, z). Substituting this into the right-hand side of Stein’s Lemma yields n n Etr(G(z)V ′ (A)) = (V ′ (z)EtrG(z) + EtrQ(A, z)). 2 2 By letting mn (z) = n1 trG(z) and defining the polynomial expectation Qn (z) = E[ n1 trQ(A, z)], the right-hand side simplifies: n n2 ′ Etr(G(z)V ′ (A)) = (V (z)m̄n (z) + Qn (z)). 2 2 To evaluate the left-hand side, we apply the matrix chain rule dG = −G(dA)G to find the derivative ( ∂Gij −(Gii Gjj + G2ij ) for i ̸= j, = ∂Aij −G2ii for i = j. Consequently,

1 = − ((trG(z))2 + tr(G(z))2 ). ∂Aij 2

X ∂Gij i≤j

d Since tr(G(z) ) = dz trG(z) = nm′n (z), 2

X ∂Gij i≤j

∂Aij

=−

n2 n E[mn (z)2 ] − E[m′n (z)]. 2 2

Therefore, E[mn (z)2 ] + V ′ (z)m̄n (z) + Qn (z) +

1 E[m′n (z)] = 0. n

Here, we assert that mn (z) ≈ m(z) + O(1/n), and so Var(mn (z)) = O(n−2 ). To prove this rigorously is quite challenging, and we refer the reader to [18, Section 8]. Consequently, keeping only first-order terms, m̄n (z)2 + V ′ (z)m̄n (z) + Qn (z) +

1 ′ m̄ (z) ≈ 0. n n

(H.1)

Remark H.1. Taking n → ∞ in (H.1) reveals that this large class of random matrix ensembles always exhibits a polynomial relation in the limiting Stieltjes transform m and its argument z that is quadratic in m. △ Observing that the first three terms comprise a bivariate polynomial in z and mn (z), we equate this with the second-order expansion about any point of interest m∗ : 2 P2 (z, m) = P (z, m∗ ) + ∂m P (z, m∗ )(m − m∗ ) + 21 ∂m P (z, m∗ )(m − m∗ )2 ,

and arrive at our correction ansatz: P2 (z, mn (z)) +

1 ′ m (z) = 0 n n

Now we consider how this affects the location of branch points. 48

(H.2)

Proposition H.1. Let (z∗ , m∗ ) ∈ C2 be a branch point of m(t, z) with m∗ := m(t, z∗ ). Define the pullback (ζ∗ , y∗ ) = ϕ−1 t (z∗ , m∗ ). Then, evaluated along the integral curve y(ζ) satisfying the initial differential equation (H.2), the point (ζ∗ , y∗ ) satisfies y∗2 − n(τ − 1)P2 (ζ∗ , y∗ ) = 0, followed by the pushforward z∗ = ζ∗ − (τ − 1)y∗−1 . Proof. From the differential equation P2 (ζ, y(ζ)) + n1 y ′ (ζ) = 0, we immediately have that dy = −nP2 (ζ, y). dζ Substituting this into the criticality condition dz dζ = 0 gives dz = 1 + (τ − 1)y −2 (−nP2 (ζ, y)) = 0, dζ from whence the result follows. A few remarks about Proposition H.1: assuming that y∗ remains finite as n → ∞ and P ≈ P2 , the condition y∗2 − n(τ − 1)P2 (ζ∗ , y∗ ) = 0 implies that lim P2 (ζ∗ , y∗ ) = 0.

n→∞

1 Furthermore, since P2 (ζ, y) = − n1 y ′ , we find that y∗2 + (τ − 1)y∗′ = 0, and so y∗′ = − τ −1 y∗2 . Applying this to the derivative of (H.2) gives ! y∗2 1 1 ′′ ′ + y∗′′ = 0. 0 = ∂ζ P2 + (∂y P2 )y + y = ∂ζ P2 + (∂y P2 ) − n τ −1 n

Therefore, taking n → ∞, we find that y∗2 ∂y P2 − (τ − 1)∂ζ P = 0. Consequently, we recover both of the original branch point conditions in the limit as n → ∞. Most significantly, however, is that we can show that the finite-size correction implies a Tracy–Widom-type scaling of O(n−2/3 ) in the left-edge. Performing a complete analysis in this respect is well beyond the scope of this paper. Instead, we can shortcut the proof by assuming a particular form of mn (zn ) where each zn is the branch point of the finite-size correction mn . In particular, if (z∗ , m∗ ) is the limiting branch point as n → ∞, we may assert that mn (zn ) ≈ m∗ + bn u(a−1 n (zn − z∗ )), for sequences an , bn and a function u. Given that mn has n poles located at each of the eigenvalue locations, and we should expect that an is the appropriate rescaling of the spectral gap as n → ∞, u should contain at least two poles: one for the eigenvalue at the branch point (e.g., the left edge), and another for the next closest eigenvalue. This assumption is sufficient to uniquely determine the limiting behavior, and we will in fact find that u has infinitely many poles to capture all of the eigenvalues. Proposition H.2. Assume that for each n, the finite-order resolvent mn satisfies (H.2) and that (z∗ , m∗ ) is a branch point of P in the sense that P (z∗ , m∗ ) = 0 and ∂m P (z∗ , m∗ ) = 0. Assume that ∂z P (z∗ , m∗ ) ̸= 0 and 2 ∂m P (z∗ , m∗ ) ̸= 0. Consider a sequence (zn , mn (zn )) defined in terms of coefficients an , bn > 0, a parameter ξ, and a function un (ξ) by zn = z∗ + an ξ, mn (zn ) = m∗ + bn un (ξ). and where an , bn → 0 as n → ∞. Assume that un (ξ) converges to a limit u(ξ) pointwise that is meromorphic with at least two poles. Then an ≍ n−2/3 , bn ≍ n−1/3 . 49

Proof. Let xn = an ξ and yn = bn un (ξ) denote the local scaling variables such that zn = z∗ + xn and mn (zn ) = m∗ + yn . Expanding P around z∗ gives P (zn , m∗ ) = P (z∗ , m∗ ) + ∂z P (z∗ , m∗ )xn + O(x2n ),

(H.3a)

∂m P (zn , m∗ ) = ∂m P (z∗ , m∗ ) + O(xn ),

(H.3b)

2 2 ∂m P (zn , m∗ ) = ∂m P (z∗ , m∗ ) + O(xn ),

(H.3c)

Since (z∗ , m∗ ) is a branch point, P (z∗ , m∗ ) = ∂m P (z∗ , m∗ ) = 0. Now define the non-zero constants A = 1 2 2 ∂m P (z∗ , m∗ ) and B = ∂z P (z∗ , m∗ ). Then substituting (H.3) into the definition of P2 gives P2 (zn , mn (zn )) = Ayn2 + Bxn + En , where En = o(xn ) + O(xn yn ) + o(yn2 ). But now m′n (zn ) =

d(m∗ + bn un (ξ)) b = n u′n (ξ), d(z∗ + an ξ) an

and so (H.2) implies that Ab2n un (ξ)2 + Ban ξ +

bn ′ u (ξ) + En = 0. nan n

(H.4)

Now let θn = max{|A|b2n , |B|an , |bn /(nan )|}. Dividing (H.4) by θn gives αn un (ξ)2 + βn ξ + γn u′n (ξ) +

En = 0, θn

where αn = Ab2n /θn , βn = Ban /θn , and γn = bn /(nan θn ). By construction, αn , βn , γn ∈ [−1, 1], and for any n, at least one of these coefficients has an absolute value of 1. By the Bolzano–Weierstrass Theorem, passing to a subsequence if necessary, these coefficients converge to finite limits α, β, γ, respectively, with at least one limit being strictly non-zero. Since En is strictly bounded by terms that vanish faster than θn , the remainder vanishes in the limit. Therefore, taking the limit as n → ∞, αu(ξ)2 + βξ + γu′ (ξ) = 0. (H.5) q β If γ = 0, then αu2 + βξ = 0. This implies that u(ξ) = ± − α ξ, which has a branch cut and is without poles, β 2 contradicting the hypotheses2 . If α = 0 then βξ + γu′ (ξ) = 0. Integrating this reveals that u(ξ) = − 2γ ξ +C

for some constant C which has no poles. If β = 0, the equation becomes αu2 + γu′ = 0. Solving this γ differential equation yields u(ξ) = α(ξ−C) for some constant C > 0. This is meromorphic, but possesses only one pole3 . Therefore, α, β, γ must all be non-zero. For this to occur, it must be the case that

Now, since b2n ≍ an , this implies that b2n ≍

b2n ≍ an ≍

bn . nan

bn 2 nbn

and so b3n ≍ n−1 . Immediately this implies that

= nb1

bn ≍ n−1/3 and an ≍ n−2/3 . It remains only to identify that the limiting u satisfies the condition of being meromorphic with at least two poles. Note that the limiting equation (H.5) is a Riccati equation, and so by applying the transformation γ v ′ (ξ) , u(ξ) = · α v(ξ) and then the linear rescaling w = −(αβγ −2 )1/3 ξ, the differential equation of v(ψ) becomes the Airy equation v ′′ (w) − wv(w) = 0. The solutions for v are linear combinations of the Airy functions Ai and Bi in w [32, 2

This represents the regime where the finite-order correction is removed entirely, and so only describes the asymptotic macroscopic square-root behavior near the edge of the spectral density. 3 This corresponds to the regime where the spectral gap has been stretched to the point where only one pole at the location of the left edge remains.

Left Edge of Spectral Density

26 K 25 K

Size n

24 K 23 K 22 K 21 K 20 K

Tracy-Widom regime λmin − λ∞ ∼ n−2/3 Exponential regime λmin ∼ e−αn Empirical left edge

10−4

10−3

λmin

Figure H.1: Simulated (black) smallest eigenvalue λmin of the diffusion model matrix example, with Tracy– Widom (blue) predicted decay at rate ≍ n−2/3 and exponential (red) predicted decay at rate ≍ e−αn for fitted α. Section 9.2(i)]. These are entire functions with infinitely many roots on the real axis where the derivative does not vanish [32, Section 9.9]. Therefore u(ξ) is meromorphic with infinitely many poles on the real axis, and satisfies the required conditions. The finite-size correction becomes particularly relevant in our diffusion model example, as free decompression predicts that the left-edge should go to zero, even though the matrix in question is strictly positive-definite. Fortunately, the decay rate in the intermediate regime appears to be well-captured by the Tracy–Widom-type scaling behavior predicted in Proposition H.2; see Figure H.1. Unfortunately, as n grows larger, the left-edge begins to obey an exponential regime instead. This regime is not captured by the finite-size correction. Instead, the behavior is better predicted due to large deviation theory, as it is known that the left tail for very large n exhibits exponential decay [40]. We do not expect that simple corrections to the free decompression operation will adequately capture this behavior at these size scales.

Appendix I

Random Matrix Models

This section describes the random matrix models and empirical datasets used in Section 4. We begin with the analytically tractable compound free Poisson and free Lévy laws in Appendix I.1 and Appendix I.2, whose Stieltjes transforms satisfy explicit algebraic relations and therefore provide controlled benchmarks for free decompression. We then discuss the Pennington–Bahri Hessian model in Appendix I.3, which is a neural-network-inspired special case of the algebraic free Lévy family. Finally, we describe the empirical Hessian and diffusion-model matrices in Appendix I.4 and Appendix I.5, which serve as large-scale data-driven examples.

I.1

Compound Free Poisson Law

The compound free Poisson law [31, p. 218] is a convenient multi-bulk testbed whose Stieltjes transform satisfies an explicit algebraic relation. It may be viewed as a natural generalization of the Marchenko–Pastur law, with the usual Marchenko–Pastur case recovered as the simplest one-parameter member of the family. This distribution admits a natural random-matrix realization [44, Section 22.5.4]. Let Z ∈ Rp×n , with p/n → λ ∈ (0, ∞), have iid centered entries with variance 1/n. Suppose Σ is deterministic positive semidefinite with empirical spectral distribution converging weakly to a compactly supported probability measure ν. Then the companion sample covariance matrix S :=

1 ⊺ X X, n

X := Σ 2 Z, 51

(I.1)

has a limiting eigenvalue distribution given by the compound free Poisson law with rate λ and jump distribution ν. When λ < 1, the companion matrix S has (1 − λ)n zero eigenvalues asymptotically, and the limiting law therefore contains an atom at the origin of mass max{1 − λ, 0}. The limiting eigenvalue distribution above can be characterized through its R-transform described below. Let ν denote the jump distribution of the compound free Poisson law. For the closed-form algebraic examples used here, we take ν to be a finite-atom probability measure, ν(dx) :=

r X

wi δti (dx),

(I.2)

i=1

where δti denotes a unit point mass at x = ti , the atoms are distinct, wi > 0, and R-transform is then Z r X x ti R(w) = λ ν(dx) = λ , wi 1 − ti w R 1 − xw i=1

i=1 wi

= 1. The (I.3)

where λ > 0 denotes the intensity, or rate. Let m(z) be the Stieltjes transform of the law. The standard relation between the Stieltjes transform and the R-transform reads; see, e.g., , [31, Theorem 12.7], z = R(w) +

1 , w

w = −m(z),

(I.4)

where the physical branch of m is selected by the Herglotz condition ℑm(z) > 0 for z ∈ C+ . Q r When ν is finitely supported as in (I.2), multiplying (I.4) by the common denominator m i=1 (1 − ti m) yields an explicit algebraic constraint P (z, m) = 0, (I.5) where P has degree r + 1 in m. In the special case r = 1, this relation reduces to the familiar quadratic equation for the Marchenko–Pastur Stieltjes transform. Example I.1. For the two-atom case r = 2 used in our experiments, P is cubic: P (z, m) := a3 (z)m3 + a2 (z)m2 + a1 (z)m + a0 (z) = 0, with coefficients

(I.6)

a3 (z) = z t1 t2 , a2 (z) = z(t1 + t2 ) + (1 − λ)t1 t2 , a1 (z) = z + (t1 + t2 ) − λ(w1 t1 + w2 t2 ),

(I.7)

a0 (z) = 1, where w2 := 1 − w1 . Hence, m(z) is a three-sheet algebraic function, providing a controlled benchmark for continuation and evolution in multi-bulk settings. △ For the numerical example using the compound free Poisson law, we use λ = 0.1, (t1 , t2 ) = (2, 5 12 ), and (w1 , w2 ) = (0.75, 0.25). This produces a cubic algebraic Stieltjes transform, so in the numerical fit we use (degm , degz ) = (3, 1) and impose only the zeroth-moment normalization constraint. The cubic case also provides a useful visualization of the sheet structure underlying algebraic Stieltjes transforms. Figure I.1 shows the three branches of m(z) separately, together with the half-plane gluings that analytically continue the physical branch across its cuts. Figure I.2 shows the same continuation geometry as a three-sheet surface. These figures illustrate the central advantage of working with the full algebraic relation P (z, m) = 0: the physical Stieltjes transform is not an isolated single-valued function, but one sheet of a larger spectral curve whose other sheets encode its analytic continuations.

I.2

Free Lévy Law

The compound free Poisson law can be generalized by adding a drift and a semicircular, or free Gaussian, component. This yields the free Lévy law νFL = νδa ⊞ νSC 2 ⊞ νCFP(λ,ν) , σ

(a) m1 on C

(b) m2 on C

(ab) m1 on C+ glued to m2 on C−

(ac) m1 on C+ glued to m3 on C−

(ba) m1 on C− glued to m2 on C+

(ca) m1 on C− glued to m3 on C+

Im(z)

−2 0

Re(z)

Im(z)

−2

Im(z)

−2 2

Re(z)

Figure I.1: Branches and analytic continuations of the cubic compound free Poisson Stieltjes transform in (I.6). The first row shows the three algebraic branches m1 , m2 , m3 on the z-plane, with branch cuts marked on the real axis. The physical branch is m1 . The lower panels show the four possible half-plane gluings involving the physical sheet: for example, (ab) denotes m1 on C+ glued to m2 on C− , while (ba) denotes the opposite gluing. Across a shared cut, the individual branches are discontinuous, but the glued half-planes form a continuous analytic continuation. Colors encode the complex value of the branch using the domain-coloring convention of Appendix A. where δa denotes the deterministic shift X 7→ X + a, SCσ2 is Wigner’s semicircle law with variance σ 2 , and CFP(λ, ν) is the compound free Poisson law with rate λ and jump distribution ν, as in Appendix I.1. In full generality, the free Lévy–Khintchine representation characterizes freely infinitely divisible distributions: when ν is a Lévy measure, the corresponding law is freely infinitely divisible, and conversely every freely infinitely divisible law admits such a representation [31, Theorem 13.8]; see also Benaych-Georges [4] for random-matrix constructions associated with freely infinitely divisible laws. In this work, we restrict to the finite-atom case (I.2), which preserves an algebraic Stieltjes transform while still producing a broad family of multi-bulk distributions. The free Lévy law also admits a natural random-matrix realization by combining the drift, Gaussian, and jump components; see Benaych-Georges [4], Burda et al. [11] for related random-matrix realizations of free Lévy and freely infinitely divisible laws. Let W be an n × n Wigner matrix with independent centered entries of variance 1/n, up to symmetry, and let S be the companion sample covariance matrix in (I.1), whose limiting law is CFP(λ, ν). Assuming W and S are independent, the additive model A := aI + σW + S has an empirical spectral distribution converging to νFL under standard assumptions. 53

Figure I.2: Three-sheet visualization of the cubic compound free Poisson spectral curve from Figure I.1. The middle sheet represents the physical branch m1 , while the upper and lower sheets represent the non-physical branches m2 and m3 . The sheets are vertically separated for visualization and are brought together along shared branch cuts to illustrate analytic continuation between branches. This geometric view explains why representing the Stieltjes transform as one sheet of a larger algebraic curve is useful: crossing a branch cut does not terminate the function, but moves the continuation to another sheet. The R-transform of νFL is the sum of the R-transforms of its free additive components: RFL (w) = Rδa (w) + RSC 2 (w) + RCFP(λ,ν) (w). σ

Here, Rδa (w) = a, RSC 2 (w) = σ 2 w, and, for the finite-atom jump distribution (I.2), RCFP(λ,ν) is given by σ (I.3). Hence, r X ti RFL (w) = a + σ 2 w + λ wi , (I.8) 1 − ti w i=1 which is a rational function. Using the relation between the R-transform and the Stieltjes transform in (I.4), with w = −m(z), we obtain the self-consistent equation −

r X 1 ti = z − a + σ 2 m(z) − λ wi . m(z) 1 + t i m(z) i=1

Qr Multiplying by the common denominator m i=1 (1+ti m) yields the implicit algebraic constraint P (z, m) = 0, where r r Y X Y P (z, m) := σ 2 m2 + (z − a)m + 1 (1 + ti m) − λm wi ti (1 + tj m). (I.9) i=1

i=1

j̸=i

For σ = 0, the polynomial (I.9) has degree r + 1 in m and reduces to the compound free Poisson relation, up to the drift shift a. When σ > 0, the additional semicircular term increases the degree to r + 2. Example I.2. For the two-atom case r = 2 and σ > 0, (I.9) becomes a quartic polynomial in m, P (z, m) = a4 (z)m4 + a3 (z)m3 + a2 (z)m2 + a1 (z)m + a0 (z),

(I.10)

(a) m1 on C

(b) m2 on C

(d) m4 on C

(ab) m1 on C+ glued to m2 on C−

(ac) m1 on C+ glued to m3 on C−

(ad) m1 on C+ glued to m4 on C−

(ba) m1 on C− glued to m2 on C+

(ca) m1 on C− glued to m3 on C+

(da) m1 on C− glued to m4 on C+

Im(z)

2 0 −2 −2

Re(z)

Im(z)

2 0 −2

Im(z)

2 0 −2 −2

Re(z)

−2

Re(z)

−2

Re(z)

Figure I.3: Branches and analytic continuations of the quartic free Lévy Stieltjes transform. The first row shows the four algebraic branches of P (z, m) = 0 on the z-plane, with branch cuts marked on the real axis. The lower panels show half-plane gluings between the physical sheet and the non-physical sheets; for example, (ad) denotes the upper half-plane of the physical sheet glued to the lower half-plane of the branch shown in panel (d), while (da) denotes the opposite gluing. Across a shared cut, the separate branches are discontinuous, but the glued half-planes form a continuous analytic continuation. Colors encode the complex value of the branch using the domain-coloring convention of Appendix A. with coefficients

a4 (z) = σ 2 t1 t2 , a3 (z) = (z − a)t1 t2 + σ 2 (t1 + t2 ), a2 (z) = (z − a)(t1 + t2 ) + (1 − λ)t1 t2 + σ 2 ,

(I.11)

a1 (z) = (z − a) + (t1 + t2 ) − λ(w1 t1 + w2 t2 ), a0 (z) = 1. In particular, for σ = 0, the leading coefficient a4 vanishes and (I.10) reduces to the cubic compound free Poisson law with drift. △ For the numerical example below, we use λ = 0.1, (t1 , t2 ) = (2, 5 12 ), (w1 , w2 ) = (0.75, 0.25), a = 0, and σ = 0.4. This produces a quartic algebraic Stieltjes transform, so in the numerical fit we use (degm , degz ) = (4, 1) and impose only the zeroth-moment normalization constraint. Figure I.3 visualizes the four sheets and their analytic continuations. As in the compound free Poisson example, the individual branches are discontinuous across their cuts, but the appropriate half-plane gluings form continuous analytic continuations on the full spectral curve. We next use this law as a controlled free-decompression benchmark. A matrix realization of size n = 24 K is generated from the free Lévy law above, and a 22 K-sized principal submatrix is used as the initial compressed matrix. The algebraic curve is fitted from the initial empirical spectrum and then evolved by free decompression to the full size. Figure I.4 shows the resulting density evolution. Panel (a) compares the initial empirical histogram with the density recovered from the fitted algebraic curve. Panel (b) shows the free-decompression flow through intermediate sizes, and panel (c) compares the final FD prediction at n = 24 K with the empirical spectrum 55

(a) Initial Density (n = 4K)

ρ(λ)

1.25

(b) Free Decompression

1.00

0.75

0.50

0.25

0.00

0.8

n = 22.0 K n = 22.5 K n = 23.0 K n = 23.5 K n = 24.0 K

1.25

Fitted Empirical Histogram

Prediction Empirical Histogram

0.6 0.015

0.4 0.2

0.00 0

0.000

0.0 0

Figure I.4: Free decompression of the free Lévy benchmark. A matrix realization of size n = 24 K is generated from the quartic free Lévy law, and a 22 K-sized principal submatrix is used as the initial compressed matrix. Panel (a) shows the initial empirical histogram together with the density recovered from the fitted algebraic curve. Panel (b) shows the FD evolution through intermediate sizes. Panel (c) compares the final FD prediction with the empirical histogram of the full matrix. The initially single-bulk density splits into three bulks, including a small rightmost component that is accurately recovered by the evolved spectral curve. of the full matrix. The initial density has one bulk, while the decompressed density splits into three bulks. The two right bulks arise from the compound free Poisson component, and the smallest rightmost bulk is resolved despite having much smaller mass and height than the dominant component. The same example also illustrates the direct evolution of spectral edges. Rather than locating edges from the reconstructed density at each time, we evolve the branch-point conditions directly as described in Appendix G.1. Figure I.5 shows the edge trajectories from n = 22 K to n = 24 K. The density field computed by free decompression is shown in the background, while the black curves show the independently evolved spectral edges. Starting from one initial bulk, two cusp events create four additional edges, producing three bulk intervals at the final size. The agreement between the black edge curves and the boundary of the nonzero-density region confirms that the edge evolution and density evolution are consistent.

I.3

Pennington–Bahri Model

The Pennington–Bahri model [36] is a random-matrix model for Hessian matrices of single-layer neural networks with ReLU activations under simplifying assumptions. In this model, the network is assumed to be trained to mean-squared error training loss n2 ε over m data points, where n2 is the number of data labels. Writing J for the network Jacobian, the Hessian decomposes into two additive components H = H0 + H1 , where H0 = JJ⊺ /m is the Jacobian component and H1 is the network Hessian term. Under the simplifying assumption that the entries of J and H1 are independent Gaussian variables, H0 and H1 have Wishart and Wigner limiting distributions, respectively. Consequently, the limiting Hessian distribution can be written as an additive free convolution νH = νH0 ⊞ νH1 . Equivalently, its R-transform is the sum of the R-transforms of the two components, RH (w) = RH0 (w) + RH1 (w) =

1 + 2εw, 1 − λw

where λ is the aspect ratio of the Jacobian. This places the Pennington–Bahri law within the finite-atom free Lévy family of Appendix I.2: it corresponds to zero drift, one compound-free-Poisson support point, and a nonzero semicircular component. In particular, its Stieltjes transform satisfies a cubic algebraic relation. To perform decompression of this model, we set ε = 0.1 and λ = 34 , and sampled a 6000 × 6000 matrix from the Pennington–Bahri law. From this matrix, a 2000 × 2000 submatrix was uniformly sampled, and 56

Evolution of Spectral Edges

24 K

n(τ )

10−1 23 K

ρ(λ, τ )

100

10−2 Spectral edge Cusp point

22 K −2

10−3

Figure I.5: Evolution of spectral edges for the free Lévy benchmark. The vertical axis is the matrix size n(τ ), and the horizontal axis is the spectral coordinate λ. Black curves show spectral edges evolved directly from the branch-point equations in Appendix G.1, while the background color shows the density ρ(λ, τ ) computed by free decompression on a fine time grid. The initial spectrum has one bulk; two cusp events create new pairs of edges, leading to three bulk intervals at the final size. The evolved edges align with the boundary of the nonzero-density region, showing consistency between direct edge evolution and full density evolution. free decompression was applied to recover an estimate of the spectral distribution of the original matrix. Figure I.6 shows the result and should be compared with Figure 1 of Pennington & Bahri [36]. This example demonstrates the potential of free decompression as a tool for estimating spectral properties of large Hessian models from smaller submatrices that can be explicitly formed or diagonalized.

I.4

Hessian of Two-Layer Network

Suppose that X ∈ Rn×p is a matrix of n inputs in Rp . Let f (w, X) ∈ Rn×m denote the network output for inputs X with weights w = [wi ]N i=1 . To investigate the spectral behavior of Hessian matrices for f at initialization, we let w∗ denote a target collection of weights comprised of i.i.d elements generated from a normal distribution N (0, σ 2 ). The choice of i.i.d elements here is to mimic standard initialization schemes, including the popular He initialization [19]. Let E ∈ Rn×m denote a collection of residual errors, where Eij ∼ N (0, 2ϵ) are also independently generated for some ϵ > 0. The Hessian matrix to be considered is given by H = [Hij ]N i,j=1 where Hij =

1 ∂2 ∥f (w, X) − Y∥2F , 2n ∂wi ∂wj

Y = f (w∗ , X) + E.

Recall that H is symmetric and therefore has real eigenvalues. However, H is not necessarily positive-definite, and indeed, often possesses negative eigenvalues corresponding to directions of descent. The fraction of negative eigenvalues is called the normalized index and is a prominent object of study in random matrix analyses of deep learning models [36]: index =

#{λi (H) < 0} , N

H ∈ RN ×N .

At scale, when the Hessian becomes so large that it cannot be stored in memory, the index becomes a prime target to be computed using free decompression. As a concrete example, we consider a two-layer feedforward neural network with a ReLU activation function where f (W1 , W2 , X) = ϕ(XW1 )W2 , W1 , W2 ∈ Rp×p , where ϕ(x) = max{0, x} is the ReLU activation function applied elementwise. Our input data X is constructed from the CIFAR-10 dataset [23]. Let Xraw ∈ R60000,1024 represent the total batch of 60,000 training images 57

(a) Initial Density (n = 2K)

0.8

(b) Free Decompression

0.8

Fitted Empirical Histogram

ρ(λ)

0.6

0.4

0.2

0.6

n = 2.0K n = 3.0K n = 4.0K n = 5.0K n = 6.0K

Prediction Empirical Histogram

0.5 0.4 0.3 0.2 0.1

0.0

0.0 −1

−1

Figure I.6: Free decompression for the Pennington–Bahri Hessian model. A matrix A ∈ Rn×n with n = 6000 is generated from the Pennington–Bahri law with λ = 3/4 and ε = 1/10, and a submatrix of size n0 = 2000 is used as the initial compressed matrix. Panel (a) shows the empirical spectrum of the submatrix together with the fitted density. Panel (b) shows the free-decompression flow from n0 to n. Panel (c) compares the final FD prediction with the empirical spectrum of the full matrix. from CIFAR-10 of size 32 × 32 rendered in grayscale and flattened. To control the dimensionality, we truncate this matrix to its n × p upper-left corner, called X̃. In our tests, we use n = 512 and p = 128. This matrix is centered by subtracting the feature-wise mean: 1 Xc = X̃ − X̃11⊺ , p where 1 is a p-dimensional vector of ones. The empirical covariance matrix of this centered data is given 1 by Σ = p−1 Xc X⊺c . To ensure the data is isotropic as per Pennington et al. [37], we apply Zero Component Analysis (ZCA) whitening. Using the eigendecomposition Σ = VΛV⊺ , we construct the ZCA whitening matrix YZCA = VΛ−1/2 V⊺ , where the inverse square root of the eigenvalue matrix Λ is bounded below for numerical stability. The final network input is the whitened data X = YZCA Xc . Finally, we choose σ 2 = 1/p, and ϵ = 0.1. Figure I.7 shows the normalized Hessian index estimated by free decompression across matrix sizes, together with the empirical index computed from direct eigendecomposition at the available sizes. Table I.1 complements this comparison by reporting runtime and scalar accuracy metrics. Since the FD output is a continuous spectral measure, whereas the empirical spectrum is discrete, we use distributional metrics that compare probability measures directly. In particular, the normalized Wasserstein distance W1 /L provides a scale-free discrepancy that accounts for the geometry of the spectral axis and does not require an additional KDE or binning convention [17, 34]. The Hessian’s index error remains small across all sizes, while the direct eigendecomposition time increases rapidly with O(n3 ). In contrast, free decompression only requires the eigendecomposition of the initial submatrix (at size n0 = 212 ), followed by a comparatively inexpensive evolution step. For this empirical Hessian example, we fit a polynomial with (degm , degz ) = (8, 4) and impose moment constraints up to degree two.

I.5

Details of the Diffusion Model Example

The diffusion model that we consider is based on a random features neural network model, which is by now well studied [24] and, in the diffusion-model setting, was introduced in George et al. [16]. A key advantage of this model is that, despite its heavy simplification, it still exhibits memorization–generalization phase transitions that are characteristic of full diffusion models with modern architectures [6, 49]. In this simplified setting, the spectral distribution of the expected activation matrix is tied to the transition between early training, generalization, and memorization behavior [10]. Our goal is to demonstrate that this spectral 58

Hessian Index (Negative Eigenvalues Fraction) Estimated Empirical

Index (%)

15.0%

10.0%

5.0%

0.0% 22 K

23 K

24 K

25 K

Figure I.7: Normalized Hessian index across matrix sizes. The curve shows the index estimated by free decompression, while markers show the empirical index computed directly from empirical eigenvalues at the available sizes. Table I.1: Comparison of direct empirical eigenvalues versus free decompression (FD) for the Hessian example. The Direct column reports the process time for computing the eigenvalues of the Hessian at each size. In the FD column, the first term is the time to compute eigenvalues of the initial 212 -sized submatrix, and the second term is the free-decompression evolution time. The index error is the absolute error in the normalized Hessian index, namely the fraction of negative eigenvalues. KS denotes the Kolmogorov–Smirnov distance between the empirical spectrum and samples drawn from the FD reconstruction. The normalized Wasserstein distance W1 /L is the Wasserstein-1 distance divided by the empirical spectral width L = λmax − λmin at each size. The last two columns report relative errors of the first moment µ1 and standard deviation σ. All accuracy metrics are multiplied by 100 and reported as percentages. Size

Process Time (sec)

Direct

FD (ours)

Index Error

W1 /L

∆µ1 /µ1

∆σ/σ

0.9 91.9 271.9 759.4 2154.4 4755.3 14023.5

0.9 + 0.0 0.9 + 4.7 0.9 + 5.4 0.9 + 7.8 0.9 + 11.1 0.9 + 15.7 0.9 + 22.2

0.0% 0.0% 0.1% 0.6% 0.8% 1.3% 0.4%

0.4% 1.2% 3.0% 2.9% 2.4% 2.1% 2.8%

0.0% 0.3% 0.7% 0.7% 0.5% 0.3% 0.3%

0.3% 3.5% 8.5% 9.0% 6.8% 2.1% 3.5%

0.3% 6.1% 0.9% 0.1% 2.0% 4.9% 13.0%

12.0

2 212.5 213.0 213.5 214.0 214.5 215.0

Accuracy Metrics

Moments Rel. Error

distribution, even when highly multi-scale and multi-bulk, can be accurately estimated from a sub-sampled model. To this end, a 64,000 × 64,000 matrix was generated by sampling m = 800 i.i.d data-points xj , with each j x ∼ N (0, Id ) for d = 100. The random features were generated by sampling 64,000 i.i.d random feature weights w ∼ N (0, Idf ), and computing the expected diffusion activation matrix " ! !⊺ # m 1 X Wxjt Wxjt √ √ U= E σ σ . m j=1 d d Here, the nonlinearity σ is taken to be the hyperbolic tangent function, and p xjt (ξ) = e−t xj + ∆t ξ, ξ ∼ N (0, Id ), denotes the Ornstein–Uhlenbeck evolved sample at diffusion time t = 0.01. The expectation in U is approximated by Monte Carlo using 1024 samples. 59

Table I.2: Atom-like behavior of the leftmost bulk in the diffusion-model spectrum. The empirical mass is the mass of the leftmost bulk detected from the histogram, while the atom-law column reports the prediction w(n) = 1 − (1 − w0 ) n0 /n from Proposition G.5, with w0 = 0.8 and n0 = 22 K. The center-of-mass column reports the empirical center of the same bulk, and its error is measured relative to the reference center after the leftmost component has separated from the neighboring bulk. Size n 22 K 23 K 24 K 25 K 26 K

Left Bulk Mass

Left Bulk Center of Mass

Empirical

Atom Law (ours)

Error

Empirical

Error

80.0% 88.8% 94.4% 97.2% 98.6%

80.0% 90.0% 95.0% 97.5% 98.8%

0.0% 1.3% 0.6% 0.3% 0.2%

1.87 × 10−3 1.80 × 10−3 1.80 × 10−3 1.81 × 10−3 1.81 × 10−3

3.7% 0.0% 0.1% 0.4% 0.4%

We fit an algebraic spectral curve P (z, m) = 0 to the empirical Stieltjes transform of the n0 = 22 K submatrix, using degrees (degm , degz ) = (7, 5) and moment constraints up to degree four. The fitted curve is then evolved by free decompression from n0 = 22 K to n = 26 K. For the edge diagnostics, we do not infer edges from the reconstructed density alone; instead, we evolve branch-point candidates directly from the spectral curve using Proposition G.1, and identify the corresponding physical spectral edges as described in Appendix G.1. The evolution detects one cusp, corresponding to the splitting of the leftmost component into a spike-like bulk and a neighboring small-eigenvalue bulk, consistent with the cusp mechanism in Proposition G.2. The spectrum of U has a pronounced multi-bulk structure over several orders of magnitude. In particular, the leftmost component is an extremely narrow, spike-like bulk carrying most of the spectral mass: it contains about 80% of the mass at n = 4000, increasing to about 98% at n = 64,000, while occupying only a very small interval near the origin. Importantly, in our computation this component is not modeled as an atom; it is treated as an absolutely-continuous bulk and evolved by the same free-decompression procedure as the rest of the spectrum. Nevertheless, its mass follows the atom-law prediction of Proposition G.5 with high accuracy. Table I.2 reports this comparison, together with the center of mass of the leftmost bulk. We also compare the locations of the bulk edges detected from empirical histograms with the edges evolved by free decompression; see Figure 6. Since the spectrum spans several orders of magnitude, edge errors are measured on the logarithmic λ-axis. At the initial size n = 22 K, the two leftmost components have not yet separated, while at larger sizes they split into a spike-like leftmost bulk and a neighboring small-eigenvalue bulk. The only edge excluded from the comparison is the left edge of the spike-like bulk. This edge is governed by finite-size extreme-eigenvalue effects rather than by the macroscopic free-decompression flow: initially it follows the Tracy–Widom-type behavior analyzed in Appendix H and illustrated in Figure H.1. The remaining edges, which correspond to the macroscopic bulk boundaries, are reported in Table I.3.

Appendix J

Implementation and Reproducibility Guide

The algorithms in this paper are implemented through the freealg.AlgebraicForm interface of the opensource Python package freealg.4 The package also includes several closed-form free-probability distributions, such as compound free Poisson and free Lévy laws, which are useful for synthetic experiments and validation. Listing 1 shows a minimal setup for a compound free Poisson example: it constructs the distribution, visualizes its spectral-curve branches, generates a large matrix realization, and then subsamples a smaller principal submatrix. We found that error induced by sampling different submatrices was negligible, being orders of magnitude smaller than the deterministic errors reported in our tables. The synthetic distribution in Listing 1 is used only to provide a controlled demonstration. In applications, the user may instead start directly from a different input matrix, as in the empirical Hessian and diffusion4

freealg is available for installation from PyPI (https://pypi.org/project/freealg), the documentation can be found at https://ameli.github.io/freealg, and the source code is available at https://github.com/ameli/freealg.

Table I.3: Errors of the evolved bulk edges for the diffusion-model spectrum. The entries compare empirically detected histogram edges with the corresponding edges evolved by free decompression, measured on the logarithmic λ-axis. At n = 22 K, the two leftmost components have not yet separated, so the entries under Bulk 2 are omitted; after the split, Bulk 1 denotes the spike-like leftmost component and Bulk 2 denotes the neighboring small-eigenvalue bulk. The left edge a1 of the spike-like bulk is excluded because it is governed by finite-size extreme-eigenvalue effects rather than by the macroscopic free-decompression flow; see Appendix H and Figure H.1. Size

Bulk 1

Bulk 2

Bulk 3

Bulk 4

22 K 23 K 24 K 25 K 26 K

0.0% — — — —

0.0% 0.3% 0.1% 0.8% 1.2%

— 0.1% 0.2% 0.7% 0.4%

— 0.3% 0.2% 0.3% 0.2%

0.0% 0.2% 0.1% 0.1% 0.0%

0.2% 0.5% 0.2% 0.6% 0.8%

0.0% 0.2% 0.2% 0.1% 0.1%

0.1% 0.1% 0.1% 0.1% 0.0%

model examples in Appendix I.4 and Appendix I.5. Given only the smaller matrix Asmall , Listing 2 fits an algebraic relation P (z, m) = 0, estimates the corresponding density, and applies free decompression to infer the density and atoms at the larger matrix size n = 8000. The same interface can also evolve densities over an array of matrix sizes and track lower-dimensional spectral features such as edges, cusps, and atom masses; the atom and edge evolution for this example is shown in Figure 1. Notebook files that reproduce the figures appearing in this work can be found in the codebase. Further details on function arguments and class parameters are available in the package documentation. All experiments were performed on a consumer grade machine, with an AMD Ryzen 7 5800X processor, 128GB DDR4 RAM, running Fedora Linux. Listing 1: Synthetic compound free Poisson setup using freealg. The distribution object is used only to generate a controlled matrix example and visualize the associated spectral-curve branches. # Install freealg with "pip install freealg" import freealg as fa import numpy as np 5

# Create an object for the Compound Free Poisson distribution (see Appendix I.1) dist = fa.distributions.CompoundFreePoisson(t=[2, 5.5], w=[0.75, 0.25], lam=0.1) # Visualize all branches and their gluing as in Figure I.1 x = np.linspace(-1, 9, 200) y = np.linspace(-3, 3, 200) dist.plot_branches(x, y) # Generate a matrix realization of size n = 8000 for this distribution A_large = dist.matrix(size=8000)

# Sample the above matrix to a smaller submatrix of size n0 = 1000 A_small = fa.submatrix(A_large, 1000)

Listing 2: Free-decompression workflow using freealg.AlgebraicForm. Starting from the smaller matrix Asmall , the code fits an algebraic spectral curve, estimates the density, decompresses to a larger matrix size, and tracks edges and cusps. # Create an algebraic-form object for the small matrix af = fa.AlgebraicForm(A_small)

# Fit the polynomial P with degm (P ) = 3 and degz (P ) = 1 (see Algorithm D.1) af.fit(deg_m=3, deg_z=1) # Estimate the density of the small matrix from P using Algorithm E.1 x = np.linspace(0, 10, 200) rho_small = af.density(x, plot=True)

# Free decompression of density and atoms corresponding to a larger matrix of # size n = 8000 using Algorithm F.1 (see Figure 3 ). rho_large, atoms_large = af.decompress( size=8000, x=x, kind=’free’, return_atoms=True, plot=True) 15

# Evolve spectral edges from the small to the large matrix size using Algorithm G.1 # Also this find cusp points where the bulks split (see Figure 1) sizes = np.linspace(1000, 8000, 100) t = np.log(sizes / 1000) edges, cusps = af.edge(t, kind=’free’)

Record · ID 155318 · SHA-256 87bc128a97b6a4d1
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.