{'fraction_non_alphanumeric': 0.040150298475803765, 'fraction_numerical': 0.0166358595194085, 'mean_word_length': 4.342723004694836, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Estimating the Spectral Density of Large Implicit Matrices

Many important problems are characterized by the eigenvalues of a large matrix. For example, the difficulty of many optimization problems, such as those arising from the fitting of large models in statistics and machine learning, can be investigated via the spectrum of the Hessian of the empirical loss function. Network data can be understood via the eigenstructure of a graph Laplacian matrix using spectral graph theory. Quantum simulations and other many-body problems are often characterized via the eigenvalues of the solution space, as are various dynamic systems. However, naive eigenvalue estimation is computationally expensive even when the matrix can be represented; in many of these situations the matrix is so large as to only be available implicitly via products with vectors. Even worse, one may only have noisy estimates of such matrix vector products. In this work, we combine several different techniques for randomized estimation and show that it is possible to construct unbiased estimators to answer a broad class of questions about the spectra of such implicit matrices, even in the presence of noise. We validate these methods on large-scale problems in which graph theory and random matrix theory provide ground truth.


Introduction
The behaviors of scientific phenomena and complex engineered systems are often characterized via matrices. The entries of these matrices may reflect interactions between variables, edges between vertices in a graph, or the local geometry of a complicated surface. When the systems under consideration become large, these matrices also become large and it rapidly becomes uninformative -and potentially infeasible -to investigate the individual entries. Instead, we frame our questions about these matrices in terms of aggregate properties that are insensitive to irrelevant details.
To this end, the eigenvalues and eigenvectors of these matrices often provide a succinct way to identify insightful global structure of complex problems. The eigenvalues of the Hessian matrix in an optimization problem help us understand whether the problem is locally convex or has saddle points. The influential PageRank algorithm (Brin and Page, 1998) examines the principal eigenvector of a graph for the most "important" vertices. The sum of the eigenvalues of an adjacency matrix is closely related to the number of triangles in a graph, and insights from the eigenvalues of the Laplacian matrix have given rise to the field of spectral graph theory (see, e.g., Chung (1997)). In dynamical systems, numerical algorithms, and control, eigenstructure helps us analyze local stability and typical behaviors (Luenberger, 1979;Bertsekas, 2016).
The problem of estimating the eigenvalue distribution of large matrices has a particularly long history in fields such as condensed matter physics and quantum chemistry, where one often wishes to estimate the density of states (DOS). The stochastic methods we develop in this paper build upon these decades of work, which we briefly review here. One of the core components of our procedure is the use of a stochastic estimator for the trace of a matrix by constructing a quadratic form with a random vector. This procedure is generally attributed to Hutchinson (1989), but it was coincidentally introduced in Skilling (1989). This latter work by Skilling sketches out an approach to eigenvalue estimation based on generalized traces and Chebyshev polynomials that strongly influences the present paper. Although Lanczos and Chebyshev polynomial expansions had been used for many years in physics and quantum chemistry (see, e.g., Haydock et al. (1972), Tal-Ezer and Kosloff (1984)), the combination of polynomial expansion and randomized trace estimation as in Skilling (1989) led to the kernel polynomial method for estimation the density of states in Silver and Röder (1994) and Silver et al. (1996), in which a smoothing kernel is convolved with the spectral density in order to reduce the Gibbs phenomenon. See Weiße et al. (2006) for a review of the kernel polynomial method and a broader historical perspective on these estimation techniques in physics and quantum chemistry. Outside

Problem Formulation
The matrix of interest A is assumed to be real, square, symmetric, and and of dimension D, i.e., A ∈ R D×D where A ij = A ji . Many of the results in this paper generalize to complex and asymmetric matrices, but the real symmetric case is of primary interest for Hessians and undirected graphs. The primary challenge to be addressed is that D can be very large. For example, the modern stochastic optimization problems arising from the fitting of neural networks may have tens of billions of variables (Coates et al., 2013) and the graph of web pages has on the order of a trillion vertices (Alpert and Hajaj, 2008). When A is this large, it cannot be explicitly represented as a matrix, but can only be indirectly accessed via matrix-vector products of the form Ax, where we can choose x ∈ R D . An additional challenge arises in the case of large scale model fitting, where the loss surface may be a sum of billions of terms, and so even Ax cannot be feasibly computed. Instead, an unbiased but randomized estimatorÂx is formed from a "mini-batch" of data, as in stochastic optimization (Robbins and Monro, 1951).
We will assume that all of the eigenvalues of A are within the interval (−1, 1), something achievable by dividing A by its operator norm plus a small constant. With this assumption, the abstraction we choose for investigating A is the generalized trace tr(F (A)), where F : R D×D → R D×D is an operator function that generalizes a scalar function f : (−1, 1) → R to square matrices. We take f to have a Chebyshev expansion γ 0 , γ 1 , . . .: where the functions T k : (−1, 1) → [−1, 1] are the Chebyshev polynomials of the first kind. These are defined via the recursion Chebyshev approximations are appealing as they are optimal under the L ∞ norm. Other polynomial expansions are optimal under different function norms, but this generalization is not considered here. We can generalize the Chebyshev polynomials to square matrices using the recursion Proposition 2.1. Let A be diagonalizable into A = U T ΛU with orthonormal U and diagonal Λ. Then the matrix Chebyshev polynomial T k (A) applies the scalar Chebyshev polynomial T k (a) to each of the eigenvalues of A.
Proof. The proposition is true by inspection for T 0 (A) and T 1 (A). We proceed inductively using the two-term recursion. Assume that the proposition is true for T k−1 (A) and T k−2 (A), then Applying the Chebyshev recursion: When f (a) has Chebyshev coefficients γ 1 , γ 2 , . . ., we thus define F (A) as (2.10) This matrix function can now be seen to be simply applying f (a) to the eigenvalues of A. Returning to the generalized trace, we see that it corresponds to first applying f to each eigenvalue and then summing them: where the λ d (A) are the eigenvalues of A. By computing the traces of the Chebyshev polynomials, we can then can choose different f to ask different questions about the spectra of A. In the proceeding sections, we will describe the construction of estimatorsŜ such that E[Ŝ] = tr(F (A)). The generalized trace is a surprisingly flexible tool for interrogating the spectrum. In particular, it enables us to look at the eigenvalues via the normalized spectral measure, which is a sum of Dirac measures: (2.12) A generalized trace tr(F (A)) can thus be seen as an expectation under this spectral measure: Through the spectral measure, one can compute empirical summaries of the matrix of interest, such as How many eigenvalues are within of zero? or How many negative eigenvalues are there?. Such questions can be viewed as particular instances of functions f where they are, e.g., indicator-like functions.
In this paper, we will take this general view and seek to construct a general purpose spectral density estimate from the Chebyshev traces of A using a smoothing kernel K(λ, λ ). If K(λ, λ ) ≥ 0 and 1 −1 K(λ, λ ) dλ = 1, then we can define the smoothed spectral densitỹ (2.14) In the generalized trace framework, we can ask pointwise questions aboutψ by constructing a λ-specific f λ (a) = K(λ, a).
We are seeking to estimate a one-dimensional density, so a dense collection of pointwise estimates is useful for visualization, diagnosis, or to ask various post hoc questions such as those above. As in other kinds of density estimation, smoothing is sensible as it leads to practical numerics and we typically do not care about high sensitivity to small perturbations of the eigenvalues. Smoothing also minimizes the effects of the Gibbs phenomenon when using Chebyshev polynomials. Various smoothing kernels have been proposed for the case where the polynomials are truncated (See Weiße et al. (2006) for a discussion.), but in this work we will not have a fixed-order truncation and so we will derive a kernel that is 1) appropriate for density estimation, 2) reduces the Gibbs phenomenon, 3) maps naturally to (−1, 1), and 4) does not depend on the order of the approximation. We will break down the task of computing estimates of generalized traces into two parts: the application of F to the matrix of interest (Section 3), and the computation of the trace of the result (Section 4). In both cases, we will use randomized algorithms to form unbiased estimates.

Randomized Matrix Power Series Estimation
In this section we discuss the construction of unbiased Monte Carlo estimates of functions applied to symmetric matrices. Two forms of Monte Carlo estimator will be used: 1) multiplications of independent random matrix variates in order to get unbiased estimates of matrix powers, and 2) Monte Carlo estimation of infinite polynomial series. These are framed in terms of randomized Chebyshev polynomials, slightly generalizing the scalar-valued power series importance sampling estimator described in Papaspiliopoulos (2011, Section 4.6). 4 The core observation in the randomized estimation of infinite series is that sums such as Equations (2.1) and (2.10) can be turned into expectations by introducing a proposal distribution that has support on the non-negative natural numbers N 0 . The estimator is then a weighted random power or polynomial of random order. We are here making the further assumption that we only have randomized estimates of the matrix and so these matrix powers must be constructed with care: the power of an expectation is not the expectation of the power. Instead, we use a product of independent estimates of the matrix, with the kth independent estimate denotedÂ k , to construct an unbiased estimate of the matrix power. For example, in the case of fitting a large statistical or machine learning regression model, each Hessian estimateÂ k might arise from small and independent subsets of data. This relieves us of the burden of computing a Hessian that is a sum over millions or billions of data. The following lemma addresses such an unbiased construction for Chebyshev polynomials.
Lemma 3.1. LetÂ 1 ,Â 2 , . . . be a sequence of independent, unbiased estimators of a symmetric real matrix A ∈ R D×D . From this sequence, construct a new sequence as follows: (3.1) ThenT k is an unbiased estimate of the matrix-valued Chebyshev polynomial of the first kind T k (A).
Proof. Unbiasedness for k ∈ {0, 1} is clear by inspection. For k > 1, we proceed inductively by assuming thatT k−1 andT k−2 are unbiased.Â k is independent ofT k−1 and therefore Although this randomized estimator is unbiased, its variance explodes as a function of k. The intuitive reason for this is that deterministic Chebyshev polynomials have a somewhat magical property: even though the leading coefficients increase exponentially with k, the terms cancel out in such a way that the sum is still within the interval [−1, 1]. However, in the case here the polynomial is random. The noise overwhelms the careful interaction between terms and results in an explosion of variance. However, it is possible to find an upper bound for this variance and identify situations where the Chebyshev coefficients decay fast enough to achieve finite variance. For a D × D matrix A, we denote its spectral norm by A and its Frobenius norm Lemma 3.2. LetÂ 1 ,Â 2 , . . . be a sequence of independent, unbiased estimators of a symmetric real matrix A ∈ R D×D and letT k be defined as in (3.1). Suppose that Then, for all k ≥ 0, and note that tr(B k ) = t k + t k−1 for all k ≥ 1 and that B k is positive semidefinite for all k ≥ 1. We can write the defining recurrence forT k as This gives the following recurrence for B k : Taking the trace of both sides gives where we have use the inequality tr(AB) ≤ A tr(B) for positive semidefinite matrices A and B, together with the triangle inequality, and the definition of α. We then see that for all k ≥ 2, it follows that t k ≤ Dα k for all k ≥ 0.
Fortunately, in the situation of interest here, the coefficients weighting each of the Chebyshev terms fall off even more rapidly than the variance increases. In particular, the von Mises smoothing introduced in Section 3.1 has subexponential tails that counteract the explosion of the noisy Chebyshev polynomials. Thus with independent estimators to form stable unbiased estimates of the weighted Chebyshev polynomials, we can now form an unbiased estimate of F (A).
Proposition 3.1. Let F : R D×D → R D×D be an analytic operator function with Chebyshev polynomial coefficients γ 0 , γ 1 , . . .. Construct the sequenceT 0 ,T 1 , . . . as in Lemma 3.1 above. Let π be a probability distribution on the natural numbers (including zero) N 0 , where π k is the probability mass function and π k > 0 if γ k = 0. Construct an estimatorF by first drawing a random nonnegative integer k from π and then Proof. This simply combines importance sampling with the Chebyshev estimator from Lemma 3.1: The variance of importance sampling estimates can be poor. In the following, we derive the proposal distribution that minimizes the variance of the estimator under the Frobenius norm, i.e., it minimizes the expected squared Frobenius norm of the difference between the estimatorF and its mean F (A). 6 Proposition 3.2. Suppose the sequence γ 0 , γ 1 , . . . is such that The Frobenius norm can be written as a trace, so that The second term is independent of π, so it suffices to minimize the first term. Under π this term can be expanded as To show that this value is minimal, we consider an alternative distribution π and apply Jensen's inequality on the convex function z → z 2 to write Thus we can see that π offers expected squared Frobenius norm less than or equal to that expected under all other possible proposals π.
We also note that an importance sampler can be constructed with a hazard rate formulation, resulting in a "randomized truncation" that is still unbiased, using an estimator with k ∼ π as before: (3.12) Although this is appealing because it makes use of the intermediate steps of the recursive computation ofT k , we have not been able to identify hazard rates that generally result in variance reduction. It nevertheless seems likely that such hazard rates exist.

Von Mises Smoothing
As discussed in Section 2, it is useful to frame the problem in terms of estimating the spectral measure ψ(λ). This is closely related to the kernel polynomial method (see, e.g., Weiße et al. (2006)) in which truncated Chebyshev approximations are smoothed with order-dependent Fejér, Jackson, or Lorentz kernels to minimize the Gibbs phenomenon in the approximation. Here we are not making fixed-order truncations, but smoothing is still desirable. We introduce the von Mises kernel to address our situation. It has support on 7 the interval (−1, 1), integrates to one, and is approximately Gaussian for large κ, but with a convenient closed-from Chebyshev expansion: where I η (κ) is the modified Bessel function of the first kind. This kernel arises from constructing a von Mises distribution on the unit circle and then projecting it down onto the interval (−1, 1), along with the appropriate Jacobian adjustment. The parameter κ is inversely related to the squared width of the kernel. Figure 1 shows several densities from different parameters, on the unit circle and projected.
Of particular interest is that the von Mises kernel has a closed-form Chebyshev polynomial expansion for any λ . Chebyshev polynomials can be constructed in three different ways. First, on the interval (−1, 1) they can be constructed via the recursion in Section 2 or as T k (x) = cos(k cos −1 (x)). They can also be computed on the unit circle in the complex plane as T k (x) = 1 2 (z k + z −k ) where z = cos −1 (x) + i sin −1 (x). Finally, they can be constructed via the phase alone, as |z| = 1 on the unit circle, i.e., z = e iπθ , so that T k (x) = cos(kθ). In the z and θ cases, the polynomials are projected onto the interval (−1, 1). See Trefethen (2013) for a comprehensive treatment of these representations.
We focus only on the phase view. The von Mises distribution is a probability distribution defined on the unit circle, with density given by (3.14) Its variance is given by 1 − I 1 (κ)/I 0 (κ). For a fixed µ and κ, the coefficients can be computed via (3.16) The first coefficient for Chebyshev polynomials does not have the factor of two and so γ 0 = 1/π. Returning 8 to the interval (−1, 1) we have For our purposes here, κ ≈ 5000 and so the ratio of modified Bessel functions is well approximated as: This approximation is computationally useful and also makes it clear that the γ k rapidly converge to zero as a function of k.

Bounding the Variance of the Importance Samples
When computing an importance sampled estimate of the von Mises kernel, a natural choice of proposal distribution is to have it be proportional to the absolute value of the coefficients in Equations (3.17) and (3.18).
Here the proposal distribution will be denoted as q rather than π to distinguish it from the constant π: where the normalizing constant is given by (3.22) Note that the 1/π factors are absorbed into Z(λ , κ). Observing that we can see that the normalization constant has upper bound In the theorem below we use this bound on Z(λ , κ) and the bound on the Chebyshev estimators from Lemma 3.2 to provide an upper bound on the expected squared Frobenius norm of the overall importance sampler, as a function of λ and κ.
Proof. We first use the law of total expectation and incorporate the bound from Lemma 3.2: The γ k and proposal are chosen to resolve to a simple form: (3.27) In particular, we get a straightforward bound on this term for k > 0 by discarding the cosine: (3.28)

Substituting this back into the overall bound and using the standard result
We can now insert the bound for Z(λ , κ) to get the overall result: This bound does not depend on k, but we can see that the variance effectively grows exponentially with κ. This can be seen directly as a bias/variance tradeoff: larger κ corresponds to less smoothing and lower bias with a large variance penalty.

Skilling-Hutchinson Randomized Estimation of Generalized Traces
Having constructed an unbiased estimator of analytic operator functions, we now focus on estimation of the trace of such matrices. The classic randomized estimator described in Hutchinson (1989) and Skilling (1989), is the following, which we will call the Skilling-Hutchinson trace estimator : Proof.
trace of a scalar

Algorithm 1 Spectral Density Estimation (Without Control Variates)
Ensure: Smoothing parameter κ, query location λ, sequenceÂ 1 ,Â 2 , . . . Require: Unbiased estimate of the κ-smoothed spectral densityψ (λ) 1: π ← [|γ 0 (λ, κ)|, |γ 1 (λ, κ)|, |γ 2 (λ, κ)|, . . .] / ∞ k=0 |γ k (λ, κ)| 2: k ∼ π 3: x ∼ N (0, I) The requirement for the random variable x is fairly weak and is satisfied by, for example, standard normal and Rademacher variates. The properties of such estimators have been studied in Avron and Toledo (2011). Finally, the following theorem makes it clear that this trace estimator can be combined with the randomized Chebyshev expansion to construct an overall unbiased estimate of the generalized trace. Proof. Following the simple argument from the proof of Proposition 4.1, we have The resulting basic algorithm that combines randomized Chebyshev polynomials, von Mises smoothing, and the Skilling-Hutchinson trace estimator is shown in Algorithm 1. For simplicity, the pseudocode shows a single sample for both the polynomial and the quadratic form.

Variance Reduction Through a Control Variate
The Skilling-Hutchinson trace estimator's variance can be reduced via the construction of a control variate. Let x ∼ N (0, I) and form the Skilling-Hutchinson trace estimator asẐ = x T Ax. A natural choice of control variate is to compute the quadratic form using a second matrix B which has a known trace. This is similar to a preconditioner in which one introduces a second related matrix that is somehow easier to deal with than the original. In this case, the new estimator iŝ Proposition 4.2. The variance of the Skilling-Hutchinson trace estimator applied to a real symmetric A is

Algorithm 2 Spectral Density Estimation (With Control Variate)
Ensure: Smoothing param κ, query location λ, sequenceÂ 1 ,Â 2 , . . ., matrix B, weighting c Require: Unbiased estimate of the κ-smoothed spectral densityψ (λ) 1: π ← [|γ 0 (λ, κ)|, |γ 1 (λ, κ)|, |γ 2 (λ, κ)|, . . .] / ∞ k=0 |γ k (λ, κ)| 2: k ∼ π 3: x ∼ N (0, I) Proof. The general expression for the variance follows from (4.6) When x ∼ N (0, I), we can use the formulas for the second and fourth moments, resulting in the simplification An alternative proof of the identity for the Gaussian case, from Avron and Toledo (2011), uses the rotational symmetry of the Gaussian. Decompose the matrix into A = U T ΛU, where U is orthogonal and Λ is diagonal. As x is a spherical Gaussian, y = Ux is also marginally spherical Gaussian. Thus (4.10) The random variates y 2 d are independently χ 2 -distributed with one degree of freedom and so have variance 2. The variance of the overall estimator is then given by twice the sum of the squared eigenvalues, which is equal to twice the squared Frobenius norm of A: Proof. For a estimatorẐ formed fromẐ and a zero-mean control variateŴ , a classic result from Monte Carlo is that the variance-minimizing weighting is given by c = C[Ẑ,Ŵ ]/V[Ŵ ], where C[Ẑ,Ŵ ] is the covariance between the original estimator and the control-variate. In the specific case of interest here, then: where C(·, ·) is the covariance between the estimator and the control variate under the distribution induced by x. The denominator is 2 tr(B 2 ) using the previous lemma. The covariance in the numerator can be computed via For a given B and weighting c, some algebra shows that the variance is reduced by When c is chosen optimally, the control variate reduces the variance by 2 tr(AB) 2 /tr(B 2 ). This is maximized in the trivial case when B = A. The first case of practical interest is B = αI for scalar α = 0. This reduces the variance by an amount that does not depend on α, as the optimal c scales the control variate appropriately. The second case of practical interest is where B is a diagonal matrix and, unsurprisingly, the variance minimizing diagonal matrix is the one that is the diagonal of A. If the Skilling-Hutchinson trace estimator is used sequentially, then the raw diagonal estimator sitting beneath it can be used to form a diagonal B, i.e., E[x T A x] = diag(A) where is element-wise product (Bekas et al., 2007). That is, the diagonal estimate from previous estimations can be used to form the control variate for later estimations. We also note that the control variate described above can be easily adapted for variance reduction in the diagonal estimator. When estimating B to be the diagonal of A, we set c = 1. Pseudocode is provided in Algorithm 2.

Empirical Evaluation of Estimator Performance
In this section we evaluate the performance of the estimator and demonstrate the feasibility of approximating the spectra of large matrices. We focus on two classes of matrices: sparse adjacency matrices of certain undirected graphs known as Kneser graphs, and dense random matrices designed to model the Hessian matrices often appearing in machine learning problems.

Kneser Graphs
Graph theory provides a natural source of large matrices with well-understood spectra. The Kneser Graph, K(n, k) is constructed from a set S with |S| = n. The estimated and true spectrum of the Kneser graph K(23, 11). This estimate was generated with κ = 1000.
The graph K(n, k) is regular, and each vertex has degree n−k k . The eigenvalues of K(n, k) are for i = 0, ..., k. For a more detailed discussion of Kneser graphs, see Godsil and Royle (2001).
Choosing n = 23 and k = 11, we get a regular graph with 1352078 and 16224936 edges and degree 12. This graph has the following eigenvalues, with associated multiplicities: Eigenvalue Multiplicity  -11  22  -9  1518  -7  24794  -5  144210  -3  326876  -1  208012  2  326876  4  245157  6  67298  8  7084  10  230  12  1   14 To better simulate the situation where we have an unbiased estimator with non-zero variance, we added Gaussian noise with variance 1/d 2 to each non-zero entry of the adjacency matrix for K(23, 11). This noise was added independently each time matrix multiplication was applied. The smoothed spectral density was estimated as described in Algorithm 1. The one exception is that, for each random vector x, we select 10000 indices k ∼ π. We obtained excellent agreement between our predictions and the known eigenspectrum of K(23, 11), as shown in Figure 2.

Random Matrix Theory
A random matrix A is a matrix whose entries A ij are random variables. The distribution of the A ij determines the type of random matrix ensemble to which A belongs. We will focus on two ensembles: the real Wigner ensemble for which A is symmetric but otherwise the A ij are independent with mean zero; and the real Wishart ensemble for which A = BB T , where the columns of B are independent and normally distributed with mean zero. We will also focus on the limiting spectral density of A, i.e., the large D limit of the empirical spectral distribution defined in Equation. 2.12, ( 5.1) Here we attach the D subscript to ψ to indicate the size of the matrix when it is relevant, with ψ ∞ (λ) being the limiting spectral density. For the random matrix ensembles considered here, and under mild conditions (such as bounded fourth moments), ψ D (λ) converges almost surely to ψ ∞ (λ). Throughout this section, validation calculations were performed with noisy versions of the Wishart and Wigner ensembles. Multiplicative Gaussian noise with variance 100/d 2 was applied independently at each matrix multiplication step, where d is the dimension of the matrix in question. As in Subsection 5.1, spectral density was estimated as described in Algorithm 1. This time, we select 100000 indices k ∼ π for each of 2048 random vectors x.

Wigner Ensemble
For concreteness, we consider a real Wigner matrix A whose elements are normally distributed, The limiting spectral density of a random matrix of this type was first computed by Wigner (1955). The result is the celebrated semi-circle law, Because ψ sc is an even function, its odd moments vanish. The even moments are related to the Catalan numbers, C k , From Stirling's formula, it is easy to see that for large k, the moments decrease as the 3/2 power, m 2k = 1 √ πk 3/2 + O(k −5/2 ) , (5.5) but they are nonzero for any finite k. In contrast, the Chebyshev expansion actually truncates at finite order. Indeed, by inspection we have, This result holds only as D → ∞; there will be finite-size corrections to this result for finite D. Moreover, for our simulations, there is additional error resulting from the variance of our estimators; nevertheless, we still observe excellent agreement with these infinite-D limiting predictions. See Figure 3.

Wishart Ensemble
Next we consider matrices A drawn from the Wishart ensemble. In particular, we consider A = BB T where B ∈ R D×N and B ij ∼ N (0, σ 2 N ). If both N and D are taken to infinity at the same rate, i.e., D, N → ∞ with D = φN for some constant φ, then the limiting spectral density can be shown to obey the Marchenko-Pastur distribution (Marčenko and Pastur, 1967), where, The Wishart matrix A is low rank if φ > 1, and the nullity equals (N − D)/N = 1 − φ, which is the source of the delta function density at 0 in that case. Note that the minimum eigenvalue is λ − , which may be larger than 0, and depends on the value of φ. The moments of the Marchnko-Pastur distribution are polynomials in φ with coefficients equal to the Narayana numbers, (5.10) where 2 F 1 is the Gauss hypergeometric function and N (j, k) are the Narayana numbers, The estimated and true spectral densities of six matrices are shown. These matrices are convex sums of Wishart and Wigner matrices, so their spectra are known in closed form to provide ground truth for the randomized estimation. Estimates were produced using 500 samples with κ = 5000. Although some estimates exhibit variance due to boundary effects, the overall differences are small.
Although relatively complicated, the Chebyshev coefficients can be written as a sum over the m k , (5.12) One reason that this expansion is more complicated than that of the Wigner case is because the support of the spectrum is in the range [λ − , λ + ], whereas the Chebyshev polynomials are a natural basis on [−1, 1].
Indeed, a much simpler expansion arises when A is scaled and shifted so that its spectrum has support on [−1, 1]. For φ < 1, lettingÃ = αA + βI D with α = 1 2σ 2 √ φ and β = − 1+φ 2 √ φ gives a spectral densityψ mp (λ), The Chebyshev coefficients can be found to equal, (5.14) Again, these results hold only as D → ∞, and again there is still additional variance from our estimators, but we still observe excellent agreement with these infinite-D limiting predictions. See Figure 3.

Wishart plus Wigner
As described by Pennington and Bahri (2017), the sum of a Wishart and a Wigner matrix can serve as a simplified model of the Hessian matrix of the loss function of machine learning models, as we now review. Suppose we consider an arbitrary (possibly non-convex) parameterizable model function f (θ) ≡ f θ (x) of some data x. If the model is evaluated with a convex loss function l(f (θ)), then the Hessian matrix has a simple form. In particular, the second derivative of the loss with respect to the parameters θ is, (5.15) where, [H l ] ij = J iα J jβ ∂ 2 l(f (θ)) ∂f α (θ)∂f β (θ) , and [H f ] ij = ∂l(f (θ)) ∂f α (θ) (5.16) and the Jacobian matrix is, Because we've assumed l to be convex, H l is PSD. A simple approximation which captures this additive structure of the Hessian is outlined in Pennington and Bahri (2017), and involves approximating H l as a Wishart matrix and H f as a Wigner matrix. In particular, we will consider the approximation, where B is standard Wigner and C is standard Wishart, as defined above, and γ is an interpolation parameter that can be interpreted as a measure of the loss value at convergence. When γ = 1, the model has converged to a high-index saddle point with a large loss value, at which point the Hessian has no contribution from the Wishart matrix C and has an equal number of positive and negative eigenvalues. The global minimum corresponds to γ = 0, at which point there is no contribution from the Wigner matrix B and the Hessian has no negative eigenvalues. As described in Pennington and Bahri (2017), the number of negative eigenvalues can be computed as a function of ≡ γ 2 2(1−γ) 2 , which measures the relative contribution of B and C to H. The details are not particularly illuminating, but the result is, where ξ ± = ξ 2 ± χ 2 and, ξ = 2 −1/6 χ 1 + χ 2 1 − 2χ 3 2 1/3 , χ 1 = 4 3 + 9 2 φ + 18 2 φ 2 , χ 2 = 2 2 + 3 φ − 3 φ 2 .

Discussion
In this work we have described a framework for estimation of spectral densities of large matrices. Our primary objective has been to construct a tool for the empirical investigation of matrices that can only be probed in implicit and noisy ways. Our overall approach has been to combine ideas that have been proposed across several different fields into an estimation framework that provides unbiased estimates with controlled variance, even when the matrices become large. We use this framework to then construct smoothed estimates of the overall density. This smoothing necessarily introduces bias, but it allows us to visualize spectra whose atomic nature would otherwise make them difficult to interrogate. In future work, we intend to use this tool in the analysis of high dimensional statistical and machine learning models, where the local geometry of the optimization surface-as captured by the Hessian and Fisher matrices-is challenging to examine.