## 2 Spectral Density Estimation
### Spectrum-Adapted Polynomial Regression/Orthogonal Polynomial Expansion
In matrix-vector notation, the vectors \(\pi_{k,M}\in\mathbb{R}^{M}\), which are the discrete orthogonal polynomials evaluated at the \(M\) abscissae, can be computed iteratively by the relation  
\[\boldsymbol{\pi}_{k+1,M}=(\text{diag}(\{x_{m}\})-\alpha_{k,M}\mathbf{I}_{M}) \boldsymbol{\pi}_{k,M}-\beta_{k,M}\boldsymbol{\pi}_{k-1,M},\;k=0,1,\ldots,M-1,\]  
with \(\boldsymbol{\pi}_{-1,M}=\boldsymbol{0}_{M}\) and \(\boldsymbol{\pi}_{0,M}=\boldsymbol{1}_{M}\). Figure 4 shows an example of these discrete orthogonal polynomials. Figure 4: Comparison of (top) the first six discrete orthogonal polynomials defined in (8), adapted to the estimated cumulative spectral density of the random Erdös–Renyi graph from Figures 1–3, to (bottom) the first six standard shifted Chebyshev polynomials with degree \(k=0\) to \(5\). In the region of high spectral density, shown in the zoomed boxes on the right, the discrete orthogonal polynomials feature more oscillation while maintaining small amplitudes, enabling better approximation of smooth functions in this region. Finally, the degree \(K\) polynomial approximation to \(f(\mathbf{A})\mathbf{b}\) is computed as  
\[p_{K}(\mathbf{A})\mathbf{b}=\sum_{k=0}^{K}\frac{\langle f,\pi_{k,M}\rangle_{d \lambda_{M}}}{\langle\pi_{k,M},\pi_{k,M}\rangle_{d\lambda_{M}}}\pi_{k,M}( \mathbf{A})\mathbf{b},\]  
with \(\pi_{-1,M}(\mathbf{A})\mathbf{b}=\mathbf{0}_{N}\), \(\pi_{0,M}(\mathbf{A})\mathbf{b}=\mathbf{b}\), and  
\[\pi_{k+1,M}(\mathbf{A})\mathbf{b}=(\mathbf{A}-\alpha_{k,M}\mathbf{I}_{N})\pi_ {k,M}(\mathbf{A})\mathbf{b}-\beta_{k,M}\pi_{k-1,M}(\mathbf{A})\mathbf{b},\ k=0,1, \ldots,K-1\ (\text{where }K\leq M-1).\]  
Before proceeding to numerical experiments, we briefly comment on the relationship between the spectrum-adapted approximation proposed in this section and the Lanczos approximation to \(f(\mathbf{A})\mathbf{b}\), which is given by [19] (Section 13.2), [23]  
Equation 3.0  
where \(\mathbf{Q}_{K}\) is an \(N\times(K+1)\) matrix whose columns form an orthonormal basis for \(\mathcal{K}_{K}(\mathbf{A},\mathbf{b})=\text{span}\left\{\mathbf{b},\mathbf{A }\mathbf{b},\ldots,\mathbf{A}^{K}\mathbf{b}\right\}\), a Krylov subspace. In (9), \(\mathbf{T}_{K}=\mathbf{Q}_{K}^{\top}\mathbf{A}\mathbf{Q}_{K}\) is a \((K+1)\times(K+1)\) tridiagonal Jacobi matrix. The first column of \(\mathbf{Q}_{K}\) is equal to \(\frac{\mathbf{b}}{||\mathbf{b}||}\). The approximation (9) can also be written as \(q_{K}(\mathbf{A})\mathbf{b}\), where \(q_{K}\) is the degree \(K\) polynomial that interpolates the function \(f\) at the \(K+1\) eigenvalues of \(\mathbf{T}_{K}\)[19] (Theorem 13.5), [55]. Thus, unlike classical polynomial approximation methods, such as the truncated Chebyshev expansion, the Lanczos method is indirectly adapted to the spectrum of \(\mathbf{A}\). The Lanczos method differs from proposed method in that \(\mathbf{T}_{K}\) and the Lanczos approximating polynomial \(q_{K}\) depend on the initial vector \(\mathbf{b}\). Specifically, the polynomials \(\{\pi_{k}\}\) generated from the three-term recurrence of the form (8)  
Equation 3.1  
with the \(\{\alpha_{k}\}_{k=0,1,\ldots,K}\) and \(\{\gamma_{k}\}_{k=1,2,\ldots,K}\) coefficients taken from the diagonal and superdiagonal entries of \(\mathbf{T}_{K}\), respectively, are orthogonal with respect to the piecewise-constant measure  
\[\mu(x)= \begin{cases}0,&x<\lambda_{1}\\ \sum_{j=1}^{i}[\mathbf{\hat{b}}(j)]^{2},&\lambda_{i}\leq x<\lambda_{i+1}\\ \sum_{j=1}^{N}[\mathbf{\hat{b}}(j)]^{2}=1,&\lambda_{N}\leq x\end{cases} ,\]  
where \(\mathbf{\hat{b}}=\mathbf{V}^{\top}\mathbf{q}_{1}=\mathbf{V}^{\top}\left( \frac{\mathbf{b}}{||\mathbf{b}||}\right)\), and \(\mathbf{approximation of \(f(\lambda)=e^{-\lambda}\). The second column examines the convergence of relative errors in approximating \(e^{-\mathbf{A}}\mathbf{b}\) for matrices with various spectral distributions, for each of the four methods.