Spectral Estimation with Free Decompression
Authors: Siavash Ameli
Chris van der Heide
Liam Hodgkinson
Michael W. Mahoney
Venue: NeurIPS 2025 spotlight ()
Keywords: Random Matrix Theory
Free Decompression
Spectral Density Estimation
Primary area: probabilistic_methods
TL;DR: We use random matrix theory to estimate the spectral density of matrices too large to fit into memory.

Abstract: Computing eigenvalues of very large matrices is a critical task in many machine learning applications, including the evaluation of log-determinants, the trace of matrix functions, and other important metrics. As datasets continue to grow in scale, the corresponding covariance and kernel matrices become increasingly large, often reaching magnitudes that make their direct formation impractical or impossible. Existing techniques typically rely on matrix-vector products, which can provide efficient approximations, if the matrix spectrum behaves well. However, in settings like distributed learning, or when the matrix is defined only indirectly, access to the full data set can be restricted to only very small sub-matrices of the original matrix. In these cases, the matrix of nominal interest is not even available as an implicit operator, meaning that even matrix-vector products may not be available. In such settings, the matrix is "impalpable", in the sense that we have access to only masked snapshots of it. We draw on principles from free probability theory to introduce a novel method of "free decompression" to estimate the spectrum of such matrices. Our method can be used to extrapolate from the empirical spectral densities of small submatrices to infer the eigenspectrum of extremely large (impalpable) matrices (that we cannot form or even evaluate with full matrix-vector products). We demonstrate the effectiveness of this approach through a series of examples, comparing its performance against known limiting distributions from random matrix theory in synthetic settings, as well as applying it to submatrices of real-world datasets, matching them with their full empirical eigenspectra.

Spectral Estimation with Free Decompression
Siavash Ameli
Chris van der Heide
ICSI and Department of Statistics University of California, Berkeley sameli@berkelln edu edu
Dept. of Electrical and Electronic Eng. University of Melbourne chris. vdh@gmaill.com
Liam Hodgkinson
Michael W. Mahoney
School of Mathematics and Statistics
ICSI, LBNL, and Department of Statistics University of Califomia, Berkeley
University of Melbourne
Ihodgkinsondunimell edu.
mmahoney@stat.c berkeley edu
Abstract
Computing eigenvalues of very large matrices is critical task in many machine learning applications, including the evaluation fog-determinants, the trace ma- trix functions, and other important metrics. As datasets continue to grow in scale, the corresponding covariance and kernel matrices become increasingly large, often reaching magnitudes that make their direct formation impractical Existing techniques typically rely on matrix-vector products, which can provide efficient approximations, if the matrix spectrum behaves well. However, in settings like distributed learning, or when the matrix is defined only indirectly, access to the full data set be restricted to only very small sub-matrices of the original matrix. In these cases, the matrix of nominal interest not even available as an implicit operator, meaning that In such settings, the matrix is "impalpable," in the sense that we have access to only masked snapshots of it. We draw principles from free probability theory to introduce a novel method of "free decompression" to estimate the spectrum of such matrices. Our method can be used extrapolate from the empirical spectral densities of small submatrices to infer the eigenspectrum of extremely large (im- palpable) matrices (that we cannot form even evaluate with full matrix- vector products). We demonstrate the effectiveness of this approach through a series of examples, comparing its performance against known limiting distributions from random matrix theory in synthetic settings, well as applying it to submatrices of real-world datasets, matching them with their full empirical eigenspectra.
impossible.
matrix-vector products may not be available.
Introduction
Eigenvalues encode essential information about matrices and are central to the evaluation of spectral invariants in computational linear algebra. Two such invariants, the log-determinant and the trace of the inverse, appear frequently in statistical modeling and Bayesian inference (Ubaru & Saad, 2017). However, these quantities require access to the full range eigenvalues, posing substantial challenges in modern large-scale settings. As datasets continue scale up, it is becoming increasingly common to encounter matrices that cannot practically formed in their entirety. This may occur, for example, if the matrices are so large that they exceed the storage capacity in local memory. Popular methods for operating with large matrices often rely on randomized low-rank approximations from Randomized Numerical Linear Algebra (RandNLA) (Drineas & Mahoney, 2016; Murray et al., 2023; Dereziński Mahoney, 2021), and can broadly be classified as either iterative subspace or sampling methods, although hybrid methods have become prominent (Ubaru et al., 2017; Saibaba et al., 2017: Chen
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

et al., 2021). Implicit methods (Adams et al., 2018) that rely on Krylov subspaces using Amoldi or Lanczos-based techniques (Ubaru et al., 2017) have been successful for estimating arbitrary spectral invariants (Gardner et al., 201 8; Potapczynski et al., 2023), and they require access only to matrix- vector products. These techniques avoid explicit matrix storage and can often converge quickly. However, their performance degrades significantly when the matrix is ill-conditioned.
The situation is particularly dire when most matrix entries are inaccessible and matrix-vector prod- ucts are unavailable or too costly. This can occur when the data used to generate the matrix is spread across a distributed system, or when the matrix itself is enormous (see e.g., Ameli et al. (2025)). Such matrices can be considered impalpable, because neither the matrix nor its matrix- vectorr products can be handled in their entirety. Despite being inaccessible, the information these impalpable matrices contain can be crucial, making methods to summarize or access them highly valuable. Approaches which subsample a few rows and columns are feasible, assuming that any observed quantities can be directly extended to the full matrix (Drineas et al., 2008). This is the strategy behind the Nyström method, which accesses a matrix via column/row sampling, and then reconstructs the matrix under linearity assumptions (Gittens& & Mahoney, 2016; Gittens al., 2016). Unfortunately, such tech- niques incur significant bias, as they target behavior dictated by large eigenvalues by design, giving less prominence to the near-singular dimensions that play crucial role in determinants and inverses.
In the setting of impalpable matrices, or when otherwise fine- fine-grained behavior of ill-conditioned matrices becomes important and implicit methods fail, there are very few tools currently available (Couillet & Liao, 2022; Ameli et al., 2025). To address this gap, we provide a novel method for approximating the empirical spectral density of extremely large (impalpable) matrices, requiring only access to appropriate submatrices. Our technique is based on tools from random matrix theory and free probability, and it applies to virtually any class of matrices. Assuming only that the underlying matrix is sufficiently large, the empirical spectral density of a (sampled) submatrix can be evolved by a partial differential equation (PDE) in terms of the dimension of the matrix to access the spectral density of the full (impalpable) matrix. Due to connections with free probability theory, which we explain below, we call this process free decompression.
To summarize, we:
(I introduce and derive free decompression, the first method for spectral density estimation of impalpable matrices, by evolving the empirical spectral density of submatrices, and provide rigorous polynomial error bounds under mild assumptions (Proposition 4);
(II) demonstrate the utility of our free decomposition method for a range of synthetic examples, social media graph data, and empirical neural tangent kernels; and
(III) provide a uree-friendly Python package, freealg,' implementing all proposed algorithms, that reproduces the numerical results of this paper.
The remainder of this paper is structured as follows:
Section 2 provides background and outlines the different classes of matrices of interest; Section provides the necessary tools order to construct our free decompression procedure; Section 4 describes its implementation and issues faced in oper- ationalization; Section provides numerical examples. We conclude with diseussion in Section 6. A summary of the notation used throughout the paper is provided in Appendix A. An overview of free decompression and its application to well-studied random matrix ensembles is presented in Appendix B.1. Discussions of the significant numerical considerations required to perform free decompression are provided Appendix C and Appendix D. Proofs of our theoretical contributions can be found in Appendices E and F. Background on the neural tangent kernel example is provided in Appendix G. Finally, Appendix H presents guide to our implementation of free decompression.
2 Impalpable Matrices
We are interested in estimating the spectral density of wide class of matrices for which operating on the full matrix is infeasible. In the most extreme case, these matrices cannot be formed in their entirety, either implicitly or explicitly. We call these matrices impalpable. We will consider the case where the impalpable matrix of interest is Hermitian. Since we are interested in the setting where we
'The source code of the Python package freealg is available at https://githh:ii com/ ame /mtee//eeaal. with documentation at https: //ameli github. 1 offfreeagg.
2

---

cannot access the full matrix, we will work with submatrices indexed by randomly sampled columns of the impalpable matrix of interest, where the submatrix consists the intersection of the impalpable matrix's rows and columns that appear in this index set. This corresponds to an assumption that the matrices are part of sequence that asymptotically free (Maïda, 2023). We will then see that it is convenient to work with tools from random matrix theory and free probability, namely the familiar Stieltjes transform (Couillet & Liao, 2022), whose favorable properties we exploit in our free decompression procedure.
2.1 Categories of Matrix Difficulty
We find it useful to classify matrices by four modes of computational difficulty: explicit, implicit, out-of-core, and impalpable, the latter of which is our focus.
Explicit Matrices. This is the most straightforward class of matrices. They can be formed explicitly as contiguous arrays and operated on in memory. Existing implementations on such matrices are highly optimized, with their runtime limited only by their computational complexity.
matrices become too large to be stored in memory. Implicit matrices are those that can be computed or estimated only in terms of matrix- vector products v + Av (Adams et al., 2018). Treating matrix A E R as implicit often has significant computational benefits, even before the memory wall (Gholami et al., 2024); for iterative methods, computation time can often be reduced, e.g.. from O(n²) operations O(n²). If the formation of the matrix-vector product is efficient, this can be further reduced to (n). Memory requirements may follow similar reductions. Implicit methods exist for most standard linear algebra operations, often based on Arnoldi- or Lanczos-baser iterations, and appear mature software packages (Potapczynski et al., 2023). One of the most popular implicit matrix algorithms for spectral function estimation is stochastic Lanczos quadrature (SLQ) (Ubaru et al., 2017); however, the performance of such methods may deteriorate for highly ill-conditioned matrices (Ameli et al., 2025).
Implicit Matrices. Difficulties begin arise
Out-of-Core Matrices. Implicit methods typically operate on Krylov iteration schemes, with error rates depending on condition numbers of the matrix-sed Bhattacharjee et al. (2025), e.g., in the case of SLQ. This can be disastrous for large and ill-conditioned dense matrices, with accurate estimates requiring a large number of matrix-vector products. Exact methods become necessary in this setting if a high level of accuracy is required, but the memory bottleneck still remains. Out-of-core matrices have the property that for some integer and operated on memory, and the full matrix can be stored locally by iterating through these blocks (Aggarwal & Vitter, 1988). Decompositions (Rabani & Toledo, 2001; Mu et al., 2014), inverses (Caron & Utard, 2002), and products of out-of-core matrices can be (slowly) computed to high precision, assuming adequate (typically very large) available external storage.
n,, any subblock of size R m= can be formed explicitly
Impalpable Matrices. The most ex- treme category of matrices includes cases that do not fit into the other categories. Here, the matrix cannot be formed in mem- ory, is either core operations, or has missing portions. Matrix-vector products assumed too pensive to compute accurately, or are ren- dered useless due strong ill-conditioning. In principle, such a matrix can be consid- ered "impalpable." of either its elements or properties cannot be operated on. Comprising the worst-case scenario, impalpable matrices offer unique challenges for linear algebra practitioners. Determining accuracy of any proposed solution is an underdetermined problem in general. The only feasible strategy is extrapolation, using available information. While this is inherently risky, the class of impalpable matrices simply offer no alternative. To our knowledge, only the FLODANCE algorithm of Ameli et al. (2025) offers valid approach for estimating log-determinants of impalpable matrices.
Access in Memory
Matrix-Vector
Any Subblock
Type
Matrix
Product
large for out-of-
Explicit Implicit
Out-of-core
X
Impalpable
X
X
Table 1: Comparison of matrix classes, according to at- tributes that can be formed and operated on memory.
the vast proportion
3

---

2.2 Examples of Impalpable Matrices
Enormous Gram Matrices. Our motivating example of an impalpable matrix is an ill-conditioned Gram matrix A (M(2), B3 )) comprised of a massive number n of points æ; and kernel that is expensive to compute. In this case, computing few subblocks of A is viable, but out- of-core methods are too expensive. Ameli al. (2025) demonstrates the challenges of computing the log-determinant of a neural tangent kernel matrix for large neural network model. Their most cost-effective strategy extrapolated the behavior of the log-determinant from smaller submatrices in accordance with scaling law. This is an example of an impalpable matrix strategy, relying on smaller portions of the matrix without matrix-vector products. The potential gains are best demon- strated with synthetic example: the covariance matrix of 2...... 2n where each 2; ~ *N(O,p. I pxp) p-dimensional. Once and become extremely large, forming this matrix is prohibitive. However, its spectral density is well-approximated by the Marchenko-Pasun law (Marčenko & Pastur, 1967), and its ge-deteemminnt over submatrices of increasing size is predictable (Nguyen Vu, 2014; Cai et al., 2015; Ameli et al., 2025; Hodgkinson et al., 2023a).
Catastrophic Cancellation. Impalpable matrices typically arise due to issues of scale, but they need not be large themselves. Indeed, such matrices can arise due to catastrophic cancellation. Some matrices can be comprised of a non-trivial product of rectangular matrices; for example, cond( A A) = cond(A) A) so representing the outer product A T (e.g. in the normal equations) in floating-point precision can incur greater numerical error (Ameli et al., 2025). Higham (2022) refers to the formation of A as a "cardinal sin" of numerical linear algebra. More products yield even larger errors. General Gram matrices exhibit similar behavior: the condition number can become so large that the matrix cannot be accurately represented in double precision, rendering even exact methods hopeless (Lowe al., 2025). Fortunately, submatrices have smaller condition numbers, so there may be a submatrix that can be represented double precision. Ideally, spectral behavior of these submatrices then be extrapolated to obtain insight into the behavior of the full matrix in meaningful way.
3
Free Decompression
Our objective is to approximate the spec- tral distribution of an impalpable matrix A, given only explicit access to a randomly sampled n X n submatrix. This seemingly intractable task can be accomplished by working with the Stieltjes transform of the submatrix, and observing that for large this Stieltjes transform varies with n in a manner that is well approximated by a PDE. This enables an approximation of the Stieltjes transform of the full matrix A, from which a corre- sponding spectral distribution is easily computed. We emphasize that the matrix A itself may be arbitrary and deterministic. Asymptotic freeness enters our model through the method sumbatrix selection, making it amenable to random matrix theoretic and free probabilistic tools.
Free decompression of random submatrix A, to a larger matrix requires:
1. estimation of its Stieltjes transform mA,
2. evolution of mA, in via equation (I);
3. evaluation of the spectral distribution of A.
3.1 Tools from Free Probability
In order to make free decompression concrete, we require tools from free probability. Let A be fixed matrix of large size. If we can approximate the spectral density of a small randomly sampled submatrix, and extrapolate its behavior to the full matrix, it should also be possible to extrapolate further to larger matrices than the one we are given. To this end, we treat A as an element in an infinite sequence of matrices Ag........ of increasing size, where A,, " For consistency, we assume that for each n, the matrix A, is the top-left n x n submatrix of A This ensures that [Anl [A)y is constant in n for fixed j,j S n. We will also assume that for any integer k > 1, the quantity [rr(AA)) converges as n co. This is equivalent assuming that the empirical spectral distributions of A, converge weakly as n > de (Voiculescu, 1991). Note that this imposes no requirements the matrix to be approximated.

---

[a] Initial Empitioal Dersity (n IK)
(b) Pree Decompression
[c) Final Empirical Density (a 32K)
Dexatty Satinate
Bendurrack Dendty (döaatt)
2.0
Regiiical Spectoon
4X
Fise Descspecau Soluikon
Denelty ((
RX
1.5
lok
82k
1.0
0.5
0.0
2
3
2
Figure 1: (a) Estimation of the Marchenko- law with ratio = 50 from eigenvalues Wishart matrices with d 50,000 degrees of freedom and size n. x n 1000 1000; (b) approximation of densities for various values of n via free decompression (Algorithm 1); and (c) comparison of our decompressed approximation of 32,000 with the known limiting analytic density for A 30
Still, we cannot predict the upcoming row and column in An+1 from A, alone, and we will need to appeal random transformations connect each element the sequence together. One way to transform A into random matrix, while preserving its eigenspectrum, is to perform similarity transform with random matrix P. recovering PAP- -1 The most inexpensive choice of random matrix P is a (Haar) random permutation matrix, so that the (random) similarity transform is equivalent to the following operation: let S,, n be the uniform distribution on the space of permutations on {1, . n. and consider the matrix A, (Ao(t)o())):=1. where a S,..
Free Compression. Due result of Nica (1993), the sequence of matrices A is asymptotically free as n + oc, in the sense of free probability theory. We shall not further discuss freee random matrices here, but we refer the interested reader to Maïda (2023) for a detailed introduction. Instead, we rely on one useful property.
m m (z) denote the Stieltjes transform of A of size n nn, given by m (2) = n - E[t(AA)) and let (2) denote the functional inverse of m(=), satisfying m(w(z)) = 2. The corresponding R-transform is then defined as R(=) w(-z) z_! Our method depends on the following fundamental theorem of Nica & Speicher (1996); see also Olver & Nadakuditi (2012, Section 7).
Theorem UNICA-SPEICIIL
For a free random mairix A of size X with R-transform R(z). the R-transform of the top-left n, Tts submatrix of A given by Rn, (2) R)=
In other words, to arrive at the R-transform for submatrix comprised of only the top-left proportion of rows and columns, one need only scale the argument by a. In free probability theory, this procedure is called free compression (Olver & Nadakuditi, 2012).
Free Decompression. Our main approach is based the observation that nothing obstructs this operation from being conducted in reverse, to glean information about a larger matrix from smaller one. For this reason, we call our procedure free decompression: starting from the R-transform of a smaller n, x submatrix Rm, (z), the R-transform of the larger n X n matrix is given by R, (2) = Rm, (=z)..
Example (COVARIANCE MATRICES). Let X E nxd have elements Xi NN 1) and let A -XXT "x* Assuming d n, then A is almost surely full rank, and is n X n Wishart- distributed with d degrees of freedom. The top-left n, n., corner of A" is also Wishart-distributee with d degrees of freedom; and so, for large n, and d, its eigenvalues approximately follow the Marchenko-Pastur law with ratio À = /d.d. The H-transform is given by Rn, (z) = (1-2n./d 2n./d) Applying free decompression, we estimate the R-transform of the full matrix be
R₂ = (1
= 1-22 n -1---)))))----
R2 100=(11-2
72 n2
This is the R-transform of a Marchenko-Pastee law with ratio = n/d, and it approximates the spectral density of Wishart-d matrix with degrees of freedom and size x This
5

---

is precisely what we expected for the full matrix A. In effect, this allows us to infer spectral properties of the A= = X E R from the submatrix A., = Xx,,, XT, E R", where the x, E IR xd are (much) shorter matrices with (much) smaller aspect ratio À (see Figure 1).
However, for more general classes of matrices, inverting the Stieltjes transform in order to compute its R- t-trransfomm is not analytically tractable. Instead, we will work with the Stieltjes transform directly.
3.2 Evolution of the Stieltjes Transform
Our approach can be made precise by appealing to Theorem 1 and operating directly on the R transform. We do this in Proposition 2, which shows that free decompression corresponds to the evolution of PDE in the Stieltjes transform, given by (1). To do so, for fixed R-transform, we let R(t, = R(ze"), ), so that a n/n, = e', and we consider the corresponding Stieltjes transform m(t, 2) and its inverse w(t, z). Then, the following proposition, proved in Appendix E, holds.
Proposition 2. The Stieltjes transform m(t, 2) corresponding to R-transforms R(ze") satisfies
dm
Öt = m + m
Since (1) is first-order quasilinear PDE, it can be readily solved using the method of characteristics. This is outlined in Proposition 3, which is proved in Appendix E. A key consequence of this proposition is that we obtain explicit solutions to (1) terms of the initial data and desired "time" t.
Proposition 3. Suppose my is analytic in an open domain Q C C, and mo(z) + for all S. Consider the initial-value problem (1) with m(0, z) mm(z). For zo, E Q there is a unique solution m(t,2) 2) whose graph can be parametrized by curves (I(T), m()))
Bzo where
= =T 2(T) = 20 mo (c- -1- m(T = mo
(2)
valid for such that 2(T) E S. Furthermore, letting Φ z, =0) -- 20 + m(=0) 11( 1), so ((t(T), which allows the explicit solution of (1) by m(t, 2) = mo((tt =))e
-==0))1(()) O, there exists an inverse function 20 = 2t,, solving a(t, o(t, 2))
Finally, Proposition 4, proved Appendix F, shows that errors in free decompression can grow, but no more than polynomially fast in the decompression ratio _.
Proposition 4. Let m be a Stieltjes transform of a density such that lIpll -0 < +oc, and let 6m,, denote the error in an approximation to starting with matrix size of ng. For 8m., the error under free decompression (1) to a ratio n n/n. there is a constant > O such that for any I >> O,
|lōmmllz S (n//n)"(l0mm, 11
3.3 Heuristic Motivation
While our method is ultimately constructed via Theorem 1, it is illustrative to consider the behavior of the Stieltjes transforms of successive matrices A,, in our sequence. In what follows, we will proceed informally, roughly following the reasoning outlined in Tao (2017). This heuristic approach to aid the reader's intuition by viewing the evolution equations in terms of the Schur complement, to motivate the full proof in Appendix E. We let e, denote the i-th column basis vector. Taking two successive matrices in our sequence A,,, we exploit properties of the Schur complement and the Woodbury matrix identity, to obtain the relationship between the trace of their inverses
tr A n+1 - 21- tr(A,, =I) -I- + e 4 =1) 2 en
(3)
e (An+1 =I) en+1
Assuming the rows and columns of A have already been randomly permuted, then the distribution of each row and column of (An41 =I)- -1 is the same. Therefore, in expectation,
n+1
1
1
=I)
en+1 =I)_1ent]
E
Te (A (An - =I) =1)- e
tr( (A - =1] =1)-
n +
n n+1

---

(a) Principal Branch on Cf and C-
[b] Principal Branch on C*2 Secondary Branch C-
0.5
0.5
Tmja]]
-0.5
-0.5
-1.0
0.5
1.5
2.0
-0.5
0.0
0.5
1.3
Re(=)
Re(=)
Figure 2: Analytic continuation of the Stieltjes transform of Marchenko Pastur distribution (with ratio x 1
50 (a) The principal branch contains branch cut along the support of the spectral density, while (b) the secondary branch continuous in this region. Curves (white) highlight the evolution of (t, over E [0, 4]. crossing the support of the spectral density.
Algorithm 1: Pscudocode for Free Decompression
Input :random submatrix A,, ER,,"" X*, of Hermitian matrix A RR"*" spectral smoothing function EigApprox; gluing function Glue; Stieltjes transform approximation Stieltjes, perturbation
Output
:Estimated spectral density values {
Compute the spectrum of A,, (1))==1
If Step I:
Compute spectrum
2 BigApprox((A)i... (S)e
If Step 2: Estimate smoothed spectral density II Step Estimate glue function I/ Step 4: Compute the Stieltjes transform
3
4 m(=) Stieltjes(S, G).
5 foreach x E do
6
Find 2g E C such that + i0 = 2z m(s.)
If Step 5: Solve characteristics
return m22z) R2 Tex
If Step Return the estimated spectral density
The same holds for e (A,+1 n+1 2I) e This is effectively Hutchinson trace estimation procedure (Bekas et al., 2007), and we note that the variance of the Hutchinson estimators be- come small for large as well KRoosta-Khorasani & Ascher, 2015). Assuming this and noting that ir( 1r(AntI =I)) = Lrr( n+1 =I) .normalizing by + 1 we deduce the difference equation
MA (=) mA, (z) = -nn=nn--mm (- + (z)mia, (z)) On
(4)
mm(e' 2) n|mA, (z) mA, (z)), (4) is well For large, provided e' for that de approximated by the differential equation the Stieltjes transform m. (1, =) given by (1).
4 Approximating the Stieltjes Transform
In principle, Propositions 2 and 3 suggest straightforward approach estimating the Stieltjes transform under free decompression, but these results disguise two main practical challenges:
1. The Stieltjes transform is singular on the support of the spectral density, requiring analytic con- tinuation way that does not disrupt the complex transport flow of (1). Without imposing additional properties, numerical analytic continuation is notoriously ill-posed (Trefethen, 2023). The second challenge is accurately approximating the initial Stieltjes transform m(2) from the submatrix eigenvalues. This well-studied problem, but as we will see, the nature of the PDE (1) necessitates accomplishing this task to a much higher degree of precision than is typically required.

---

4.1 Analytic Continuation
The Stieltjes transform m is typically defined on the upper half-plane C+ = {2 3(2) > 0} and extended to the lower half- C via the Schwarz reflection principle, m(z) = m(z). This defines the principal branch, which exhibits a branch cut along the real axis on I := supp(p). However, the curves ((t, traced by the PDE dynamics descend from C+ into C crossing the real axis precisely through I, where the principal branch becomes discontinuous. This is formalized in the following proposition.
Proposition 5. Let the curve & + ((1,2) =) be the solution of charucteristic curve given in Propo- sition 3. Then, for each z E C+, there exists t" > O such that Bo(t", 2) = 0. In particular. Ro(t" 2) E
Proposition 5 is proven in Appendix C.I, and this crossing necessitates defining a new branch of the transform-equal to the principal branch in C+ -hat extends analytically into C and remains holomorphic across the interior 7". The existence of such continuation is guaranteed by the principle of analytic continuation: since m is holomorphic in C+ and has singularities I, it can be extended across the cut into C We refer this extension the secondary branch.
Figure illustrates this continuation for the Marchenko- approaches were considered to construct the secondary branch, and are discussed in Appendix D. The most successful approach, and the one considered in experiments to follow, introduces gluing function G(z) defined in (C.5), which compensates for the discontinuity across I", In Appendix C, we show that this function can be accurately approximated by low-degree rational function (a Padé approximant). When the analytic form of the Stieltjes transform is known, this construction is exact, as shown in Table B.3. In empirical settings when the density is estimated from the eigenvalues of submatrix- the same approximation still enables accurate continuation.
law, showing both branches. Multiple
In the sequel, all references to the Stieltjes transform m implicitly refer to this secondary, analytically continued branch.
4.2 Approximating the Empirical Spectral Density
In light of Proposition 5, free decompression cannot applied-eve even in a weak sense- to the empir- ical Stieltjes transform, m A, (z) Ei - Instead, we must work with the Stieltjes transform of smoothed approximation to the empirical spectral distribution. For the synthetic examples consider, the spectral densities of interest exhibit square-root behavior at the edges of their support, as shown in Table 1. making Chebyshev polynomials natural basis for approximation.
However, for real datasets, which may exhibit more complex or irregular spectral behavior, is often preferable to use the more flexible class of Jacobi polynomials. These introduce two tunable hyperparameters (a, 3) and recover Chebyshev polynomials when (a, 3) = ( $). To fit these polynomials, we found that Galerkin projections of kernel density estimates were often necessary, using Gaussian or Beta kernels. Once the spectral density is approximated in the Jacobi basis, the Stieltjes transform is computed directly using Gauss-Jacobe quadrature (Shen et al., 201 1, Section 3.2.2). Details of these procedures and a full algorithmic description of our implementation can be found in Appendix D. Pseudocode for the full free decompression method is provided in Algorithm 1.
5 Numerical Examples
We now present numerical examples that demonstrate the utility of our free decompression method. To begin, we consider a number of synthetic examples for random matrices whose spectral densities and their corresponding Stieltjes transforms have known analytic expressions. We then consider two evaluations on real-world scagesscale datasets: Facebook Page-Page network (Leskovec & Krevl, 2014; Rozemberczki et al., 2021); and the empirical Neural Tangent Kernel (NTK) corresponding to ResNet50 (He et al., 2016) 6) trained on CIFAR 10 (Krizhevsky, 2009). All experiments were conducted on consumer-grade device with an AMD Ryzen07 5800X processor, NVIDIA RTX 3080, and 64GB RAM. Hyperparameters used in each experiment are listed in Appendix H.
8

---

[a] Laplacian- Facebook Page-Page
(b) Neural Tangent Kernel
CIFAR-10 10
2.0
21
21
1.5
=
22K
2
((A)
0.0
0.0
0.5
1.0
1.5 2.0 4
10 12 14
A
Figure 3: (a) Empirical spectral densities (dashed) of the symmetrically normalized Laplacian matrix of the SNAP Facebook dataset, compared against estimates obtained via free decompression from an initial submatrix of size 2 (solid). (b) Empirical spectral densities (dashed) of submatrices of the log-NTK matrix of size 50,000, computed from the CIFAR- 10 dataset using a ResNet-50 model, and estimates obtained via free decompression from an initial submatrix of size 20 (solid).
Synthetic Experiments: We consider five families of distributions that commonly arise spectral
distributions in random matrix theory: the Wigner (1955) semicircle law (free Gaussian); Marčenko & Pastur (1967) (free Poisson); Kesten (1959) and McKay (1981) (free binomial); Wachter (1978) (free Jacobi): and the general free-Meixner family (Saitoh & Yoshida, 2001; Anshelevich, 2008) (free analogs of the classical Meixner (1934) laws). In Table B.I, we list, for each law, its absolutely continuous density p(x) supported on [X_, X+], and any point masses. Immediately following, Table B.2 summarizes the simplest madomm-atrie or combinatorial constructions realizing these laws in practice. Figure depicts our free decompression procedure for the Marchenko-Pasee case. We sampled X < R
1000, d 50,000, with each X, NN 1). The left panel shows the empirical density of A = fxxT well as our density approximation using Beta kernels projected onto 50-degree Chebyshev polynomials. The center panel shows range of densities estimated by free decompression. The right panel compares our freely decompressed approximation of the spectral distribution of F==x= n 32,000, with the known exact Marchenko-Pase Pastur law with ratio The principal and secondary branches are compared in Figure 2. The remaining models
32 A = 50 are treated in Appendix I. with comparison of the expected log-determinants and distributional distances.
SNAP Facebook Dataset: We use the publicly available Facebook Page-Page network from SNAP (Leskovec& & Krevi, 2014; Rozemberczki et al., 2021), where nodes represent verified Facebook pages (e.g., politiiii,,,,, musicians) and edges indicate mutual likes. This graph undirected with 22, 270 nodes and 171,002 edges. We study the eigenvalue distribution of the symmetrically normalized Laplacian matrix of the graph (see e.g.. (Mahoney, 2016, 16, Definition 4)), which yields a compact spectrum in [0, 2]. To avoid artificial spikes caused by degree-1 leaf nodes (which concentrate eigenvalues at À = 1), we apply a small Erdös-Rényí Rényi perturbation graph (n.. p) at the connectivity regime where p = Iog(n) (Huang & Landon, 2020) to the adjacency matrix of the graph before constructing the Laplacian. Comparison of the empirical spectral densities of randomly sampled submatrices with their free decompression approximations shown in the left panel of Figure 3.
Neural Tangent Kernel: Finally, we consider the spectrum of an empirical NTK derived from a ResNet50 network trained on CIFAR 10. The NTK is a well-studied object in machine learning (Jacot et al., 2018; Novak et al., 2022); see Appendix G for background. This example is particularly challenging where naïve approaches fail. The empirical NTK is known to be highly ill-conditioned for well-trained multi-class models (Ameli et al., 2025). When the NTK is constructed in blockwise fashion, classification models with e classes often exhibit a low-rank structure, with approximately 1/c of the eigenvalues being extremely small. To address this, we project the NTK matrix onto its non- null eigenspace, removing the near-zero eigenvalues, and then reconstruct Hermitian matrix in the reduced space. Subsamples are drawn from this reduced matrix. Specifically, we took À to be the full 55,560 x 55,50 5500 matrix (corresponding to 5556 images), and reduced this to 50,000 x 50,000 matrix A after removing the null component. Submatrices of dimension 1024, 2048, 4096, 8192, 16,382,

---

Table 2: Comparison of process time for direct computation of the spectral density versus free decompression (FD) on the NTK dataset. In the FD column, the first term is the time to compute eigenvalues of the initial 2 sized submatrix (common to all rows) and the second term is the FD step. Metries (columns 4-8) compare the true spectral density with the FD approximation, and are reported as distributional distances (total variation, Jensen-Shannon, Kolmogorov-Sminat relative errors of the first and second moments, all multiplied by 100 to give percentages. Values are the mean over 20 randomly sampled initial submatrices with standard deviations in parentheses.
and
Size
Process Time (sce)
Distributional Distances (x 100)
Moments Rel. Error (x100)
Direct
FD (ours)
TV
J5
KS
A/1//1
A
2
3.6
3.6 1.6+000 0.0
0.00%
(100.00)
0.00% (1.0.00)
0.00% (1.0.00)
0.00% (1.0.00)
0.00% (1.0.00)
2"
10.2
0.6
1.72%
(1.0.77)
7.60% (11199)
0.48% (10.10)
0.05% (100.3)
0.09% (100.44)
2
50.9
0.6
2.06%
(+0.35)
4.67% (=0.67)
(=0.12)
0.01% (=0.01)
0.02% (#:0.02)
2
358.9
0.6
3.24% (+0.56)
6.30% (=0.51)
1.18
0.01% (+0000)
0.02% (+:0.01)
2
2820.2
3.6+0.7
4.33%
(+0.87)
7.55% (=0.80)
1.76% (=0.48))
0.01% (+0.01)
0.03% (+0.01)
2 20451.2
3.6
3.6+-0.8
0.8
5.16% (11.18)
7.96% (11.11)
2.51% (10084))
0.02% (.10.01)
0.05% (100.2)
50K
67331.1
3.6 1 0.8
5.94%
(11.48)
8,33% 1141)
3.02% (:1.14)
0.17% (100033
0.49% (1.0.66
and 32, 768 were then sampled, and their empirical spectral distributions are shown as dashed lines in the right panel Figure Free decompression approximations based on the Stieltjes transform of the 1024 x 1024 submatrix (see Figure G.1) are shown as solid lines. Table 2 reports the corresponding process times and errors in total variation (TV) distance, Jensen-Shannon Shannon (IS) distance, Kolmogorov- Smirnov (KS) distance, and the relative error in the first and second moments.
6 Conclusion
In modern settings, it is becoming increasingly common that the spectrum information of matrix required, even when the matrix too large fit into memory, is otherwise unavailable. These matrices comprise the category of impalpable matrices, offering unique challenges to numerical linear algebra practitioners. We have outlined free decompression, novel method of approximating the spectral density of large impalpable Hermitian matrix from the spectrum of randomly sampled submatrix. Our initial implementation of this method shows promise a tool extract accurate information about the spectrum of extremely large matrices that would otherwise be unavailable. We hope that this work motivates other researchers to further develop these and similar tools further unlock information in impalpable matrices that is otherwise out of reach.
Limitations. While our method imposes few restrictions on the large matrix of interest, the quality of approximation depends strongly on asymptotic estimates of the spectral density of A,, While promising for future research, performing free decompression accurately remains highly challenging, with Proposition 4 highlighting the need for accurate spectral density estimation and analytic continu- ation. Without general approach to estimate and operate on the Stieltjes transform on log-scale, that is, estimating + m(e") )) highly ill-conditioned matrices still remain out of reach. At present, the implementation of our method is sensitive to the quality of approximation of the gluing function; and, while our approach works flawlessly in simple scenarios, it becomes more delicate with real data. Nevertheless, we stress that despite these limitations, no competing alternative exists for estimating the spectral density of impalpable matrices.
Acknowledgments and Disclosure of Funding
LH is supported by the Australian Research Council through a Discovery Early Career Researcher Award (DE240100144). MWM acknowledges partial support from DARPA, DOE, NSF, and ONR. This work was supported in part by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research's Applied Mathematies Competitive Portfolios program under Contract No. 1232-
10

---

References
Adams, R. P., Pennington, J., Johnson, M. J., Smith, J., Ovadia, Y., Patton, B., Saunder- son, J. (2018). Estimating the spectral density of large implicit matrices. arXiv preprint arXiv:1802.03351.
Aggarwal, A. & Vitter, Jeffrey, S. (1988). The input/output complexity of sorting and related problems. Communications the ACM, 31)),, 1116-1127.
Ameli, S., van der Heide, C., Hodgkinson, L.,. Roosta, F., & Mahoney, M. W. (2025). Determinant estimation under memory constraints and neural scaling laws. In Forty-second International Conference on Machine Learning.
Anshelevich, M. (2007). Free Meixner states. Communications in Mathematical Physics, 276(3), 863-899.
Anshelevich, M. (2008). Orthogonal polynomials with resolvent-type generating function. Trans- actions of the American Mathematical Society, 360(8), 4125 4143.
Baker Jr., G. & Graves-Morris, (1996). Padé approximants, volume of Encyclopedia of Mathematics and its Applications. Cambridge University Press.
Bekas, C., Kokiopoulou, E., & Saad, Y. (2007). An estimator for the diagonal of a matrix. Applied Numerical Mathematics, 57, 1214-1229.
Bhattacharjee, R., Jayaram, R., Musco, C., Musco, C., & Ray, A. (2025). Improved spectral density estimation via explicit and implicit deflation. In Proceedings of the 2025 Annual ACM-SIAM Symposium on Discrete Algorithms (SODA) (pp. 2693 2754).: SIAM.
Brezinski, C. (1996). Extrapolation algorithms and padé approximations: a historical survey. Applied numerical mathematics, 20(3), 299-318.
Cai, T. T., Liang, T., & Zhou, H. H. (2015). Law of log determinant of sample covariance matrix and optimal estimation of differential entropy for high-dimensional Gaussian distributions. Journal of Multivariate Analysis, 137, 161-172.
Caron, E. & Utard. G. (2002). Parallel out-of-core matrix inversion. In Proceedings 16th International Parallel and Distributed Processing Symposiun (pp. 6-pp).: IEEE.
Chen, T., Trogdon, T., & Ubaru, S. (2021). Analysis of stochastic Lanezos quadrature for spectrum approximation. M. Meila & T. Zhang (Eds.), Proceedings of the 38th International Conference on Machine Learning, volume 139 of Proceedings of Machine Learning Research (pp. 1728- 1739).: PMLR
Chizat, L., Oyallon, E., & Bach, F. (2019). On lazy training differentiable programming. H. Wallach, H. Larochelle, A. Beygelzimer, F. d' d'Alché-Buc, E. Fox, & R. Garnett (Eds.), Advances in Neural Information Processing Systems, volume 32: Curran Associates, Inc
Couillet, R. & Liao, Z. (2022). Random Matrix Methods for Machine Learning. Cambridge : University Press. https:///henp. http:////hnuuaiaooo github io/book/.
Dereziński, M. & Mahoney, M. W. (2021). Determinantal point processes in randomized numerical linear algebra. Notices of the AMS, 68(1), 34-45.
Drineas, P. & Mahoney, M. W. (2016). RandNLA: Randomized numerical linear algebra. Commu- nications of the ACM, 59, 80-90.
Drineas, P., Mahoney, M. W., Muthukrishnan, S. (2008). Relative-error CUR matrix decomposi- tions. SIAM Journal on Matrix Analysis and Applications, 30, 844-881.
Dubbs, A. & Edelman, A. (2015). Infinite random matrix theory, tridiagonal bordered Toeplitz matrices, and the moment problem. Linear Algebra and its Applications, 467, 188-201.
Frommer, A. & Schweitzer, M. (2016). Error bounds and estimates for Krylov subspace approxima- tions of Stieltees matrix functions. BIT Numerical Mathematies, 56(3), 865-892.
11

---

Gardner, J., Pleiss, G., Weinberger, K. Q., Bindel, D., & Wilson, A. G. (2018). GPyTorch: Blackbox matrix-matrix Gaussian process inference with GPU acceleration, Advances in neural information processing systems, 31.
Gholami, A., Yao, Z., Kim, S., Hooper, C., Mahoney, M. W., & Keutzer, K. (2024). Al and memory wall. IEEE Micro, 44(3), 33-39.
Gittens, A., Devarakonda, A., Racah, E., Ringenburg, M., Gerhardt, L., Kottalam, J., Liu, J., Maschhoff, K., Canon, S., Chhugani, J., Sharma, P., Yang. J., Demmel, J., Harrell, J., Krish- namurthy, V., Mahoney, M. W., & Prabhat (2016). Matrix factorizations at scale: A comparison of scientific data analytics in Spark and C+MPI using three case studies. In 2016 IEEE International Conference Big Data (Big Data) (pp. 204-213).
Gittens, A. Mahoney, M. W. (2016). Revisiting the Nyström method for improved large-scale machine learning. Journal of Machine Learning Research, 17(117), 1-65.
Golub, G. H. & Meurant, G. (1997). Matrices, moments and quadrature II; how to compute the norm of the error in iterative methods. BIT, 37(3), 687-705.
Golub, G. H. & Meurant, G. (2010). Matrices, Moments and Quadrature with Applications. Princeton University Press.
Hamilton, R. S. (1982). The inverse function theorem of Nash and Moser. Bulletin (New Series) of the American Mathematical Society, 7(1), 65 222.
He, K., Zhang, X., Ren, S., & Sun, J. (2016). Deep residual learning for image recognition. In Proceedings of the IEEE conference computer vision and pattern recognition (pp. 770-778).
Higham, N. (2022). Seven Sins Numerical Linear Algebra. https://nhi gham com 2022/10// 11/seven- s sins- nuf-numert wwf--umrriaallleeeeeeeee waw---naaa--l Accessed: 2025-01-29.
Hodgkinson, I.., der Heide, C.,, Roosta, F,, & Mahoney, M. W. (2023a). Monotonicity and double descent in uncertainty estimation with Gaussian processes. In Proceedings of the 40th International Conference Machine Learning (pp. 13085-13117).
Hodgkinson, I.,, van der Heide, C., Salomone, R., Roosta, F., & Mahoney, M. (2023b). The inter- polating information criterion for overparameterized models. arXiv preprint ar arXiv:2307.07785wv.
Hodgkinson, L.,, der Heide, C., Salomone, R., Roosta, F,, Mahoney, M. W. (2023c). A PAC-Bayesian perspective on the interpolating information criterion. In Advances in Neural Information Processing Systems, Mathematics of Modern Machine Learning Workshop.
Huang, J. & Landon, B. (2020). Spectral statistics of sparse Erdös-Rém graph Laplacians. Annales de I I'Institut Henri Poincaré, Probabilités et Statistiques, 56(1). 120-154.
Immer, A., Korzepa, M., & Bauer, M. (2021). Improving predictions of Bayesian neural nets via local linearization. In A. Banerjee & K. Fukumizu (Eds.), Proceedings of The 24th Interna- fional Conference on Artificial Intelligence and Statistics, volume 130 Proceedings of Machine Learning Research (pp. 703-711):: PMLR.
Immer, A., Van Der Ouderaa, T. A., Van Der Wilk, M., Ratsch, G., & Schölkopf, B. (2023). Stochastic marginal likelihood gradients using neural tangent kernels. In A. Krause, E. Brunskill, K. Cho, B. Engelhardt, Sabato, Scarlett (Eds.), Proceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research (pp. 14333-14352): PMLR.
Jackson, D. (1912). On approximation by trigonometric sums and polynomials. Transactions of the American Mathematical Society, 13(4), 491-515.
Jacot, A., Gabriel, F,, Hongler, C. (2018). Neural tangent kernel: Convergence and generalization in neural networks. In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, & Garnett (Eds.), Advances in Neural Information Processing Systems, volume 31: Curran Associates, Inc.
12

---

Kesten, H. (1959), Symmetric random walks on groups. Transactions of the American Mathematical Society, 92(2), 336-354.
Kim, S., Park, S., Kim, K.-S., & Yang, E. (2023). Scale-invariant Bayesian neural networks with connectivity tangent kernel. In The Eleventh International Conference on Learning Representa- tions.
King, F. W. (2009). Hilbert Transforms. Encyclopedia of Mathematics and its Applications. Cambridge University Press.
Krizhevsky, A. (2009). Learning multiple layers of features from tiny images. Technical Report TR-2009, University of Toronto.
Lenczewski, R. (2015). Limit distributions of Gaussian block ensembles. Acta Physica Polonica B, 46(9), 1833-1850.
Leskovec, J. & Krevl, A. (2014), SNAP Datasets: Stanford large network dataset collection. stanford
http://snap. http://snapp.taaford.. edu/data.
Longoria, I. S. A. & Mingo, J. A. (2023). Freely independent coin tosses, standard Young tableaux, and the Kesten-McKay law. The American Mathematical Monthly, 130(1), 35-48.
Lowe, D., Kim, M., & Bondesan, R. (2025). Assessing quantum advantage for Gaussian process regression. arXiv preprint arXiv:2505. 22502.
Mahoney, M. W. (2016). Lecture notes on spectral graph methods. arXiv preprint arXiv: 04845.
Maïda, M. (2023). Strong convergence of the spectrum of random permutations and almost- Ramanujan graphs. Séminaire Bourbaki, (pp. 75e).
Marčenko, V. A. & Pastur, L. A. (1967). Distribution of eigenvalues for some sets of random matrices. Mathematics of the USSR-Showe SSooni,, 1(4), 457.
McKay, B. D. (1981). The expected eigenvalue distribution of large regular graph. Linear Algebra and its Applications, 40, 203-216.
Meixner, J. (1934). Orthogonale polynomsysteme mit einer besonderen gestalt der erzeugenden funktion. Journal The London Mathematical Society-second Series, (pp. 6-13).
Meurant, G. (2006). The Lanczos and Conjugate Gradient Algorithms: From Theory to Finite Precision Computations (Software, Environments, and Tools). USA: Society for Industrial and Applied Mathematics.
Mezzadri, F. (2007). How generate random matrices from the classical compact groups. Notices of the American Mathematical Society, 54(5), 592-604.
Mu, X., Zhou, H.-X.,, Chen, K., & Hong, W. (2014). Higher order method of moments with parallel out-of-core LU solver on GPU/CPU platform. IEEE Transactions on Antennas and Propagation, 62(11). 5634-5646.
Murray, R., Demmel, J., Mahoney, M. W., Erichson, N. B., Melnichenko, M., Malik, O. A., Grigori, L., Luszezek, P., Dereziński, M., Lopes, M. F., Liang, T., Luo, H., & Dongarra, J. (2023). Randomized Numerical Linear Algebra: A perspective on the field with an eye to software. arXiv preprint XrXiv22302.11477
Nguyen, H. H. & Vu, V. (2014). Random matrices: Law of the determinant. The Annals of Probability, (pp. 146-167).
Nica, A. (1993). Asymptotically free families of random unitaries in symmetric groups Pacific Journal of Mathematics, 157(2).
Nica, Speicher, R. (1996). On the multiplication of free n-tuples of noncommutative random variables. American Journal of Mathematics, 118(4), -32.
13

---

Novak, R., Sohl-Dickstein, J., & Schoenholz, S. S. (2022). Fast finite width neural tangent kernel. In International Conference on Machine Learning (pp. 17044): PMLR.
Olver, Nadakuditi, R. R. (2012). Numerical computation of convolutions in free probability theory. arXiv preprini arXiv: 220..
Potapczynski, A., Finzi, M., Pleiss, G., & Wilson, A. G. (2023). Cola: Exploiting compositional structure for automatic and efficient numerical linear algebra. In A. Oh, T. Naumann, A. Glober- son, K. Saenko, M. Hardt, & S. Levine (Eds.), Advances in Neural Information Processing Systems, volume 36 (pp. 43894-43917): Curran Associates, Inc.
Rabani, E. & Toledo, S. (2001). Out-of-core SVD and QR decompositions, In PPSC.
Koosta-Khorasee F. & Ascher, U. (2015). Improved bounds sample size for implicit matrix trace estimators. Foundations of Computational Mathematics, 15(5), 1187-1212.
Rozemberezki, B., Allen, C., & Sarkar, R. (2021). Multi-scale attributed node embedding. Journal of Complex Nenworks, 9(2), cnab014.
Saibaba, A. K., Alexanderian, A., & Ipsen, C. (2017). Randomized matrix-free trace and log- determinant estimators. Numerische Mathematik, 137(2), 353-395.
Saitoh, N. & Yoshida, H. (2001). The infinite divisibility and orthogonal polynomials with constant recursion formula in free probability theory. Probab. Math. Statist., 21, 159-170.
Shen, J.. Tang, T., Wang, L.-L. (2011). Spectral Methods: Algorithms, Analysis and Applications. Springer Publishing Company, Incorporated, 1st edition.
Tao, T. (2017).
Inverting the Schur complement, and large-dimensional https://terrytao. wordpress com htom/20177099/ 20//0//1711////// compuer-oommemmetttttomommeeeeee and-- Marre-rrre--a- dimensional selfanddtsetllnn patterns/. Accessed: 2025-05-13
Gelfand-Tsetlin patterns.
inverting-
Trefethen, L. N. (2023). Numerical analytic continuation. Japan Journal of Industrial and Applied
Mathematics, 40(3), 1587-1636.
Ubaru, S., Chen, J., & Saad, Y. (2017). Fast estimation of (a)) via stochastic Lanczos quadrature. SIAM Journal on Matrix Analysis and Applications, 38(4), 1075-1099.
Ubaru, S. & Saad, Y. (2017). Applications of trace estimation techniques. In International Conference on High Performance Computing in Science and Engineering (pp. 19-33).: Springer.
Vakili, S., Bromberg, M., Garcia, R., shan Shiu, D., & Bernacchia, A. (2022). Uniform general- ization bounds for overparameterized neural networks.
Voieulescu, D. (1991). Limit laws for Random matrices and free products. Inventiones mathematicae, 104(1), 201-220.
Wachter, K. W. (1978), The strong limits of random matrix spectra for sample matrices of indepen- dent elements. The Annals of Probability, 6(1), 1-18.
Wigner, E. (1955). Characteristic vectors of bordered matrices with infinite dimensions. Annals of Mathematics, 62(3), 548-564.
Wilson, J., van der Heide, C., Hodgkinson, L., & Roosta, F. (2025). Uncertainty quantification with the empirical neural tangent kernel. arXiv preprint arXiv:2502.02870.
14

---

Appendices
Contents
Nomenclature
15
Benehmark Matrix Ensembles
15
B.1 Overview of Random Matrix Ensembles
15
B.2 Free Decompression of Closed-Form Benchmark Distributions
18
C Considerations for Solving the PDE Using Characteristics
20
C.I Challenge: Crossing the Cut
20
C.2 Holomorphic Continuation of Stieltjes Transform
Methods Estimating Stieltjes Transform
23
D.1 Estimating Stieltjes Transform using Lanczos Method
24
D.2 Estimating Stieltjes Transform using Jacobi Polynomials
25
D.3 Estimating Stieltjes Transform using Chebyshev Polynomials ................................
28
Proofs of Propositions 1 and
30
F Stability Analysis
31
G Background Material for the Neural Tangent Kernel
33
Implementation and Reproducibility Guide
34
Appendix A
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.
Appendix B Benchmark Matrix Ensembles
We evaluate method on suite of benchmark matrix ensembles with known limiting spectral distributions. Appendix B.1 reviews the underlying matrix models and associated laws, while Ap- pendix B.2 presents decompression results these closed-form distributions.
B.1 Overview of Random Matrix Ensembles
We review several canonical spectral laws arising from random matrix and combinatorial models, which serve as ground-truth distributions in our benchmarks. Concretely, we consider the Wigner (1955) semicircle law (free Gaussian), Marčenko & Pastur (1967) (free Poisson), Kesten (1959) and McKay (1981) (free binomial), Wachter (1978) (free Jacobi), and the general free-Meixner family (Saitoh & Yoshida, 2001; Anshelevich, 2008) (free analogs of the classical Meixner (1934) laws).
TTii definition differs by factor of - -1/m from the more common definition 1 v. plx) dy-
x-y
15

---

Table A.I: Common notations used throughout the manuscript.
Context
Symbol
Description
Free Probability
P
Absolutely continuous empirical spectral density, + R Compact support interval of density
I
{A_, >+ ç R
m
Stieltjes transform of density, m(z) JR pl3)
y-z dy, C 1-
H(
Hilbert transform" 2 of density, H[p (a:) p. V. f plv)
dy. -
y-
R
Voiculescu's R-transform, R(2)) w(-2) with ((m(2))
Free
Size of the target (large) matrix
Decompression
n2
Size of the sampled (small) sub-matrix, n2 Decompression scale, t log (.,,
t
Po: mi
Density and Stieltjes transform at
Notation
Upper (lower) half complex plane
R, 3
Real and imaginary part complex variable
V.
Cauchy principal value
In Table B.1, we list, for each law, its absolutely continuous density p(z) supported on = [A_, }. and any point masses. Immediately following, Table B.2 summarizes the simplest random-matrix combinatorial constructions realizing these laws practice.
Table B.1: Spectral laws used in our analytic benchmarks, along with their tree-probability analogs, absolutely continuous density functions p(z), compact support intervals [A_, and any point- mass atoms.
Distribution
Free Corresp.
Abs. Cont. Density p(æ)
Support A
Number of Atoms
2 =
Wigner semicircle
Free Gaussian
None
v(, -x)
Marchenko- Pastur
Free Poisson
(1 vX)
(1 (1
dv 4/d - 1)
Kesten-McKay
Free Binomial
+2V0-1
None
27idd -x²)
(a () 2)(( A.
vb Vola /ola+bb-
Wachter
Free Jacobi
x == 0,1
2nx(1 x²)
ev 16 - (u a)}
Meixner
Free Meixner
a 2vb
At most two
2r((1 c)x² aca 6e
All the above distributions arise arise-up to affine rescaling- from the free Meixner class measures (Anshelevich, 2007), whose aathooonal-polynoonne recursion coefficients stabilize after the first level. Equivalently, the associated Jacobi matrix takes the fordered-Toeplitz form
ag
30
30
31
J
9
31
an a1, B. = B1 (n 2 1).
See Dubbs & Edelman (2015, Table 2) for the specific (a,, B1) giving each law. In particular, setting ag 0, 30 1 yields the standard zero-mean, unit-variance normalization.
16

---

Table B.2: Empirical realizations of the spectral laws listed in Table B.1, The matrix X is assumed have i..d. standard normal entries.
Distribution
Matrix Combinatorial Model
Parameters
Wigner semicircle
Gaussian orthogonal ensemble J (X + x
2/n pn = (llg((n)
Adjacency matrix of Erdös-Rénýí graph G(n,
Marchenko-Pastur
Sample covariance (Wishart) X X E Rexd
X = 5
Kesten-McKay
Haar-orthogonal Hermitian sum Σ (0; + +
2k d 2 2
Projection model d POODDO .ongori & Mingo, 2023)
Adjacency matrix of random d-regular graph
d 2 2
Wachter
Generalized eigenvalues (S), S + 22), S; = x
d1
d2
Arises MANOVA problems
Meixner
Bordered Toeplitz tridiagonal with Jacobi coefficients a B Block-Gaussinn ensembles (Lenczewski, 2015)
a aq
B
Since the Jacobi parameters in the above become constant beyond the second level, the continued- fraction expansion the Stieltjes transform
m(2)
2 00
$
01
01
becomes eventually periodic, exhibiting a repeating tail that we encode by
m(2) 2 c0 T6
where
T= 2 01 88
Eliminating yields the quadratic equation
Q(=)m(=)² - P(=)m(=)++ +1=
where
38
P(z) = 2(= ao)
(z a1);
8
30
(z ao) - -a1) + 16
Q(2) = (= - an)
37
Note that Q22) is at most quadratic (degenerating to linear when 38 -B7)
and P is at most linear.
Solving for (2) yields
m(=) P(z) + -(=))-40(2
(B.1)
2Q(x)
where the branch of the square root is chosen such that m22) defines Herglotz map C+ + C+. The
corresponding Hilbert transform (x) = m(x + i0) inside the support interval I = [A_, X is
thus the rational function
P(x)
Hpp)())
IEI.
(B.2)
2Q(x)'
Table B.3 lists these polynomials P and the corresponding free cumulant R-transforms for every distribution in Table B.I.
17

---

Table B.3; Stieltjes, Hilbert, and R-transforms for the spectral laws listed in Table B. The Stieltjes and Hilbert transforms are constructed from the polynomials P()) and ((2) according to (B.1) 1) and (B.2).
Stieltjes and Hilhert Transforms
Distribation
P(=)
((+)
R T-Transform
Wigner semicircle
4
T
Marchenko-PPane Pastur
1 - -
Az
--Xe
(2 d)z
==
-d +1/1+44
Kesten-McKay
d
d
(a b-0+b-2)-
=(1-2)
-(a +)+++++2- +
2(a 6)2 +
Wachter
a b
a a
2z
Meixner
2)
(1-0)= c)=² 3 + be
c 1-c)
1
(1
aa)²
4b(1
c)2²
2z
(A] Initial Empirical Densiry (a - IK)
(b) Free Decemprrssion
(e) Final Empirrion Demsity (a 32K]
Bardwack Duunty
0.0125
Eupiniel Specteun
Frre Desopeasoon Saltion
Deaelty ((
16K
0.0076
322
0.0025
-00
200
200
400
400
200
200
400
-400
200
200
400
Figure B.1: 1: (a) Estimation of the Wigner semi-cirele law for A = where as matrix with standard normal entries for 1000; (b) approximation of densities for various values of n via free decompression (Algorithm 1); and (c) comparison of our decompressed approximation of n = 32,000 with the known limiting analytic density.
B.2 Free Decompression of Closed-Form Benchmark Distributions
In addition to the Marchenko-Pastur case discussed earlier (Figure 1), free decompression was tested on matrices with Wigner, respectively. shown in Figures B.1 to B.4, pression, we make use of the R-transform. Let Ra (z) denote the R-transform for the Meixner distribution with parameters a, b, cas given in Table B.3. A few calculations show that for any a > it holds (.,zz((z)) =
Wachter, and Meixner limiting spectral distributions, as identify the Meixner distribution under free decom-
(z) where
b(a) ba 2 - c(a))
a(a) aa,
c(a) c+a(1-c)
1-c
Consequently, performing free decompression on a Meixner distribution with parameters a, b, by factor a = n/nn, will yield Meixner distribution with parameters a(n/ng), b(n/n,) and c(n/n,,)
For the Marchenko-Pastur density we also compared the spectral samples of our free decompression estimates with those taken from the analytic density, in order to measure sensitivity the initial sample that was decompressed. We used three measures of accuracy: total variation, Jensen Jensen-Shannon Shannon distance, and degeemeermsssss.. To this end, free decompression was performed on ten different matrices of size n 1000, each with d = 50,000 degrees of freedom, to obtain approximations of the spectral density corresponding to 32,000 dimensional expansions. Eigenvalues were then sampled from the corresponding distributions via quasi-Monte Carlo, and compared to samples from the analytic density with aspect ratio À 50 The average distance total variation was found to be 0.2%, while the average
distance was 1.867% The mean dog-determinant was
18

---

[a] Initial Empitioal Dersity (n 1K)
(b) Pree Decompression
[c) Final Emppáicaal Density (a 32K)
3.0
Uexatry Satinate
Heneimatk Dundty (kaue))
2.5
Regiiica Spectorn
Piie Decsupreekn Solnioo
4X
Denalto ((A)
2.0
16k
1.5
X2X
0.5
-2
2
-2
0
2
0
Figure B.2: Estimation of the Kesten-No Kesten-McKay law from the model A A
(0, E
O (Longoria
& Mingo, 2023, Section 5), with d 26 2 and size 32,000, where O are Haar- Haar-orth matrices generated via Mezzadri (2007). (a) Spectral density from a subsample of size n, 1000. (b) Densities for various via free decompression (Algorithm 1). (c) Comparison of decompressed and exact density at n = 32,000.
(a) Initaal Empirica Deusily (+ = 1K)
(b) Free Decompression
(c) Final Empirical Density = 32K)
10
Frrghtiaal Spectrom
Foe Desompeasos Solutkan
Density ((A)
8
16%
0.2
0.4
0.5
1.0 0.2
0.6 0.8
1.0
0.2
0.4
0.6
1.0
Figure B.3: (a) Estimation of the Wachter distribution with ratios a 80, 50, from eigenvalues of Wishart matrices with d, 80,000 and d2 50,000 degrees of freedom, respectively, and size n 1000; (b) approximation of densities for various values of n via free decompression (Algorithm 1); and (c) comparison of our decompressed approximation of n 32,000 with the
known limiting analytic density for 32 32
50
(8) Initial Empirical Density (u 1K)
(b) Free Decompression
(c) Final Empirical Density (n = 32K)
0.4
Bendurack Dumity (Beart)
Empitical Specieum
Fisc Descapecsus Soluion
Deaait ((A)
141
0.2
12K
0.1
-00
-10
10 20
-20
100
0 10
20
-20
-10
0
10
20
Figure B.4: (a) Estimation of the Meixner distribution with a 0.1, b 4, c 0.6 from 1000 quasi-Monte Carlo samples; (b) approximation of densities for various values of via free decompression (Algorithm 1); and (c) comparison of our decompressed approximation of n 32,000 with the known analytic density given by the Meixner distribution with z 0.566, b x 253, c ~ 0.210.
computed be 13514.88 with a standard deviation of 1.503, versus the baseline of 13353894. This corresponds an absolute relative error of 1.78% with standard error of 0.01%.
19

---

Appendix C Considerations for Solving the PDE Using Characteristics
Our task is to compute mnt, 2) for a grid of points 2 E C+ and times t E r]. given the initial con-
dition m)(0= =) = (z). By Proposition 3, one has the explicit representation mn(1, =) = -mo
where the "initial label" 2g satisfies the characteristic relation
20
(C.1)
mo (2o)
Equivalently, each characteristic is straight line emanating from 20 with direction vector -mp(20)
To evaluate m(t, =) at given 2), one must solve the nonlinear equation (C.1) for the unknown initial point =o (tt, =). Although this root-finding is numerically tractable, it raises two interlock- ing issues. First, the map + @0t, 2) departs the upper half-plane and crosses the real axis through the spectral support 1 := supp(pn). where my is discontinuous. Second, evaluating my on the lower half-plane requires a holomorphic extension of mg across the cut I.
In the remainder of this section, we address these points in two stages. In Appendix C.I, we prove that every trajectory @tt, 2) indeed intersects IR precisely on I. In Appendix C.2, we construct a second-sheet continuation of mg. holomorphic across I, via an additive Ahemann-Hilber (gluing) ansatz.
C.1 Challenge: Crossing the Cut
Consider fixed target point = E Ct. We seek to determine the initial points Zi = +(t, z) that, under the flow of the characteristic equation (C.1), pass through z at time t. Define the curve C + an by =1( ott, =), given implicitly as the solution of (C.1) for fixed This curve traces the locus of initial conditions 2g whose corresponding characteristic trajectories intersect the point 2 at time t.
Unlike the characteristic curves themsellvs-wwine are straight lines this curve (t) is highly nonlinear, beginning at =0(0) z. In the following propositions, we show that this curve descends from its initial point x C+, eventually crossing into the lower half-plane C and in fact intersects the real axis precisely on the support of p.
Proposition C.1 (Characteristic curve exiting Ct). Let Po be an absolutely continuous probability density (with respect Lebesgue measure) supported on a compact interval IR. Assume there
exist a compact subinterval I' C. int I and a constant co > O such that Pa(z) 2 Co for all x E I'.
Let ma(2) denote the Stieltjes transform of Pa: which defines a Herglotz map C+ b st. For each E C+ and t >> 0, define @(t) = olt, =) E C implicitly by
(tt) mo(
(C.2)
Then there exists a finite time 1" 1" (x) > such that 3p((t) 0 0 St<t and So((") = 0.
Proof. (et = e' and define w(2) - 12) so that C+ + -t Fix = ztiy e Ct. Write ((t) = ut) + iv(() and w((((t)) = (t) + iB(t) with 3(t) > 0.
First, we show that when the density Po has atoms, (t) does not asymptote the edges of I as t + co. Taking real parts in (C.2) gives
u(t) + (t))a(t).
If u(t)) + n" E d1 and v(t) + 0t. then by Plemelj- mn + 10) E R is finite (no atom at 2"), hence a(t) + a Rw(u" + 10) -1/mo(u" + i0) * 0. which forces u(t) + ((t) + u" + a(t)a" a contradiction since a(t) - co. Thus the trajectory cannot converge to endpoint while u(t) > 0. This rules out the pathological possibility that the trajectory merely approaches the boundary without ever entering the support. Consequently, we may assume the curve eventually moves inside the interior of I, where the following estimates apply.
Next, we show is bounded below. By the Plemelj boundary values, for u E we have
Po (w)
lim Swou + iv) -
|mg (u + 20)
20

---

Since Po 2 > on I' and (2 + :0) is continuous on I', there exist to > 0 and > 0 such that (u(uuiu 2 c for all I' and 0 < v S vo. Hence, once u(t) E I' with u(t) S vo. we have 3(t) 2
Finally, taking imaginary parts in (C.2) gives
v(t) (()
As + 0o, we have a(t) + oo while {tt) remains bounded below by e whenever u(t) E I' and v(t) - vo. Thus, u(t) + oc, so by continuity there exists first time t" where o(t = c
Proposition C.2 (Crossing on the support). Suppose the hypotheses of Proposition C.I hold, and define the curve >> Φ(t, E C implicitly by (C.2) for fixed a E C Let 1. >> O be the first time such that Sø(t") = 0, whose existence is guaranteed by Proposition C.I. Then the real part of the crossing point lies in the support of the density: Rø(t") E E.
Proof. Let x" = x(t") ) be the real-valued point where curve & intersects the real axis. Suppose, for contradiction, that x & I. Since I is compact, there exists open interval J R containing x². On such an interval J. the Stieltjes transform mz( admits real boundary values. Indeed, by the Plemelj -onell--okkottki formula.
lim mo( + ic) Hx)o](x)) impo(>)) I EER. c-0
(C.3)
and since Po(x:) for all E the imaginary part vanishes and mo(#) is real-valued on J. In particular, my (ax) E RR. Substituting s((") x into (C.2), we find that the right-hand side of the equation is real, so z E R, contradicting the assumption that 2 E C+. Therefore, the assumption & is false, and we conclude x²
Since the trajectory C enters the lower half half-planc planee we must evaluate the initial field mo both across the real line at points on supp(po), and inside C However, by definition, the Stieltjes transform
pPo(x) m²(2) - - dz,
is only defined for ((2) > O, and exhibits jump discontinuity across the real axis precisely on supp(Po). This creates an analytic obstruction: evaluate (o(t, 2)) for arbitrary we require a holomorphic continuation of my to the lower half-plane that is smooth across the cut I. This continuation constructed in the next section.
C.2 Holomorphic Continuation of Stieltjes Transform
Let mmm(2 (z) denote the Stieltjes transform corresponding to density p supported compact interval I C R (such as the field m²(2) and the density Po in our case). By definition, this function Herglotz map C+ + Ct holomorphic on C+ and admits boundary values on R via the Sokhotski-Plemeii theorem:
lim m m (x + ic) Hxp)(x) imp(x))
-
natural extension of m² to CT is the Schwarz reflection
m (2) m mz)) 3(2) < 0.
This extension holomorphic on R but has discontinuity across J. Specifically, the discontinuity arises from jump in the imaginary part:
lim 3(m (x (x +
3(m* (x ii)) = 2Tp((x).
Thus, 3(m is discontinuous on I, while it vanishes and remains continuous on R I. We also note that the real part, R(m H[pl (I:), is continuous across the entire real axis.
Our goal is construct second analytic continuation of m denoted m that is holomorphic across the interior of the cut I. That is, m agrees with m on C+ and continues analytically to C
21

---

through I, at the expense of developing discontinuity across R \ I. This setup forms an additive Riemann-Hilbenn problem: find two branches m and m such that
lim m + (x + ie) lim m (x ie), a NEI. -
e-0
We propose the following structure:
m (=), 1-m (z) G(2),
3(z) 0, 3(2) < 0,
(2)
(C.4)
where G(z) is glue function chosen to cancel the jump across /. In particular, we require that
3(G((z + 10)) =
x E
(C.5)
R(G(x + :0) = (x),
x EEI.
The first condition in the above ensures the continuity of the imaginary part on the whole of R::
lim (mm (æ ie)) 3(m (æ ie)) πT(x)) πT(z)) 0,
c-
and the second condition guarantees continuity of the real part across
lim R(m (xr + ie)) R(m (x ie)) O, e-
while allowing discontinuitics on \ I, which are harmless for our purposes. In this way, m+ and m form two complementary Riemann sheets: m holomorphic C+ UR (7) but discontinuous across I, while m is holomorphic in C U I, with branch cut on IR \
To characterize G. note that G(z) = m (z)+ + m (z). and it must be analytic in C On the real axis inside the support, this glue function satisfies RG(x) 2H [x]),, and thus inherits key structural properties from the Hilbert transform. In order to design an efficient functional form for G, we first study the behavior of H[p) itself in Proposition C.3, which guides the minimal rational approximation ansatz in the next step.
Proposition C.3 (Properties of the Hilbert transform). Let be bounded, piecewise- density supported on a connected compact interval I [}_, X4] C TR. Define the Hilbert transform
Hipp(x) = P.v. (/m
p(() t-x dt,
where p.v. denotes the Cauchy principal value. Then H(pp and its derivative with respect to satisfy:
1. For x E (-00, H(pp(z) > 0 and (zp) > 0.
2. For x E [A= the equation Hzp](z) = 0 has an odd (hence at least one) number of solutions x" E (A_,44).
3. For x E (14,00) (x) < 0 and H(pp (x) 0.
4. As - too, (x) -1/x.
Proof. Step 1. Outside the support. If & [A_, A41. differentiation under the integral sign valid:
H(p)(()) = (rr d+
p(t)
dt >
(t (t-
Hence, H[p] is strictly increasing on (-o0, (A4, A direct sign check of the integrand confirms that Hxp((x > 0 for æ A_ and H(pl(z) 0 for \4. establishing items (1) and (3). Furthermore, the asymptotic statement in item (4) follows directly from (King, 2009, Eq. 4.111). which gives H[p[(a) ~ - pdr) /xr as |ar + oo; since dd this is -1/a.
Step Existence of a zero. From steps (1), we have H(p
x X_, and H[p) < as X4. Since H[p) is -aattttuuu-gganteeeede by by the boundedness of p the intermediate value
22

---

theorem ensures that there exists at least one root xx' E (4_ such that (xpp(x") = 0. Moreover, every simple zero flips the sign of Hi[pl, so starting with positive sign on the left and finishing with negative sign on the right forces an odd number of sign-changes and hence an odd number of interior zeros.
The properties in Proposition C.3 provide clear picture of the behavior of the Hilbert transform of a probability density with compact support; see Figure C.I for an example.
Hzp)(z)
Figure C.I: A schematic example of the behavior of the Hilbert transform as prescribed by Proposi- tion C.3.
We now seek practical approximation for the real part of the glue function G(ær) on the support interval J. Motivated by the fact that many classical spectral laws-such as the free Meixner family of distributions in Appendix B.1- admit closed-form Hilbert transforms of rational form, we adopt a Padé-type ansatz:
P(x)
G(x)
(C.6)
Q(z)*
with P and Q polynomials of minimal degree. In the classical examples of Appendix B.I, this rational structure is not merely an approximation but exact, shown in (B.2), with P at most linear and Q at most quadratic (see Table B.3).
Even in more general settings-wheet where the density p is empirical lacks closed-form expression- the Padé approximant remains attractive as provides a compact, flexible representation with inter- pretable parameters. Once commit this functional form, the Hilbert transform to constrain the numerator and denominator more effectively. Proposition C.3 provides key analytic features of [[p that guide this design. Namely, P should have at least zero inside 1, while Q should have zeros in C or within I, though it may have zeros in IR \ T. Minroring these features, we parametrize the Padé approximant as
leverage structural properties of
P(z) d G(z) Q(2)
Σ
Ti 2 2 aj
which is a Padé approximant of degree [q + 1/g]. or [q/9] if or 4 1/g] if = d We constrain the parameters to ensure there at least one zero in / and no pole a, inside 1.
Appendix D
Methods of Estimating Stieltjes Transform
Given probability density p(x) supported on an interval / C IR, its Stieltjes transform is defined by
p(x)
m(z)
dx, 2EC\T.
2 x
The density p can be recovered from the boundary behavior of m(()) via the inversion formula:
1 lim 3m(z + ie). क c-
p(x)
Accurately estimating m(x) is therefore essential for recovering the underlying spectral density. Given Hermitian matrix A E =x" the Stieltjes transform of its empirical spectral distribution is
23

---

defined by
1 m(=) - Lr (A - =1) -I)--
=ECVL
A direct way to compute m(=) is
m(z) 1
where A are the eigenvalues of A. However, because the characteristic curves must pass through the roots of m (Proposition C.2), applying free decompression to this empirical version of m will remove its poles. necessary consider approximations of m that correspond Stieltjes transforms for smooth densities.
In this section, describe three numerical strategies for estimating the Stieltjes transform of Hermitian matrix. The most widely used approach in the literature is based on the Lanczos algorithm (Appendix D.1), which constructs continued fraction approximation to the resolvent. However, we find that this method often suffers from slow convergence and numerical instability, due in part its reliance on random initial vector and the large number of Lanezos steps required achieve acceptable accuracy. As such, it is often ill-suited for tasks that demand high-precision estimates of m(z).
To overcome these limitations, we present two alternative methods. The first, described in Ap- pendix D.2, expands the spectral density in a Jacobi polynomial basis, offering flexibilityy in modeling various edge behaviors. The second, described in Appendix D.3, specializes this approach to Cheby- shev polynomials, for which derive closed-form expression for the Stieltjes transform of each basis element. This enables direct and stable evaluation via power series. Together, these methods offer more accurate and robust alternatives to Lanczos-bats approaches, especially in regimes where high-resolution spectral information is required.
D.1 Estimating Stieltjes Transform using Lanczos Method
The key idea the Lanczos method is to approximate m(z) by constructing a rational function specifically, continued fraction- reflects the spectral properties of A. This is achieved via a three-term recurrence relation satisfied by orthogonal polynomials associated with the spectral measure (Golub & Meurant, 2010; Meurant, 2006). The recurrence is obtained from the Lanczos algorithm, which builds a tridiagonal (Jacobi) matrix capturing the action A on low-dimensional Krylov subspace. In brief, the procedure is follows.
Given a normalized starting vector vy R", ||o|| = 1, the Lanczos algorithm constructs an orthonormal basis (Vo, UI...., 1 for the Krylov subspace
(A, vn) = span {ro: Au A AZ vo.... A A
and produces the symmetric triadiagonal Jacobi matrix
0gg 8
31
01
B2
T
Bp-1
ap-1
which satisfies
AV, V,
where V, : (Vo, up-1] nxp, This identity captures the three-term recurrence relation satisfied by the orthogonal polynomials associated with the spectral measure of A. The Stieltjes transform m²)) can then be approximated by the first element the resolvent of Tp:
m,(x) = =1) en.
24

---

where e = (1,0, ...... 0) E I" denotes the first standard basis vector. The above expression is equivalent to the continued fraction
mp (2)
0g
C1
B7
ap
This approximation is the [0/p] Padé approximant of ml(2), and typically exhibits convergence with increasing p, especially when m=)) is analytic away from the real axis (Golub Meurant, 1997).
However, the effectiveness of this approximation in practice depends on several implementation details, which we briefly summarize below.
The choice of the starting vector is not critical, as any unit-norm vector suffices. In practice, randomly chosen unit vector often leads convergence, while theory, the approximation m2))) is independent of this choice when p However, for larger values of numerical stability may require full or partial reorthogonalization during the Lanezos process to maintain orthogonality of the Krylov basis vectors (Meurant, 2006).
typical stopping criterion is based monitoring the difference between successive approximants and terminating when
(=) (1)p-(()) <,,
for a fixed test point E CT: see (Frommer & Schweitzer, 2016; Ubaru et al., 2017) for detailed analyses of convergence rates and practical recommendations.
The continued fraction (2) can be evaluated efficiently either through recursive computation or, equivalently, by solving the tridiagonal system (T, =I)æ e₁ The overall computational cost (pm²) for dense matrices or (nz()) for sparse matrices (where nnz( (A) is the number of non-zero elements A), plus Opp") for reorthogonalizatio.. If the algorithm run completion with = n and numerical orthogonality is preserved, the approximation yields the exact spectrum of A, and thus fully recovers m(2).
D.2 Estimating Stieltjes Transform using Jacobi Polynomials
Jacobi polynomials provide a highly flexible basis for approximating spectral densities, particu- larly when the distribution exhibits distinctive edge behavior-suce as square-root singularities or vanishing derivatives at the boundaries of the support. This family of orthogonal polynomials can represent a wide range of distributions encountered in random matrix theory (e.g., distributions in Appendix B. and more general beta-like distributions.
Let p(r) be probability density with support the interval [A_, Through affine change of variables, we map the domain to the canonical interval [-1, 1 via
2x - (x+ + +
(x)
1 >
The density is then expanded in the Jacobi polynomial basis:
p(x) Sh w (o.0) (((x) (1(x))
(D.1)
where p(o..) is the Jacobi polynomial of degree l and parameters a, > -1, orthogonal with respect the weight
w (a.)) (1 x) (1
(1-
0, -1
This expansion is especially effective when p(x) ~ (x) near the endpoints, aligning the basis with the singularity structure of p. and accelerating the convergence.
25

---

Projection. Given set of eigenvalues {A, we first construct smoothed density estimate p(z). e.g., via Gaussian beta kernel density estimation. Mapping to the variable t, we approximate the expansion coefficients using the orthogonality of Jacobi polynomials:
Σj Etttt))P ( (t,)
[
(1)
(D.2)
where { grid points closed-form expression
[-1, 1] and || |? denotes the squared norm, which admits
2
T(k + + 1) I(A +
|
2k + a + +
[(A: + 1) T(k + a +
In the above, I is the Gamma function.
Stieltjes transform in Jacobi basis. Once the coefficients, t/e are obtained, the Stieltjes transform
m(z) //^
p(xx
dx,
x
can be approximated via the expansion (D.1). First, by mapping and
1
x(t) 2 (A+ + A-) ) +
(A+ - x_)r,
t=|-1,1]
2
1
<(u) (A+ + }
(A+ \_),
2
we arrive at the transformed expression
K
w (a,3) (t) dt, u ue we 1-1,1.
m(u) ~
th
X
1-u
k=
The inner integral in the above defines the so-called Jacobi function of the second kind,
Q
(w)
["] w (a.9) (t)
(a.0)
(t)
dt,
t t-
and the Stieltjes transform becomes
m(z(u)) ~
2 x X_
Σ
the Qa
(0. (u))
k=
Gauss-Jacobi quadrature. In practice, we do not evaluate (f (a. directly via numerical integra- tion. Instead, we compute each using dedicated Gauss-Jacobs quadrature rule (Shen al., 2011, Section 3.2.2) tailored to the weight w (a,9)
pla (t,)
(a.9)
where (4,, w are the quadrature nodes and weights for the given weight function. The number of nodes n, depends the degree k typically chosen ny. = max(k 1, Evaluating each mode Qla.9) (u) independently using its own quadrature is critical for numerical stability, especially for large k: This separation of summation and integration avoids artifacts that arise when attempting cendratuur-integran the full series at once.
Choosing parameters. Jacobi polynomials provide flexible basis for modeling spectral densities with various edge behaviors, particularly those of the form p(x)) (A4 - xx) (x - >_) Different choices of and 8 yield classical orthogonal polynomials:
26

---

Damping Gibbe Occillations
Rmpirial draaity
Thocehhoal BarcheskoPPatuu dusily
1.5
Appozimatio wilhout Jukee daming
Appossimation Jadosn duming
1.0
0.5
0.0
15
2.
3.0
Figure D.1: Effect of Jackson damping on approximating the Marchenko-Paun Pastur density from the Gram matrix XXXT, where X is d random matrix with 1000 and n/d 0.9. Parameters were intentionally selected to amplify Gibbs oscillations for illustrative purposes; typical cases exhibit much milder artifacts.
a recovers Legendre polynomials.
a gives Chebyshev polynomials of the second kind, well-suited for densities with square-root vanishing edges, such as p(x) ~ v(>+ - x)(x X_), which are common in free probability models (see Table B.1).
a
yields Chebyshev polynomials of the first kind, appropriate for densities with inverse square-rood root singularities at the edges.
This parameterization aligns the orthogonal basis with the singular structure of p(z), often resulting in exponential convergence with respect the number of terms K. In practice, find that the method is not highly sensitive to precise values of a and B, as long as their qualitative effect on edge behavior appropriate.
Gibbs oscillations. Expansions in orthogonal polynomials exhibit Gibbs oscillations, espe- cially for large truncation orders K. These oscillations can be suppressed by applying damping to the coefficients via filter Is. resulting in the smoothed expansion:
K
g Pl ((a..)) (x))
p(x)
k=0
Several damping schemes exist, including Jackson (1912), Lanczos, Fejér, and Parzen filters. In this work, we use Jackson damping, defined as
(K - k + cos
Tk K+I
sim
Tk
cot
K+I
K-I
0, K.
K+1
Figure D.1 shows how Jackson damping improves the stability of the approximation.
Regularization. modifying (D.2) as
To further stabilize the coefficients v/. we introduce Tikhonov regularization by
[ p(t)
2k + +1| ||
with the penalty term 7/ chosen as
k:
The
K
where 7 > 0 is small regularization parameter. This form penalizes higher-order coefficients more strongly, resulting smoother approximations without significantly biasing the shape or integral of the estimated density.
27

---

Preserving positivity and mass. Finally, to enforce positivity of p(z) over its support and ensure p(x)) dr 1, we perform secondary constrained optimization over 4... initialized at the projected values, subject to constraints enforcing unit mass and pointwise non-negativity over fine grid.
D.3 Estimating Stieltjes Transform using Chebyshev Polynomials
Constructing the principal branch of the Stieltjes transform numerically is straightforward: one can estimate m on the real axis by letting m(r + =0+) = Hzp](z) + +imp(z) and estimating the Hilbert transform using FFT techniques. Then, m can be extended to C* by convolution with the Poisson kernel. The lower-half plane is obtained by the Schwarz reflection principle. However, constructing the secondary branch of the Stieltjes transform numerically is ordinarily very challenging, as it is tantamount to numerical analytic continuation.
Our primary mechanism for estimating the Stieltjes transform of density with compact support is the following lemma. Here, we let j(z) = (2 2 denote the Joukowski transform, and let
((z) =
denote the inverse Joukowski transform, where take the branch of the square root to have positive imaginary part for C+. The Chebyshev polynomials of the second kind Ug (xx) defined for x E [-1, 1] by the recurrence relation
(++1(0) - (x)
from U2(2)) = and (1) 2x.
Lemma D.1. The Stieltjes transform of U2(z) x² on [-1,1] is given by Tz(2) for any E (0) UN.
Proof. Let m(2) (2) Since m C+ - C+ is analytic, by the Herglotz-Novvaninna
representation theorem, there exist constants C, D > 0 and a Borel measure fa such that
m(2) = C Dz + / A-z du(A).
1
Since J(in) - 0 as + oo, C D = 0, and so m is the Stieltjes transform of measure If
cos 0 for 0 E r]. since j(e cos 0, j(e
and J(x) =e ee Since Uj (cos 0) sin 0
sin((k: + 1)0), it follows that
U2 (x x² 3(c ((++1))))))) - - ((
Since J is analytic on C+ this, in turn, implies that
du(x) dr 6-+0
1 3[m([ + +18)] Uu(())
lim
and the result follows.
The significance of Lemma D.1 is that we can approximate a spectral density as a Chebyshev series, and immediately read off the corresponding Stieltjes transform. To see this, suppose that spectral density p on [-1, 1] has the Chebyshev expansion
p(x) = V x
(x)).
By Fubini's Theorem,
E Ex [x Σ 4/ 1 U
20
π x² dz
2 4/k
28

---

and so 1ie can be estimated by Un., 2 Ei= U2(A,). The corresponding Stieltjes transform is then m=) =z)) where
A(2)
11 M+1
k=0
Noting that [](z)
for E Ct [J(z)] S and so J(2) will lie in the domain of convergence of the power series for A. However, this series is unlikely to converge for 2 E C and is therefore insufficient for estimating the secondary branch. The solution is to analytically continue outside the domain of convergence its power series by appealing to Padé approximants. Computing Padé approximants algebraically is expensive task, so appeal to Wynn's c-algorithm to rapidly evaluate the Padé approximant for A, given estimates of coefficients 4/x (Brezinski, 1996). For sequence of partial sums (2) = Ek= thaz
s+1 2+1 = 0, 0,1,..., the
(n) e-algorithm constructs tableau {Ek (z)]n,k by
22) = 0, co (m) (z) =
(((+t)) 1 (x)())= = 66-1 (x) + (Ek (2) -6
(m) (=))
22 (z) provides the Padé approximant for A with numerator and denominator degrees of k + n and k, respectively (Baker Jr. & Graves-Morris, 1996, Eq. (4.10)) In particular, Nuttal's Theorem (Baker Jr. & the analytic continuation of A22) regardless of whether the partial sums A, converge.
The coefficient (m)
1996, Theorem 6.5.3) shows that 22n (z) converges in measure
Algorithm 2: Padé Chebyshhey Free Decompression
Input
:Hermitian matrix E R=*";
Maximum Chebyshev polynomial order K;
Evaluation set C IR;
Subsampling parameter n. (with 1 1 n, n));
Perturbation where < ô <
:utput:Estimann spectral density values {p(x))xex-
// Step Generate eigenvalues from random submatrix
1 Extract a random submatrix À from A and compute its eigenvalues (Ai)==1-
2 x_ + min A A 8n
// Estimate of the left endpoint of the support Estimate of the right endpoint of the support
3 A max()
8n-1
// Step 2: Estimate Chebyshev coefficients
4 for k: = to K do
2n,k T( Jn E U& (A:))
I/ Step 3: Compute the Padé approximant
6 Using Wynn's s algorithm, obtain the Padé approximant A(2) for the power series
Ek Vink 1+11
// Step Form the Stieltjes transform
7 m(z) + π A(J(M..))))
// Step Solve characteristic curves
8 foreach x E X do
Using Newton's method, find Z E C such that
1
+ i0
2x
n
m
2g
// Step 6: Return the estimated spectral density n,
10 return { 1 mm2z) n aaxx
29

---

To extend this procedure to an arbitrary support interval [A_,24h we can make use of conformal map
MA .4 (2) 2z - ++)
1
x A_
In this case, the spectral density is given by
((x) 2v (2+ - x.)
where the coefficients of the estimates and the Stieltjes transform becomes
Lnnk Tn(A+ ) (n();),
m(z) (M(
Appendix E Proofs of Propositions 2 and 3
Proof of Proposition 2. By construction. R(t, =22 RR(0,2) 2). In terms of the inverse Stieltjes transform, this implies
w(t,e-'2) (2)6-2))(1-ee3) 1 (0,
Differentiating both sides in / gives
dw (1, 21=) zE
-w Cw
e² =0
(1, - 2)
(E.1)
0z
Since m(t,w(t, 2)) = 2, dz 82
8m aw and 3m 0t
+
8m Bw Uz si
= 0. Applying these relations to (E.1) gives
è' dm
(1, w(, 6 22)) ze
(t, ((, e '2))
Ot
2 Oz
Now evaluate this expression at = where = = = so that wll, (z) Then
dm (t, 2) m(t, 2) + m(t Oz
dm (t, 2) =0
Öt
Rearranging this expression implies the result.
Proof of Proposition 3. Let + (t(7), x()) be a characteristic curve of the quasilinear PDE (1).
By the chain rule,
2m(t(r),,(T)) @m, 2) dt(r) Ōt dr
ami(,, 2) dz(T)
+
(E.2)
dr
dz dr
Matching coefficients in (E.2) with the PDE (1) yields the system of ODEs
dt(T) 1,
t(0) = 0,
(E.3a)
dr
dz(T)
= -m(r)
z(0) = 2u-
(E.3b)
dr
dm(r)
dr -m(r), m0) = mo (2
(E.3c)
Solving (E.3) straightforward, yielding (2).
By standard ODE existence and uniqueness theorems, for each 2g there is unique characteristic curve + =(T), m(()) near - = Since dt dr = and dz # O if mp(2)) * 0, we can dr (locally) invert t and (20) + 2. This yields unique C² function m(t, 2) solving (1) with m(0,2 2) mz).).
30

---

Finally, the condition m²(2) # O for all 2 E Ω ensures that m(T) * O in neighborhood of r = 0, so the term m in (1) is well-defined. Uniqueness follows from the uniqueness of solutions to (E.3) with given initial data. We see that Proposition 3 gives explicit solutions to (1) in terms of the initial data and desired time t.
It remains to show the existence of ø. Define ( 2, 2n) == 2g mz0) (1-e'') 2. If 1f * 0, by the implicit function theorem (Hamilton, 1982), there exists an inverse function 2g = ot, solving (t, 2, ø(t, =)) = O, which allows the explicit solution of (1) by m(t, 2) mo( =]}e This completes the proof.
Appendix F Stability Analysis
To analyze the stability of m(t,2) =) over time i, we perturb the PDE (1) by replacing with + 8m, and derive PDE governing the evolution of the perturbation field 8m(t, z). We then quantify the impact of this perturbation using an appropriate functional norm.
A natural choice for this purpose is the 2,2 -.--oomm along horizontal slices in the upper half-plane, defined by
[[m[tt. + iy)) = [f Im(t, x + iy)l dæ,
where ry == {== x iy R} is the horizontal line at height y O. This norm is particularly suitable for our analysis since we are ultimately interested in evaluating mtt, x + i6) slightly above the real axis, in order to recover the spectral density from the Sticltjes transform.
Moreover, this space includes the boundary values of Stieltjes transforms m(t, x + i0*) whenever p(t,,) E T.2 (R), since can show that
||m(t, =0 =0==00==) M22,,,- (1 (1 m lel(,, )l22
We now establish an upper bound on the growth of this norm over time, culminating in a Grönwall- type inequality. This is formalized in the following proposition.
Proposition F.1. Suppose p(l,-) E I. (B), and let m((, 2) be the Stielljes transform defined on the domain C+ Let dm be small functional (Gâteaux) perturbation of m (so [|ōm]] < 1) satisfying the Bernstein inequality
[m,,mm(t, (r,,) S M]]mm(t,))lee20 2))l) 3ET
(F.I)
which means on each horizontal slice, 8m(tt in) is M-band -imited (Paley (Paley-Wiener space), ensur- ing sufficient decay of the perturbations. Then, for every >> O, it holds that
d
16m(()))2((rr)
(-2 (t)) [16m(t)) (T)))
(F22
dt
where
R m,, m(1, 2)²
1
Cy(1) = 3 sup
+ 2M sup =ET
(F.3)
ZET
ml, 2)
Proof. Let m + ôm, where 5m E (T 2 c+) (the usual Hardy space) is small pertur-
bation with the boundary condition 8m(t, 2) at |z| + xc. Plugging in into the PDE (1) and subtracting the unperturbed PDE gives the perturbation equation
ddm -mm
(mm + dm) 8m
m + &m
dz
m 0=
Expanding the nonlinear term using first-order approximation and rearranging, obtain
06m 1 38m dz 8m
dm 0m
+ O((mm(²).
Ot
m
m Oz
The linearized equation above is essentially a transport-type PDE, with the first, second, and third terms on the right-hand side representing the transport term, damping term, and perturbation term, respectively.
31

---

To obtain the energy, we multiply both sides of the above by dm (the complex conjugate of dm), integrate over r, y and consider its real part, yielding
R / dbm dm dx = 32 [ JR Ot
1 dôm 8m dz
/
(1+
dm
R
|ôm] dz.
(F.4)
dz
m 0z
R
The term on the left-hand side of (F.4), after taking y>o can be expressed as
| 0dm
R dm dæ
| (T,,
(F.5)
R Ot
dt
Since dm is holomorphic, we have 0,Im = d, on every stripe yy Writing m -! + with real valued a, we have
1
0am m
32 (0,8m) dm
R
5mm² - 63 ((0,6m) dm) + 2z
2
Integrating over E R. the third term the right hand side of the above vanishes. As for the other terms, we apply Sauchy-Schwand and the Bernstein bound (F.1) to obtain
R {
06m m 0=
(
3.m 2
1 |8m]] dx. m
8m S
sup xET
+ MM sup sET,
(F.6)
Also, the second term the right-had side of (F.4) satisfies
1 3m RR / ( 1+ m 0z |ōm| dx s -
0_m² m
|ôml dx + sup sel
|5ml dr.
(F.7)
Putting all together by substituting (F.5), (F.6), and (F.7) into (F.4), we obtain
6m 1
(²(r,)) <
1 +
sup
R
m
+ M sup zEr,
3 1 m
1)
18md da
dt
2 zET,
m
(-1 (-1++c,(4) C2)) (1mml P
This completes the proof.
Remark E.I To complete Proposition F.1, verify that the constant Ct(() is finite for every
0 < y < co. Recall the representation the Stieltjes transform,
p(t,
mlt, 2) = - E - dg,
(F8)
for 2 x + C+ with y 0. Then, since p(t,xx) is supported compact interval of diameter
D = diam(supp(p)). we have
[m(4, =)] 2m((,=) = / E(((
de 2
y
pll,E) = D2 + y
(E x)² 1 D² + 10
J
It follows that
D2 +
sup mmm(t 2) sup |m(t, (0 1 =E,,
(F99
=E,,
y
On the other hand, from the derivative of (F.8):
pll,5) dE,
0_m(t, _ (E - 2)²
SO
10m((1,))] 1 p(1,E)
1
dg
< y
(E - =²²² + y²
32

---

Combining the above with (F.9), we obtain
sup R 2gm m(t, z)²
ag (L,( 2) m(t, =)²
D2 + y² (D² + )
sup =EI,
(F.10)
=ef; 2
y/
y/
y²
Substituting the uniform bounds (F.9 and (F.10) into (F.3) confirms that Cy(1) < for all 0 <ソ< 00.
We note that this bound is not defined at y 0, which is expected: y - 0*, the Stieltjes transform approaches the real axis where it develops discontinuities, and quantities like 2_m may diverge. However, in both our theory and implementation, we evaluate m(t, z) slightly above the real axis (e.g., at 2 = x + for fixed < 6 < 1), where all terms remain finite. In practice, the observed error propagation is significantly milder than this worst-case upper bound suggests.
A
Proposition F.1 then implies that the perturbation 8m satisfies an exponential growth bound in time:
() S 11ōmoll22(r,))"
with exponent v = -2+ Cy < 0o for any fixed y > O. But since e = n translates to polynomial growth in n,
the exponential bound n,
[16m(())((rr) S {16mmlll::(r,) 11g
Appendix G Background Material for the Neural Tangent Kernel
The neural tangent kernel (NTK) is an important object in machine learning (Jacot et al., 2018) 8), and has become a prominent theoretical and practical tool for studying the behavior of neural networks both during training and inference (Novak et al., 2022). The Gram matrix associated to the NTK has been used in lazy training (Chizat et al., 2019), and as tool for uncertainty quantification and estimation (Immer et al., 2021; Wilson et al., 2025), exploiting the deep connections between neural networks and Gaussain processes. More recently, it has been shown to arise naturally as a measure of model quality in approximations of the marginal likelihood (Hodgkinson et al., 2023b; Immer et al., 2023), as a tool for quantification of model complexity (Vakili et al., 2022), and in estimation of generalization error in PAC- Bayes bounds (Hodgkinson et al., 2023c; Kim et al., 2023).
While it was first derived in the context of neural networks, the quantity well-defined for a broader class of functions. For clarity will consider the setting where the underlying function is C². However, this can be relaxed in practice to include, as an important example, networks with ReLU activations. For such a function fe X ae parameterized by 0 R², we define the (empirical) NTK to be
xo( (f(())) Jo (J(
(G.1)
Here, Jo(fe(z)) E Rdxp is the Jacobian of fe with respect to its (flattened vector of) parameters 0, evaluated at the data point æ. Since (g x²) Idd ddxd evaluating the NTK data points produces a fourth-order tensor of shape (n, n, d, d). For computational convenience, this tensor is typically reshaped into woodimensionaal block matrix of size nd x nd, where each (ilooc contains the d x d matrix k(x4,
The fact that the NTK arises in variety of contexts makes it natural object of study. For well- trained models, the Gram matrix is nearly singular, making accurate characterization of its spectral properties crucial for estimating quantities used in downstream tasks (Ameli et al., 2025). This also motivates considering the spectrum on a logarithmic scale. In Section 5, we analyze the empirical NTK of ResNe150 model trained on the 10 dataset, 10-class classification task. In this case, 10% of the eigenvalues are extremely small (~ 10 10-14, reflecting near-singular behavior likely tied to the model's ''s class structure.
To examine the scaling behavior of the non-singular bulk spectrum, we project the NTK matrix onto its non-null eigenspace, removing the low-rank null component. We then randomly sample orthogo- nal eigenvectors within this reduced space and reconstruct corresponding Hermitian matrix, which is subsequently permuted to generate random subsets for use in our free decompression procedure.
33

---

(a) Principal Bracch on Cf and C-
[b) Principal Branch on C*2 Secondary Branch C-
2.0
1.0
1.0
0.5
0.5
[m(a]
-0.5
-055
-1.0
-1.0
-1.5%
20
-20
10
Re(=)
Re(s)
Figure G.I: Analytic continuation of the Stieltjes transform of density estimate of the spectral distribution of the log-NTK matrix described in Section 5. (a) The principal branch contains branch cut along the support of the spectral density, while (b) the secondary branch is continuous this region.
The Stieltjes transform of the approximated spectrum of this matrix is depicted in Figure G.1, which shows its analytic continuation using the principal branch (left), well as the secondary branch (right).
Appendix H Implementation and Reproducibility Guide
We developed Python package, freealg," which implements our algorithms and enables repro- duction of all numerical results in this paper. A minimal example using the freealg. FreeForm class is shown in Listing 1. This example generates the results in Figures and 2, illustrating the Marchenko-Pass Pastur distribution with parameter 3/ starting from matrix of size n, 1000 and decompressing to matrix of size 32,000.
Hyperparameters used for the numerical experiments appearing in this document be found in Table H. 1 (notation summarized in the caption). Notebook files that reproduce the figures appearing in this work can be found in the codebase. Further details function arguments and class parameters are available in the package documentation.
Table H.l: Hyperparameters used the numerical examples. Notation: h smoothing kernel band-
width; (a, () Jacobi shape parameters; K Jacobi polynomial order; regularization strength; Padé (p, 9) numerator/denoominato degrees; 8 imaginary shift for Stieltjes evaluation
Smoothing
Jacobi Polynomials
Stieltjes Transform
Example
Reference
Kemel
(a, ())
K
?
Damping
Padé (p. g)
8
Marchenko-Pare Pastur
Figure
Beta
3 10
(2.
50
0
Jackson
(1,1)
10
Wigner semicircle
Figure B.1
Beta
10-
(1)1)
50
0
Jackson
(1,1)
10
Kesten-McKay
Figure B.2
Beta
10
Jackson
(1,2)
10
Wachter
Figure B.3
Beta
10
(1.
Jackson
(1,2)
Meixner
Figure B.4
Beta
10
(1.1)
Jackson
(7,8)
Page-Page
Figure (a)
Beta
10
(0,0)
100
Jackson
(4,3)
NTK
Figure (b)
Beta
10
(2, 2)
50
10
Jackson
(3. 2) 10
ffreealg is available for installation from PyP1 http:p////pppoorrgrrrreeet/teeee the doc-
umentation can be found at https: tttpps://mmell h/tp///mmll..////////eee/eeeeeee wwwthib.o//freaall.......... and the source code is available at https:/ /ttp://gttttt com/ameli//reeal.......a
34

---

Listing 1: A minimal usage example of the freealg package.
Instalt freealg weth "pip install freealg"
import freealg as fa
Create object for the Carchenko-Pastsw distribution with the parameter À mp fa.distributions MarchenkoPastur (1/50)
Generate matrio of size nt, 1000 corresponding to this distribution A (pppatrix(ssze-11000
Create free- object for the matrin within the support 1 = [A_, if Fa.FreeForm(A, support (mm_1am_m, mp.lam_p))
Fit the distribution using Jacobs polynomials of degree K 20, with Also fet the glue function esa Pade of degree [p/q] weth p p= 9= 15 pai if.ditt(mettodd=ajaoo"o,, K-20, alpha-0.5, beta-0.5, reg-0.0, damp- jackson', pade_p-1, pade_q-1, optimizer="le', plot-True)
Es Estimate the empirical spectral denséty p(x), similar Figure 1(a)
rho ff. density(ppoo==Tru))
20
Es imate the Hilbert transform (x)((x)
hilb TI.hilbert(ppoteee (plotTrree)
Es Estimate the Stéeltfes transform m(z) (doth branches m and m ), simslar to Figure (aa
m1, (f.atielljes((lootrruu)
Decompress the spectral density corresponding larger matrin size 2 nay
similor Figure 1(c) rhoglarge, x ff.decompress plot-True)
35

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: Proofs of main theoretical claims are provided in the appendix. Numerical evidence is provided in the body of the paper (with additional expermiments and data in the appendices).
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification:
limitations section is provided.
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (c.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The au- thors should reflect on how these assumptions might be violated practice and what the implications would be.
The authors should reflect the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
applicable, the authors should discuss possible limitations of their approach to ad- dress problems of privacy and faimeess.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community, Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
36

---

Answer: [Yes]
Justification:
An extensive appendix provides all proofs and relevant assumptions
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proofs can either appear in the main paper the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should complemented by formal proofs provided in appendix supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
Justification:
Code is also provided anonymized repository.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
the contribution is dataset and/or model, the authors should deseribe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or the contribution a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, other means that are appropriate to the research performed.
While NeurIPS does not require releasing code. the conference does require all sub- missions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) If the contribution is primarily a new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is new model (e.g., a large language model), then there should either be a way to access this model for reproducing the results or a way to re- produce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors welcome describe the particular way they provide for reproducibilityy In the case of closed-source models, it may be that access to the model is limited in some way (e.g., registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
37

---

reproducible, they Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
Code is provided in an anonymized repository. Where possible, datasets are made available, with the exception the full NTK model required for model verification (due to the nature of the work and space constraints).
Justification:
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips. cc/
solosoo for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https: / /nips.
/Coaaee/dooSSSbbiccc onPolicy) for more details.
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data,
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permirted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: Hyperparameters reported in the main document. All other details either previously reported or contained in the provided codebase.
Guidelines:
The answer NA means that the paper does include experiments.
The experimental setting should be presented in the core of the paper to a level of detail that is necessary appreciate the results and make sense of them.
The full details can provided either with the code, appendix, or as supplemental material.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: A study of sensitivity to different subsamples is provided in the supplementary material.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
38

---

The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report 1-sigma error bars, but one should state it. The authors should prefer- ably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not show in tables figures symmetric error bars that would yield results that are out of range (e.g. negative error rates).
error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification: Compute resources used were reported in the main document.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount compute required for each of the individual experimental runs as well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make into the paper).
9. Code ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics https: atttps://neurii /tt.rps..p.pp//
Answer: [Yes]
Justification: The authors have reviewed the NeurIPS Code of Ethics and believe that conform to this code.
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (c.g.. if there special consid- eration due to laws or regulations in their jurisdiction..
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification:
The authors believe there is societal impact of the work performed.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
39

---

If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact spe- cific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when technology is being used intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
The paper poses no such risks
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines restrictions to access the model or implement- ing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification:
The software provided the anonymized repository was authored by the paper's s authors. Datasets used were open access or computed by the authors from open access datasets.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (c.g., CC-BY 4.0) should be included for each asset.
40

---

For scraped data from a particular source (e.g., website), the copyright and terms of service of that souree should be provided.
assets are released, the license, copyright information, and terms use in the pack- age should be provided. For popular datasets, paperswi whccode com dataseta has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's s creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: Software and documentation are provided at an anonymized repository Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moodl as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification:
This is theoretical and computational study
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contri- bution of the paper involves human subjects, then much detail as possible should be included in the main paper.
According to the NeurlPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals equivalent for research with human subjects
Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Question:
Answer: [NA]
Justification:
This is theoretical and computational study
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
41

---

Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Does the paper describe the usage of LLMs if is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, or formatting purposes and does not impact the core methodology, scientific rigorousness, originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification: LLMs were not used as a core part of this work.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, or non-standard components.
Please refer to our LLM policy (https:/ /ttpar/////ui...eeeeeee 20/C0nneeeeencce/22255////////// for what should or should not be described.
42


Reviews:
[{"review_id": "EP5HFhdJxm", "reviewer": "Reviewer_gtu3", "summary": "In this paper, the authors consider a problem of spectral estimation, where the goal is to estimate the spectrum of a large matrix from that of a small submatrix. The authors propose an algorithm based on the (time) evolution of the Stieltjes transform. Several mathematical properties of the proposed algorithm are proved, and the performance of the algorithm is tested numerically by synthetic experiments and also several datasets.", "questions": "(1) Theoretically, to obtain the spectral density from the Stieltjes transform, one can use the Stieltjes inversion that computes $(1/\\pi) \\lim_{\\epsilon \\to 0} \\textrm{Im } m(x+i\\epsilon)$. I wonder if one can avoid extending the Stieltjes transform to the lower half plane $C^-$ by applying the Stieltjes inversion.\n\n(2) Some random matrices have outlier eigenvalues (for example, due to the presence of underlying low-rank structure), which are crucial in understanding the properties of the random matrices. Are these outlier eigenvalues ignored when approximating the empirical spectral density?\n\n(3) Minor comment: The equation number is missing in line 211.", "limitations": "The limitations of the work are well addressed. It seems that there would be no negative societal impact of the work.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "Strengths: \nThe problem considered in this paper is of great importance, and the proposed algorithm is novel. The idea is simple, but it has high potential to be developed further since the Stieltjes transform is well-studied object in random matrix theory. The writing is clear in general.\n\nWeaknesses: \nSeveral assumptions, especially one on the freeness, are rather too strong.", "was_revised": true, "final_justification": "I thank the authors for the detailed replies. I would keep my scores.", "rebuttal": "We thank the reviewer for their thoughtful review and constructive feedback. Below we address the specific questions and concerns raised.\n\n---\n### Clarification on Freeness Assumptions\n\nWe believe there may be a misunderstanding of the freeness assumption. We do *not* assume that the data matrix is free with a fixed or arbitrary companion matrix. Rather, we rely on a well-established result from free probability: the *asymptotic freeness* of the pair $(\\mathbf{E}, \\mathbf{P_\\sigma A P_\\sigma^{-1}})$, where $\\mathbf{E}$ is the projection onto the top-left  $n_s \\times n_s$ block, $\\mathbf{P}_\\sigma$ is a Haar-uniform random permutation matrix, and $\\mathbf{A}$ is *any* deterministic Hermitian matrix with uniformly bounded operator norm. \n\nNica (1993, p.5) shows that the random permutation induces asymptotic freeness in the large-$n$ limit. This setting is in fact one of the cleanest regimes in which freeness can be rigorously justified and is widely used in free-probability-based methods, requiring no additional structural assumptions on the data matrix beyond a standard bounded-norm condition.\n\n---\n### Analytic continuation and the Stieltjes inversion formula\n\nThe reviewer is correct that the Stieltjes inversion formula can be used to obtain the spectral density $\\rho(t, x)$ from the Stieltjes transform $m(t, z)$, and this is exactly how we extract the density in the final step of our method (see line 10 of Algorithm 2). \n\nHowever, the analytic continuation of $m_0(z)$ into the lower half-plane is not related to this inversion step, but is instead required *earlier* in the pipeline: when solving the PDE that evolves $m(t, z)$ from the initial data $m_0(z)$.\n\nIn our method, the evolution is performed via the method of characteristics, where the evaluation point $z_0$ flows across the support of the spectral density and enters the lower half-plane. Since $m_0(z)$ is discontinuous across the support, its values in the lower half-plane cannot be obtained from upper half-plane limits alone. For this reason, the analytic continuation plays a crucial role in solving the PDE, rather than in the inversion itself.\n\nMore broadly, we spent significant effort exploring several alternative approaches to avoid analytic continuation. In particular, we tried solving the PDE directly with standard finite‐difference and finite‐element schemes: these worked for small values of $t$ but became numerically unstable as $t$ grew. This is because the solutions close to the real line that are used in this inversion formula necessarily rely on information from the lower half-plane. This difficulty is made clear by the form of the solutions in the method of characteristics, and is quantitatively encoded in Proposition 4. Consequently, and despite the need for analytic continuation, the method of characteristics has been the most effective method we have tried in practice. We are also investigating distribution-specific solvers that may bypass this difficulty in certain cases, although such methods are currently limited to special families and remain an area of ongoing work.\n\nWe have added these clarifications to the text.\n\n---\n### Outlier eigenvalues and low-rank structures\nWe also agree with the reviewer that the presence of (almost) low-rank structure and other outlier eigenvalues (such as heavy-tailed behavior) are crucial to understanding the properties of many matrices of interest. Indeed, doing so is a strong motivator for our own study and development of this method. \n\nOur current work focuses on modeling the *bulk of the spectrum*, using the absolutely continuous part of the empirical spectral density. Outliers in the form of isolated eigenvalues without mass (e.g., the top eigenvalue of a finite-rank perturbation) are indeed ignored in our current setting, as they correspond to *measure-zero* features that vanish under the Cauchy integral defining the Stieltjes transform, and therefore do not enter the PDE.\n\nHowever, *spikes with non-zero mass* (that is, Dirac components in the spectral density) can in principle be handled. The free decompression PDE preserves such components and governs the evolution of both their location and mass. We did not incorporate outliers or spikes into the current work, as their proper treatment involves additional modeling considerations that warrant a dedicated future investigation.\n\nWhile we have preliminary results in these directions, we believe the present paper already addresses a broad and important class of matrices and is a natural first step in tackling this problem. We will add a short discussion of these extensions in the revised version.\n\n---\n\nWe would again like to thank the reviewer for their constructive feedback. We would be happy to provide any additional clarification the reviewer might find helpful.\n\n---\n\n#### Reference\n\n[1] Nica, A. (1993). Asymptotically free families of random unitaries in symmetric groups. Pacific Journal of Mathematics, 157(2)."}, {"review_id": "oo635EeBXD", "reviewer": "Reviewer_BF8q", "summary": "This paper introduces \"free decompression,\" a novel method for estimating the spectral density of large impalpable matrices. The method exploits tools from free probability theory, specifically evolving the Stieltjes transform of a small submatrix via a partial differential equation to estimate the spectrum of the full matrix. The paper provides strong theoretical foundations, polynomial error bounds, and justifies the approach on synthetic examples and real-world datasets.", "questions": "1) Proposition E.1 in the appendix seems to be one of the critical results and I have several questions here. line 754, what do you mean by first order approximation (is this similar to Taylor/Laurent type expansion and in that case what is the region of validity: radius of convergence, ie. an upper bound on $|\\delta m|$)? I can clearly see $m(z)$ is meromorphic, therefore $m$ does have a finite Laurent expansion near the poles (in fact they all seem to be simple poles), however, it is not clear to me what exactly is $\\delta m$ and how are you bounding the higher order errors here. \n2) line 772: could you prove in detail why $C_y <\\infty$ please for all $y$ (in particular why those supremums are finite)? I want to understand how large $C_y$ can get as this is essentially the exponent $\\nu.$\n3)  How sensitive is the method to the choice of submatrix? Are there better strategies than random sampling? \n4) These two questions are not deal-breakers and I understand they could be well out of reach for now: Is there a (non-trivial) lower bound on $\\delta m_n$ of polynomial (or poly-log) in $\\dfrac{n}{n_s}$? Finally, what about non-hermitian matrices that do occur frequently in practice?", "limitations": "yes", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 2, "significance": 4, "originality": 4, "strengths_and_weaknesses": "Strengths:\nAs already mentioned, the paper is thoeretically very strong, including detailed formal propositions and proofs for the evolution equation and error bounds. Unlike purely theoretical work, the authors provide a complete implementation with detailed algorithms and a Python package, making their method more accessible. The experimental validation includes multiple synthetic distributions and challenging real-world applications (for instance Facebook network). Finally, the paper is very well-written with good motivation and helpful visualizations. This paper overall makes a sound contribution to an important problem in numerical linear algebra.\n\nWeaknesses:\nThe authors emphasizes that despite the limitation of their proposed work to real-world application, there is no competitive alternative. However, such a strong statement probably needs a bit more justification. Could you provide some non-trivial warning signs (i.e. classes of matrices where free decompression is expected to fail except obvious highly-ill conditioned ones) for practitioners and how do the alternative methods compare in such cases? Apart from these, there are a few more things I would like to understand clearly before recommending a full acceptance (please see the questions).", "was_revised": true, "final_justification": "This is a very good paper and I am satisfied with authors response. I am willing to increase my rating if AC deems it necessary.", "rebuttal": "Thank you for the careful review of our paper. We are glad to receive your positive appraisal, and we have provided responses to your questions and concerns below.\n\n---\n\n### On the Scope and Comparison with Alternatives\n\nRegarding competitive alternatives, the issue is not that we view existing methods as uncompetitive, but rather that we could not find any direct alternatives at all. That is, algorithms that accept a matrix as input and return an approximation of the spectral density of a much larger matrix of the same type, containing the original matrix in its upper-left corner. To our knowledge, the only existing approach is to appeal to limiting cases-e.g., approximating the spectral distribution by a Marchenko-Pastur law. While this can be appealing in some settings, we have found that real datasets are often ill-conditioned and heavy-tailed, requiring more flexible approximations.\n\nSince our method approximates the spectral density, functionals can be estimated directly via quadrature. This makes the FLODANCE algorithm (Ameli et al. 2025) a relevant indirect competitor. A detailed comparison is provided in our response to reviewer *EuBy*, to which we kindly refer the reviewer.\n\nIn our internal experiments, we have also identified several challenging matrix classes. The reviewer is correct that the most prominent of these is pathological ill-conditioning, particularly in low-rank matrices. Another difficult case arises when the spectrum has outliers or large spikes caused by structural features of the matrix. Finally, distributions with disjoint support (i.e., nonzero on multiple separate intervals) can also pose difficulties. Each of these cases requires specialized tools that are currently under active development and will be the focus of future work.\n\n---\n\n### Clarification on First-Order Approximation\n\nTo analyze the stability of the nonlinear PDE, we consider a perturbation of the solution of the form\n$$\\tilde{m}(t, z) = m(t, z) + \\delta m(t, z),$$\nwhere $\\delta m$ is the first-order variation of $m$. This is not a Taylor or Laurent expansion in $z$, but rather a first-order (in $\\delta m$) *functional variation*, corresponding to a linearization of the PDE. Substituting $\\tilde{m}$ into the nonlinear PDE and retaining only the terms linear in $\\delta m$ yields a linear PDE for $\\delta m$—the equation obtained by taking the Gateaux derivative of the original PDE operator, a standard technique in PDE stability analysis. This linearized PDE governs how perturbations to the initial condition $m_0$ propagate over time.\n        \nImportantly, this analysis is carried out for $z$ in the upper half-plane $\\mathbb{C}^+$, where $m(t, z)$ is a Herglotz function (i.e., holomorphic with positive imaginary part) and has no poles—hence, it is not meromorphic. (Poles arise only on the second Riemann sheet used during the characteristic continuation into the lower half-plane.) As such, notions like Laurent expansion, radius of convergence, and pole structure do not apply, since the Stieltjes transform is analytic on $\\mathbb{C}^+$. The perturbation $\\delta m$ is likewise defined on $\\mathbb{C}^+$ and assumed to remain small in the $\\|\\cdot\\|_\\infty$ norm. Proposition E.1 then quantifies how such perturbations evolve, bounding their growth via a Gronwall-type estimate.\n\n---\n\n### Clarification on Boundedness of $C_y(t)$ in Proposition E.1\n\nWe thank the reviewer for pointing this out, as we agree that further clarification is needed. In Remark E.1, we briefly sketched why $C_y(t)$ in Proposition E.1 remains finite for all $0 < y < \\infty$, but we agree that the argument was too coarse and should be made explicit. The key bounds follow from standard properties of Stieltjes transforms of compactly supported densities.\n\nRecall from equation (E.3) that $C_y(t)$ is composed of two terms:\n$$\\sup_{z \\in \\Gamma_y} \\left| \\Re \\frac{\\partial_z m}{m^2} \\right|, \\quad \\text{and} \\quad \\sup_{z \\in \\Gamma_y} \\left| \\Im \\frac{1}{m} \\right|.$$\nWe now verify that each of these terms is finite for all $0 < y < \\infty$. Let $m(t,z) = \\int_{\\mathbb{R}} \\frac{\\rho(t,\\xi)}{\\xi - z} \\mathrm{d}\\xi$ with $z = x + i y \\in \\mathbb{C}^{+}$. Then:\n\n- Since $\\rho(t,x)$ is supported in a compact interval $[\\lambda_{-}, \\lambda_{+}]$ of diameter $D = \\lambda_{+} - \\lambda_{-}$, we have\n$$|m(t,z)| \\geq \\Im m(t,z) = \\int_{\\mathbb{R}} \\frac{y \\rho(t,\\xi)}{(\\xi - x)^2 + y^2}  \\mathrm{d}\\xi \\geq \\frac{y}{D^2 + y^2} \\int \\rho(t,\\xi) \\mathrm{d}\\xi = \\frac{y}{D^2 + y^2}.$$\n- It follows that $|m(t,z)^{-1}| \\leq \\frac{D^2 + y^2}{y}$, so $|\\Im m(t,z)^{-1}|$ also has this bound.\n- For the derivative, note\n$$\\partial_z m(t,z) = \\int_{\\mathbb{R}} \\frac{\\rho(t,\\xi)}{(\\xi - z)^2} \\mathrm{d}\\xi,$$\nso\n$$|\\partial_z m(t,z)| \\leq \\int_{\\mathbb{R}} \\frac{\\rho(t,\\xi)}{(\\xi - x)^2 + y^2}  \\mathrm{d}\\xi \\leq \\frac{1}{y^2}.$$\n- Combining the above, we obtain\n$$\\left| \\frac{\\partial_z m(t,z)}{m(t,z)^2} \\right| \\leq \\frac{1}{y^2} \\left( \\frac{D^2 + y^2}{y} \\right)^2 = \\frac{(D^2 + y^2)^2}{y^4},$$\nwhich is finite for every $y > 0$.\n    \nThis confirms that $C_y(t) < \\infty$ for all $y > 0$. We will be happy to include this derivation in the revised appendix for completeness.\n\nNote that this bound is not defined at $y = 0$, which is expected: as $y \\to 0^+$, the Stieltjes transform approaches the real axis where it develops discontinuities, and quantities like $\\partial_z m$ may diverge. However, in both our theory and implementation, we evaluate $m(t, z)$ slightly above the real axis (e.g., at $z = x + i \\delta$ for fixed $0 < \\delta \\ll 1$), where all terms remain finite. In practice, the observed error propagation is significantly milder than this worst-case upper bound suggests.\n\n---\n\n### Sampling Considerations\n        \nThe decompression equation (Theorem 1) holds on the assumption that the submatrix is sampled from a Haar-uniform permutation of the full matrix. That is, the rows and columns are first permuted uniformly at random, and then the top-left $n_s \\times n_s$ block is selected. This ensures that the submatrix and its complement are asymptotically free in the large-$n$ limit.\n\nAny form of deterministic, importance-weighted, or leverage-score sampling would break this freeness assumption. Extending the theory to account for such alternatives would require a new analysis to show that the submatrix remains (approximately) free with respect to full matrix, which is an open problem to our knowledge. Therefore, uniform random sampling is not an arbitrary design choice; it *makes the theory valid*.\n        \nFrom a practical standpoint, we find the method is robust to the submatrix sampling. To quantify this, we repeated our experiment on the NTK matrix using *twenty* independent Haar-uniform submatrices, each of size $n_s = 2^{10}$. For each batch, we applied free decompression to reconstruct the spectral density at increasing matrix sizes $n = 2^{11}, 2^{12}, \\dots, 2^{15}, 50000$, and compared the result against the exact density at that scale.\n        \nWe measured the deviation between the estimated and true densities using multiple metrics, including total variation (TV), Jensen–Shannon (JS) divergence, and Kolmogorov–Smirnov (KS) statistic, which take values in the range $[0, 1]$. For each metric, the *mean and standard deviation across batches* are reported below. These results demonstrate that the method is *practically insensitive* to the specific submatrix; even a single random batch provides a sufficiently accurate spectral estimate.\n        \n| Size $n$ $\\quad$ | TV ($\\times 100$) $\\quad$| JS ($\\times 100$) $\\quad$| KS ($\\times 100$) $\\quad$| $\\Delta \\mu_1 / \\mu_1$ ($\\times 100$) | $\\Delta \\mu_2 / \\mu_2$ ($\\times 100$) |\n|-----------|-------------------|-------------------|-------------------|-------------------|-------------------|\n| $ 2^{10}$ |   $1.08 \\pm 0.08$ |   $3.90 \\pm 0.55$ |   $0.42 \\pm 0.05$ | $0.33 \\pm 0.20$ | $0.31 \\pm 0.21$ |\n| $ 2^{11}$ |   $1.07 \\pm 0.16$ |   $3.89 \\pm 0.76$ |   $0.42 \\pm 0.07$ | $0.13 \\pm 0.12$ | $0.13 \\pm 0.13$ |\n| $ 2^{12}$ |   $2.01 \\pm 0.30$ |   $4.69 \\pm 0.60$ |   $0.70 \\pm 0.12$ | $0.06 \\pm 0.04$ | $0.07 \\pm 0.04$ |\n| $ 2^{13}$ |   $3.18 \\pm 0.51$ |   $6.23 \\pm 0.44$ |   $1.15 \\pm 0.19$ | $0.02 \\pm 0.02$ | $0.02 \\pm 0.02$ |\n| $ 2^{14}$ |   $4.26 \\pm 0.84$ |   $7.46 \\pm 0.80$ |   $1.72 \\pm 0.45$ | $0.01 \\pm 0.01$ | $0.01 \\pm 0.01$ |\n| $ 2^{15}$ |   $5.12 \\pm 1.20$ |   $7.88 \\pm 1.15$ |   $2.48 \\pm 0.82$ | $0.03 \\pm 0.01$ | $0.04 \\pm 0.02$ |\n| $50$K     |   $5.91 \\pm 1.52$ |   $8.26 \\pm 1.47$ |   $3.00 \\pm 1.13$ |  $0.19 \\pm 0.03$ | $0.07 \\pm 0.07$ |\n\n---\n\n### Regarding a possible lower bound\n        \nA universal *lower* bound for $\\Vert \\delta m(t)\\Vert_{L^2 (\\Gamma_y)} $ cannot be given: the perturbation $\\delta m(0,\\cdot)$ can be chosen arbitrarily small (or identically zero), so the sharp bound is simply $\\Vert\\delta m(t)\\Vert_{L^2 (\\Gamma_y )} \\ge 0$.\n        \nAny positive lower bound would depend on additional information about the initial perturbation, not on the PDE itself. For stability and error-propagation purposes the relevant result is therefore the upper (growth) estimate provided in Proposition E.1.\n\n### Non-Hermitian Matrices\n        \nOur current framework is designed specifically for Hermitian matrices. Extending this framework to non-Hermitian matrices (whose spectra lie in the complex plane) is an interesting and important direction. However, doing so would require substantial modifications, including different analytic machinery and stronger structural assumptions. We consider this a valuable but nontrivial extension, and leave it to future work.\n\n---\n\nWe appreciate the reviewer’s constructive feedback and would be happy to provide further clarification or make additional improvements based on any remaining questions or suggestions."}, {"review_id": "PTDvhb1LBF", "reviewer": "Reviewer_EuBy", "summary": "This paper proposes an algorithm to estimate the spectral density of a very large matrix.\n\nThey consider the extreme case of matrices which are too large to be formed in their entirety, either implicitly or explicitly. They call such matrices *impalpable*, and differentiate them from other \"less difficult\" classes of matrices in Section 2. They work under the assumption that there is an underlying $N \\times N$ matrix of which we want to estimate the spectral density, and that one can access an $n \\times n$ submatrix of it which is obtained from the intersection of the row and columns corresponding to a randomly sampled an index set. They observe that this assumptions correspond to the assumption that these sub-matrices form a sequence (indexed by $n$, and assuming $N$ is infinite if I have understood correctly) which is asymptotically free.\n\nThe key idea of this paper is that this assumption provides access to tools from random matrix theory and free probability theory to extrapolate the spectral density of the $n \\times n$ sub-matrix to the full $N \\times N$ matrix.\n\nSection 3 introduces these tools. Namely, the Stieltjes transform, its inverse, the R-transfrom, and a theorem of Nica and Speicher which relates the R-transform of a free random matrix to that of a fixed submatrix. This theorem, which tells you how to glean information about the spectral density of a random submatrix of a large full matrix, is known as free compression. The authors propose performing this operation is reverse, to glean information about the spectral density of the full matrix from a smaller random sub-matrix. They call this *free decomposition*.\n\nDirectly implementing this operation involves inverting the Stieltjes transform in order to compute the R-transform. This is not generally analytically tractable. The bulk of the technical contribution of this paper is to propose a numerical method to perform free decompression directly using the Stieltjes transform.\n\nSection 3.2 shows that free decomposition can be written as the evolution of a PDE in the Stieltjes transform (Prop. 1), that this can be easily solved (Prop. 2), but that the error may grow \"polynomially fast\" (Prop. 3), suggesting the need to high precision computation. Section 3.3 provides a very nice informal, intuitive explanation of the previous section.\n\nSection 4 explains some practical challenges of using Propositions 1 and 2 to estimate the Stieltjes transform under free decompression and proposes some practical numerical solutions. Algorithm 1 on page 7 provides psuedo-code for the full proposes algorithm. \n\nSection 5 provides some numerical experiments on synthetic and real data, comparing the time and accuracy of the proposed method.\n\nThe authors provide a Python package to implement their algorithm.", "questions": "- For the task of computing log-determinants of impalpable matrices, how does your method compare to the FLODANCE algorithm of Ameli et al. (2025)?\n- What impact do you imagine this algorithm will have on machine learning? Can you give an example of a machine learning algorithm where this algorithm could replace an existing one to improve its efficiency?", "limitations": "Yes. The limitations are discussed briefly yet comprehensively.", "rating": 6, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 4, "originality": 4, "strengths_and_weaknesses": "This paper is highly innovative and the authors should be commended for their level of creativity in pulling all these ideas together to provide a solution to a problem which is both highly practically important and which currently has very few tenable solutions. Despite the fact the computing the spectral density of very large \"impalpable\" matrices is an extremely important problem, before this paper, there were effectively no algorithms to do this.\n\nThere are three main innovations which lead to the algorithm: 1) the framing of the matrix as a free random matrix and subsequently framing estimating its spectral density as \"free decompression\"; 2) describing \"free decompression\" as a PDE; and 3) coming up with the necessary continuations and smoothings required to numerically solve the PDE. To me, this is an extremely elegant approach and I find it difficult to find fault in it.\n\nThe empirical experiments are extensive and compelling.\n\nI believe the algorithm proposed in this paper has the potential to be highly impactful in the field of numerical linear algebra, and may open the possibility of new breakthrough which were previously impossible without such an algorithm.\n\nThe paper is extremely well-written. Though technical, the core ideas are presented in an accessible manner and the reader is guided through all of the ideas and innovations in a very clear and helpful way.", "was_revised": true, "final_justification": "I continue to strongly support this paper.", "rebuttal": "We thank the reviewer for their enthusiastic and thoughtful feedback. We are delighted that you found the contributions and presentation of our work compelling. Below, we address your comments and provide additional details that we hope you will find helpful.\n\n---\n\n### Computing Log-Determinants\n\nTo address the reviewer’s question, we conducted a direct numerical comparison between free decompression and the FLODANCE algorithm (Ameli et al., 2025) on the same input matrix used in our experiments (log-NTK from CIFAR-10 with ResNet50). Starting from an initial sub-matrix of size $n_s = 2^{10}$, we estimated the log-determinant $\\ell_n = \\log\\det(\\mathbf{A}_n)$ at increasing matrix sizes, up to $n = 50{,}000$, and report the relative errors against ground truth in the table below.\n\n\n| Size $n$           |  $2^{10}$ |  $2^{11}$ |  $2^{12}$ |  $2^{13}$ |  $2^{14}$ |  $2^{15}$ |     $50$K |\n|--------------------|-----------|-----------|-----------|-----------|-----------|-----------|-----------|\n| Free Decompression | $0.001$%  | $0.007$%  | $0.018$%  | $0.016$%  | $0.022$%  | $0.027$%  | $0.193$\\% |\n| FLODANCE           | $0.006$%  | $0.283$%  | $0.984$%  | $2.180$%  | $5.401$%  | $5.128$%  | $4.412$\\% |\n\nAt this range of matrix sizes, free decompression yields consistently more accurate estimates. We attribute this to a key difference in assumptions: FLODANCE relies on the stationarity of a log-normal stochastic process governing the determinants of growing sub-matrices. This assumption typically becomes valid at much larger scales; in the FLODANCE paper, the authors mitigate early-stage deviations by introducing a burn-in period and applying the method only after the stochastic process stabilizes. In contrast, free decompression does not rely on such stationarity and remains accurate even at smaller sizes, provided the spectrum of the initial matrix is large enough to give an accurate approximation of the spectral distribution.\n\nWhile FLODANCE is well suited to ultra-large matrices where stationarity becomes apparent, free decompression offers a complementary approach that does not rely on such assumptions and remains accurate without requiring burn-in or stochastic modeling. We also note that, using the free decompression PDE, functionals of the distribution can be approximated by ODEs. Consequently, one can derive a continuous-time ODE for the log-determinant $\\ell_n$. These types of problems provide rich avenues for future work. \n\n---\n\n### Impact and Applications\n    \nThe ability to reconstruct the full eigenvalue spectrum of a large matrix enables the approximation of a broad class of linear algebraic quantities that appear throughout modern machine learning. We categorize these quantities below according to the *type of spectral functionals* they involve, rather than by application domain, since the same mathematical structures arise in diverse problems.\n    \n#### Matrix Functionals\n    \nFree decompression (FD) approximates spectral quantities of the form\n$$f(\\mathbf{A}_n) = n \\int f(\\lambda) \\rho(\\lambda)  \\mathrm{d}\\lambda,$$\nwhere $\\rho(\\lambda)$ is the empirical spectral density of $\\mathbf{A}_n$. This framework encompasses a wide class of spectral functionals, including log-determinants $\\log \\det(\\mathbf{A}_n + \\sigma^2 \\mathbf{I})$, traces of inverse powers such as $\\mathrm{tr}(\\mathbf{A}_n + \\sigma^2 \\mathbf{I})^{-1}$ and $\\mathrm{tr}(\\mathbf{A}_n + \\sigma^2 \\mathbf{I})^{-2}$, and more generally, any analytic function $f$. These expressions arise in numerous applications, such as:\n\n- Gaussian-process marginal likelihood and predictive variance,\n- Determinantal Point Process (DPP) training and sampling,\n- Laplace-approximated Bayesian inference (e.g., for GLMs and neural networks),\n- PAC-Bayes generalization bounds in deep learning,\n- Kalman filtering and EM-based state-space models.\n    \nThis category also includes spectral entropy measures such as the Von Neumann entropy,\n$$S_{\\mathrm{VN}} = - n \\int \\lambda \\log (\\lambda)  \\rho(\\lambda) \\mathrm{d}\\lambda,$$\nwhich appears, for instance, in graph neural network pipelines as a measure of graph complexity or regularization. Extensions to Rényi or Tsallis entropies are similarly supported. Other notable examples include heat kernel traces $\\mathrm{tr} e^{-t\\mathbf{A}_n}$, log-determinants of Laplacian minors (e.g., spanning-tree counts), and spectral divergences (e.g., KL or Wasserstein) used in domain adaptation. All of these reduce to integrals of the form $\\int f(\\lambda) \\rho(\\lambda)  \\mathrm{d}\\lambda$, and hence fall naturally within the scope of FD. Once $\\rho(\\lambda)$ is recovered, such functionals can be computed efficiently via one-dimensional quadrature.\n    \n\n#### Quadratic Forms\n    \nFD can also be used to approximate quadratic forms involving (powers of) inverse matrices:\n$$\\boldsymbol{y}^{\\intercal} (\\mathbf{A}_n + \\sigma^2 \\mathbf{I})^{-1} \\boldsymbol{y}.$$\nWhile such expressions typically require a full solve, FD can efficiently estimate *expected* values or averages over random probes. This is particularly useful when one requires many such forms (e.g., in kernel ridge regression, influence functions, or dropout variance estimation). For a specific fixed vector $\\boldsymbol{y}$, FD is not directly applicable in its standard form, though natural extensions will allow such quantities to be estimated within a modified framework, and are the focus of ongoing work.\n    \n\n#### Condition Number Estimation\n    \nThe condition number $\\kappa(\\mathbf{A_n}) = \\lambda_{\\max}/\\lambda_{\\min}$ can be approximated by extrapolating the spectral edges of $\\rho(\\lambda)$. FD is effective at recovering the *bulk support* of the spectrum, which can provide reliable bounds in many practical scenarios. However, isolated outlier eigenvalues may not be accurately captured, so combining FD with one or two power iterations is advisable when precise extremal estimates are required. This is particularly relevant in diagnostics of ill-conditioning and edge-of-stability analyses.\n\n---\n    \nWe would be happy to include these categories and examples as a short section in the appendix if the reviewer finds this exposition helpful. We welcome any further questions or suggestions the reviewer may have and would be happy to revise the paper accordingly."}, {"review_id": "YUBADjCiAz", "reviewer": "Reviewer_GMts", "summary": "This paper studies algorithms that attempt to estimate certain statistics of large matrices that could not even be written down or accessed via matrix-vector product. In this impalpable model, only a small sub-matrix could be accessed by the algorithm, and the main task studied in this paper is to estimate the spectral density. The major idea is a technique called free decomposition, which is essentially the R-transform from random matrix theory. In short, the R-transform of the top-left block could be connected to the R-transform of the whole matrix. Note that this standard result does not restrict the size of the whole matrix, therefore one could view this as an algorithm for the impalpable matrix setting. Of course, computing the R-transform requires Stieltjes transform which the authors shown that could be approximated via a first-order PDE, and they further develop an algorithm to solve this PDE. Finally, experiments are performed to validate the efficiency of the proposed algorithm for large matrix.", "questions": "Can you provide an overview of algorithm 1, and subroutines used in it?\n\nThere are compilation errors in the pdf, e.g.,\n\n* Line 211, \"given by (??)\".", "limitations": "Yes.", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 2, "originality": 2, "strengths_and_weaknesses": "Strengths:\n\n1. Studying algorithm in the impalpable matrix setting is interesting, as modern day machine learning operates on larger and larger matrices. As authors noted, random matrix theory and free probability would be a crucial tool as it studies the asymptotic behavior of random matrices, which could be \"inverted\" to develop a tool in this setting.\n\n2. The main contribution, in my opinion, is the algorithm developed in section 4 to practically approximate the Stieltjes transform. Several ideas, including finding a different branch of the transform and using Jacobi polynomials to estimate the empirical spectral density, are quite interesting. This could lead to more practical algorithms for free probability and RMT.\n\nWeaknesses:\n\n1. The presentation, especially the first 9 pages, could be significantly improved. Many notations are not defined, e.g., $L^\\infty, L^2(\\Omega)$ and subroutines in algorithm 1 are mostly left unexplained. This is particularly confusing for algorithm 1, for example, what is EigApprox? What is Glue? How are these two used for computing the Stieltjes transform? Without reading the appendix and digging into the details, it's very hard to grasp. Also, it feels to me that 3.3 is a bit redundant, as it serves as an intuition for the proofs in appendix C, which are not present in the main text anyway. Finally, discussions in section 4 are highly technical, as a reader with some knowledge about RMT, I found it quite hard to understand and follow, in particular section 4.1. In its current form, I feel the paper targets the audience with deep expertise in RMT, not the general audience of NeurIPS.\n\n2. On the other hand, this paper could be treated as an elementary application of RMT and free probability (in particular the R-transform), therefore I feel its theory is not very novel. However, as I've stated above, I'm only a reader with some knowledge about RMT, so reviewers with more expertise in this area should weigh in and judge the particular contribution of the paper. I'm willing to change my ratings based on the opinions from reviewers with more expertise in RMT.", "was_revised": true, "final_justification": "After reading authors' detailed responses and other reviewers' comments, I believe authors' contribution, even from an algorithmic RMT perspective is quite significant. Also, authors proposed reasonable clarifications for their algorithm, and I hope they could incorporate them in their revision.", "rebuttal": "We thank the reviewer for their thoughtful review and constructive feedback. We hope the following clarifications fully address the concerns raised and improve the overall clarity and contribution of the manuscript.\n\n---\n\n### Clarity and Notation Improvements\n\nWe thank the reviewer for their comments regarding notation and presentation. To address these concerns, we have added a new Appendix A that includes:\n \n- A comprehensive nomenclature table defining all commonly used variables and symbols, including previously undefined terms such as $L^2(\\Omega)$ (the space of square-integrable functions) and $L^\\infty$ (the space of essentially bounded functions).\n- A brief overview of key functional transforms from free probability (Stieltjes transform, Hilbert transform, and R-transform), aimed at readers without prior expertise in random matrix theory.\n        \nIn the main text, Section 4 now begins with a short paragraph outlining the structure of Algorithm 1. This includes clarifying the purpose of each component and distinguishing which subroutines are essential to understanding the core algorithm, versus those discussed in greater detail in the appendix.\n        \nSpecifically:\n\n- **EigApprox** refers to the numerical approximation of the smoothed spectral density of the input matrix. Although seemingly routine, this step requires high precision, as the accuracy of the entire pipeline depends on it. We use Chebyshev or Jacobi polynomial expansions to ensure spectral accuracy.\n- **Glue** refers to a procedure inspired by Riemann–Hilbert techniques. It is used to reconstruct the second-sheet branch of the Stieltjes transform from its boundary values, effectively gluing together the analytic continuations across the support interval.\n        \nRegarding Section 3.3, we have decided to retain it. Although its role is primarily intuitive, other reviewers found it helpful for navigating the technical content that follows. We believe it supports accessibility rather than introducing redundancy.\n        \nWe hope these revisions improve the clarity and accessibility of the manuscript for a broader readership, including those less familiar with RMT.\n\n---\n\n### Clarifying Theoretical and Practical Contributions\n        \nWe appreciate the reviewer’s candid reflection and their openness to reconsider their rating. While the foundational properties of the R-transform are indeed classical, our work builds upon these to develop a novel pipeline for estimating the spectral density of large matrices from small submatrices—an inverse problem that, to our knowledge, has not been previously realized in this form.\n\nOur theoretical contributions go beyond a direct application of known results. In particular:\n\n- We develop a nonlinear PDE that characterizes the evolution of the Stieltjes transform under matrix growth. While this PDE can be formally derived from classical properties of the R-transform, we believe its explicit formulation, analysis, and practical use as a reverse compression method—what we term free decompression—are new to the literature.\n- Proposition 4 reveals that every characteristic trajectory necessarily crosses the real axis, which creates a bottleneck in the numerical solution. This insight motivated our use of analytic continuation, an essential step for solving the PDE that we believe is novel in this context.\n- To achieve this continuation, we introduce a numerical procedure inspired by scalar Riemann–Hilbert problems, and leverage gluing two branches and series acceleration techniques (e.g., Wynn’s $\\epsilon$-algorithm). These techniques, while classical on their own, are, to our knowledge, applied here for the first time.\n- Appendix E provides a new quantitative analysis (Proposition E.1) showing how errors in the initial spectral approximation propagate along characteristics, complementing the closed-form characteristic solution and polynomial stability bound in Proposition 3.       \n- The precision approximation of the Stieltjes transform via Chebyshev or Jacobi polynomial expansions is engineered for accuracy and scalability. The reviewer also acknowledged this as one of the strengths of the paper.\n        \nBeyond theory, we also emphasize the practical realization of these ideas. We demonstrate accurate decompression on challenging real-world datasets, including the Facebook graph Laplacian and neural tangent kernels (NTKs), and release an open-source implementation to facilitate further use by the community. To our knowledge, this is the first instance of a free-probability-based decompression method—at any scale—being successfully applied to real-world matrices of this size and complexity.\n        \nWe hope this clarifies that the paper is not merely an application of known theory, but offers new insights, tools, and practical algorithms that significantly expand the scope of free probability in modern data applications.\n\n---\n\nWe sincerely thank the reviewer for their time and feedback, and would be happy to provide any additional clarification they might find helpful."}]

Author rebuttal:

Meta-review: This paper studies the problem of estimating the spectrum of a very large "impalpable" matrix from a small submatrix. They introduce a new method called free decompression to perform this task, based on tools from free probability. 

The ideas are novel and interesting, backed up by theoretical foundations, and experiments are performed to validate the efficiency of the proposed algorithm. Many in the NeurIPS audience may not have the background to appreciate the theoretical foundations of the paper, but many would appreciate the problem and woud be interested in learning about these new tools.
Decision: Accept (spotlight)