## 2 Spectral Density Estimation  
The _cumulative spectral density function_ or _empirical spectral cumulative distribution_ of the matrix \(\mathbf{A}\) is defined as  
\[P_{\lambda}(z):=\frac{1}{N}\sum_{\ell=1}^{N}\mathds{1}_{\{\lambda_{\ell}\leq z \}}, \tag{5}\]  
where \(\mathds{1}_{\{C\}}=1\) if statement \(C\) is true and \(\mathds{1}_{\{C\}}=0\) otherwise. The _spectral density function_[36] (Chapter 6) (also called the _Density of States_ or _empirical spectral distribution_[37] (Chapter 2.4)) of \(\mathbf{A}\) is the probability measure defined as  
\[p_{\lambda}(z):=\frac{1}{N}\sum_{\ell=1}^{N}\mathds{1}_{\{\lambda_{\ell}=z\}}.\]  
[PERSON] et al. [9] provide an overview of methods to approximate these functions. In this work, we use a variant of the Kernel Polynomial Method (KPM) [38; 39; 40] described in [9; 41] to estimate the cumulative spectral density function \(P_{\lambda}(z)\) of \(\mathbf{A}\). Namely, for each of \(S\) linearly spaced points \(\{\tilde{\zeta}_{i}\}_{i=1}^{S}\) between \(\underline{\lambda}\) and \(\overline{\lambda}\), we estimate the number of eigenvalues less than or equal to \(\xi_{i}\) via stochastic trace estimation [42; 43]. Let \(\mathbf{x}\) denote a Gaussian random vector with distribution  
Figure 1: Degree 5 polynomial approximations of the function \(f(\lambda)=e^{-\lambda}\) of the graph Laplacian of a random Erdős–Renyi graph with 500 vertices and edge probability 0.2. The discrete least squares approximation incurs larger errors in the lower end of the spectrum. However, since the eigenvalues are concentrated at the upper end of the spectrum, it yields a lower approximation error \(||f(\mathbf{A})-p_{5}(\mathbf{A})||_{2}\).  
denote a sample of size \(J\) from this distribution, and \(\tilde{\Theta}_{\xi_{i}}\) denote a Jackson-Chebyshev polynomial approximation to \(\Theta_{\xi_{i}}(\lambda):=\openone_{\{\lambda\leq\xi_{i}\}}\)([PERSON], 1966; [PERSON], 1966). The stochastic trace estimate of the number of eigenvalues less than or equal to \(\xi_{i}\) is then given by  
\[\eta_{i}=\text{tr}\big{(}\Theta_{\xi_{i}}(\mathbf{A})\big{)}=\mathbb{E}[ \mathbf{x}^{\top}\Theta_{\xi_{i}}(\mathbf{A})\mathbf{x}]\approx\frac{1}{J} \sum_{j=1}^{J}\mathbf{x}^{(j)\top}\tilde{\Theta}_{\xi_{i}}(\mathbf{A})\mathbf{ x}^{(j)}. \tag{6}\]  
As in ([PERSON], 1966), we then form an approximation \(\tilde{P}_{\lambda}(z)\) to \(P_{\lambda}(z)\) by performing monotonic piecewise cubic interpolation ([PERSON], 1966) on the series of points \(\big{\{}\big{(}\xi_{i},\frac{\eta_{i}}{N}\big{)}\big{\}}_{i=1,2,\ldots,S}\). Analytically differentiating \(\tilde{P}_{\lambda}(z)\) yields an approximation \(\tilde{p}_{\lambda}(z)\) to the spectral density function \(p_{\lambda}(z)\). Since \(\tilde{P}_{\lambda}(z)\) is a monotonic cubic spline, we can also analytically compute its inverse function \(\tilde{P}_{\lambda}^{-1}(y)\). The spectrum-adapted methods we propose in Section 3 utilize both \(p_{\lambda}(z)\) and \(\tilde{P}_{\lambda}^{-1}(y)\) to focus on regions of the spectrum with higher eigenvalue density when generating polynomial approximations. Figure 2 shows examples of the estimated cumulative spectral density functions for eight real, symmetric matrices \(\mathbf{A}\): the graph Laplacians of the Erdos-Renyi graph (gnp) from Figure 1, the Minnesota traffic network ([PERSON], 1966) (\(N=2642\)), and the Stanford bunny graph (S