## 2 Spectral Density Estimation
### Spectrum-Adapted Polynomial Regression/Orthogonal Polynomial Expansion  
A second approach is to solve the weighted least squares polynomial regression problem  
\[\min_{p\in\mathcal{P}_{K}}\sum_{m=1}^{M}w_{m}\left[f(x_{m})-p(x_{m})\right]^{2}, \tag{7}\]  
where the abscissae \(\{x_{m}\}_{m=1,2,\ldots,M}\) and weights \(\{w_{m}\}_{m=1,2,\ldots,M}\) are chosen to capture the estimated spectral density function. We investigated several methods to set the points (e.g., linearly spaced points, Chebyshev points on the interval \([\underline{\lambda},\overline{\lambda}]\), Chebyshev points on each subinterval \([\xi_{i},\xi_{i+1}]\), and warped points via the inverse of the estimated cumulative spectral density function as in Section 3.1) and weights (e.g., the analytically computed estimate \(\tilde{p}_{\lambda}\) of the spectral density function, a discrete estimate of the spectral density function based on the eigenvalue counts in (6), the original KPM density of states method based on a truncated Chebyshev expansion [9] (Equation 3.11), or equal weights for warped points). Without going into extensive detail about these various options, we remark that choosing abscissae \(\{x_{m}\}\) that are very close to each other, which may occur when using points warped by the inverse of the estimated density function, may lead to numerical instabilities when solving the weighted least squares problem. In the numerical experiments, we use \(M\) evenly spaced points on the interval \([\underline{\lambda},\overline{\lambda}]\) (i.e., \(x_{m}=\frac{m-1}{M-1}(\overline{\lambda}-\underline{\lambda})+\underline{ \lambda}\)), and set the weights to be \(w_{m}=\tilde{p}_{\lambda}(x_{m})\). To solve (7), we use a customized variant of MATLAB's polyfit function that solves, again via QR decomposition, the normal equations of the weighted least squares problem:  
\[\boldsymbol{\Psi}_{\mathbf{x}}^{\top}\mathbf{W}\boldsymbol{\Psi}_{\mathbf{x}} \mathbf{c}=\boldsymbol{\Psi}_{\mathbf{x}}^{\top}\mathbf{W}\mathbf{y},\]  
where \(\boldsymbol{\Psi}_{\mathbf{x}}\) is the Vandermonde matrix associated with the points \(\{x_{m}\}\), \(\mathbf{W}\) is a diagonal matrix with diagonal elements equal to the weights \(\{w_{m}\}\), \(\mathbf{y}\) is a column vector with entries equal to \(\{f(x_{m})\}\), and \(\mathbf{c}\) is the vector of unknown polynomial coefficients. Once the coefficients \(\mathbf{c}\) of the optimal polynomial, \(p_{K}^{*}\) are attained, \(p_{K}^{*}(\mathbf{A})\mathbf{b  
An alternative way to view this weighted least squares method [52] is as a truncated expansion in polynomials orthogonal with respect to the discrete measure \(d\lambda_{M}\) with finite support at the points \(\{x_{m}\}\), and an associated inner product [53] (Section 1.1)  
\[\langle f,g\rangle_{d\lambda_{M}}=\int_{\mathbb{R}}f(x)g(x)d\lambda_{M}=\sum_{ m=1}^{M}w_{m}f(x_{m})g(x_{m}).\]  
The \(M\) discrete monic orthogonal polynomials \(\{\pi_{k,M}\}_{k=0,1,M-1}\) satisfy the three-term recurrence relation [53] (Section 1.3)  
\[\pi_{k+1,M}(x)=(x-\alpha_{k,M})\pi_{k,M}(x)-\beta_{k,M}\pi_{k-1,M}(x),\;k=0,1, \ldots,M-1, \tag{8}\]  
with \(\pi_{-1,M}(x)=0\), \(\pi_{0,M}(x)=1\), \(\beta_{0,M}=\sum_{m=1}^{M}w_{m}\),  
\[\alpha_{k,M}=\frac{\langle\pi\pi_{k,M},\pi_{k,M}\rangle_{d\lambda_{M}}}{ \langle\pi_{k,M},\pi_{k,M}\rangle_{d\lambda_{M}}},\;k=0,1,\ldots,M-1,\text{and }\beta_{k,M}=\frac{\langle\pi_{k,M},\pi_{k,M}\rangle_{d\lambda_{M}}}{\langle \pi_{k-1,M},\pi_{k-1,M}\rangle_{d\lambda_{M}}},\;k=1,2,\ldots,M-1.\]  
Given the abscissae \(\{x_{m}\}\) and weights \(\{w_{m}\}\), the three-term recursion coefficients \(\{\alpha_{k,M}\}_{k=0,1,\ldots,M-1}\) and \(\{\beta_{k,M}\}_{k=1,2,\ldots,M-1}\) can also be computed through a stable Lanczos type algorithm on an \((M+1)\times(M+1)\) matrix [53] (Section 2.2.3), [54].