Spectral Estimation with Free Decompression

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.

---

Clarity and Notation Improvements

We thank the reviewer for their comments regarding notation and presentation. To address these concerns, we have added a new Appendix A that includes:

• 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).
• 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.

In 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.

Specifically:

• 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.
• 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.

Regarding 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.

We hope these revisions improve the clarity and accessibility of the manuscript for a broader readership, including those less familiar with RMT.

---

Clarifying Theoretical and Practical Contributions

We 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.

Our theoretical contributions go beyond a direct application of known results. In particular:

• 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.
• 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.
• 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.
• 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.
• 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.

Beyond 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.

We 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.

---

We sincerely thank the reviewer for their time and feedback, and would be happy to provide any additional clarification they might find helpful.

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.

---

Computing Log-Determinants

To 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.

At 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.

While 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.

---

Impact and Applications

The 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.

Matrix Functionals

Free decompression (FD) approximates spectral quantities of the form $f(\mathbf{A}_n) = n \int f(\lambda) \rho(\lambda) \mathrm{d}\lambda,$ where $\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:

• Gaussian-process marginal likelihood and predictive variance,
• Determinantal Point Process (DPP) training and sampling,
• Laplace-approximated Bayesian inference (e.g., for GLMs and neural networks),
• PAC-Bayes generalization bounds in deep learning,
• Kalman filtering and EM-based state-space models.

This category also includes spectral entropy measures such as the Von Neumann entropy, $S_{\mathrm{VN}} = - n \int \lambda \log (\lambda) \rho(\lambda) \mathrm{d}\lambda,$ which 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.

Quadratic Forms

FD can also be used to approximate quadratic forms involving (powers of) inverse matrices: $\boldsymbol{y}^{\intercal} (\mathbf{A}_n + \sigma^2 \mathbf{I})^{-1} \boldsymbol{y}.$ While 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.

Condition Number Estimation

The 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.

---

We 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.

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.

---

On the Scope and Comparison with Alternatives

Regarding 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.

Since 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.

In 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.

---

Clarification on First-Order Approximation

To analyze the stability of the nonlinear PDE, we consider a perturbation of the solution of the form $\tilde{m}(t, z) = m(t, z) + \delta m(t, z),$ where $\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.

Importantly, 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.

---

Clarification on Boundedness of $C_y(t)$ in Proposition E.1

We 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.

Recall from equation (E.3) that $C_y(t)$ is composed of two terms: $\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|.$ We 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:

• Since $\rho(t,x)$ is supported in a compact interval $[\lambda_{-}, \lambda_{+}]$ of diameter $D = \lambda_{+} - \lambda_{-}$, we have $|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}.$
• It follows that $|m(t,z)^{-1}| \leq \frac{D^2 + y^2}{y}$, so $|\Im m(t,z)^{-1}|$ also has this bound.
• For the derivative, note $\partial_z m(t,z) = \int_{\mathbb{R}} \frac{\rho(t,\xi)}{(\xi - z)^2} \mathrm{d}\xi,$ so $|\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}.$
• Combining the above, we obtain $\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},$ which is finite for every $y > 0$.

This confirms that $C_y(t) < \infty$ for all $y > 0$. We will be happy to include this derivation in the revised appendix for completeness.

Note 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.

---

Sampling Considerations

The 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.

Any 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.

From 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.

We 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.

---

Regarding a possible lower bound

A 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$.

Any 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.

Non-Hermitian Matrices

Our 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.

---

We 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.

We thank the reviewer for their thoughtful review and constructive feedback. Below we address the specific questions and concerns raised.

---

Clarification on Freeness Assumptions

We 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.

Nica (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.

---

Analytic continuation and the Stieltjes inversion formula

The 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).

However, 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)$.

In 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.

More 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.

We have added these clarifications to the text.

---

Outlier eigenvalues and low-rank structures

We 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.

Our 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.

However, 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.

While 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.

---

We 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.

---

Reference

[1] Nica, A. (1993). Asymptotically free families of random unitaries in symmetric groups. Pacific Journal of Mathematics, 157(2).