## 1 Introduction  
Efficiently computing \(f(\mathbf{A})\mathbf{b}\), a function of a large, sparse Hermitian matrix times a vector, is an important component in numerous signal processing, machine learning, applied mathematics, and computer science tasks. Application examples include graph-based semi-supervised learning methods [1; 2; 3]; graph spectral filtering in graph signal processing [4]; convolutional neural networks/deep learning [5; 6]; clustering [7; 8]; approximating the spectral density of a large matrix [9]; estimating the numerical rank of a matrix [10; 11]; approximating spectral sums such as the log-determinant of a matrix [12] or the trace of a matrix inverse for applications in physics, biology, information theory, and other disciplines [13]; solving semidefinite programs [14]; simulating random walks [15] (Chapter 8); and solving ordinary and partial differential equations [16; 17; 18]. References [19] (Chapter 13), [20; 21; 22] survey different approaches to this well-studied problem of efficiently computing  
\[f(\mathbf{A})\mathbf{b}:=\mathbf{V}f(\mathbf{A})\mathbf{V}^{\top}\mathbf{b}, \tag{1}\]  
where the columns of \(\mathbf{V}\) are the eigenvectors of the Hermitian matrix \(\mathbf{A}\in\mathbb{R}^{N\times N}\); \(\mathbf{A}\) is a diagonal matrix whose diagonal elements are the corresponding eigenvalues of \(\mathbf{A}\), which we denote by \(\lambda_{1},\lambda_{2},\ldots,\lambda_{N}\); and \(f(\mathbf{\Lambda})\) is a diagonal matrix whose \(k\)th diagonal entry is given by \(f(\lambda_{k})\). For large matrices, it is not practical to explicitly compute the eigenvalues of \(\mathbf{A}\) in order to approximate (1). Rather, the most common techniques, all of which avoid a full eigendecomposition of \(\mathbf{A}\), include (i) truncated orthogonal polynomial expansions, including Chebyshev [23; 24; 25] and Jacobi; (ii) rational approximations [21] (Section 3.4); (iii) Krylov subspace methods such as the Lanczos method [23; 26; 27; 28; 29]; and (iv) quadrature/contour integral methods [19] (Section 13.3). Our focus in this work is on polynomial approximation methods. Let \(p_{K}(\lambda)=c_{0}+\sum_{k=1}^{K}c_{k}\lambda^{k}\) be a degree \(K\) polynomial approximation to the function \(f\) on a known interval \([\underline{\lambda},\overline{\lambda}]\) containing all of the eigenvalues of \(\mathbf{A}\). Then the approximation \(p_{K}(\mathbf{A})\mathbf{b}\) can be computed recursively, either through a three-term recurrence for specific types of polynomials (see Section 3 for more details), or through a nested multiplication iteration [30] (Section 9.2.4), letting \(\mathbf{x}^{(0)}=c_{K}\mathbf{b}\), and then iterating  
\[\mathbf{x}^{(l)}=c_{K-l}\mathbf{b}+\mathbf{A}\mathbf{x}^{(l-1)},\;l=1,2, \ldots,K. \tag{2}\]  
The computational cost of either of these approaches is dominated by multiplying the sparse matrix \(\mathbf{A}\) by \(K\) different vectors. The approximation error is bounded by  
\[||f(\mathbf{A})-p_{K}(\mathbf{A})||_{2} =\max_{\ell=1,2,\ldots,N}|f(\lambda_{\ell})-p_{K}(\lambda_{\ell})| \tag{3}\] \[\leq\sup_{\lambda\in[\underline{\lambda}]}|f(\lambda)-p_{K}( \lambda)|. \tag{4}\]  
If, for example, \(p_{K}\) is a degree \(K\) truncated Chebyshev series approximation of an analytic function \(f\), the upper bound in (4) converges geometrically to \(0\) as \(K\) increases, at a rate of \(\mathcal{O}\left(\rho^{-K}\right)\), where \(\rho\) is the radius of an open Bernstein ellipse on which \(f\) is analytic and bounded (see, e.g., [31] (Theorem 5.16), [32] (Theorem 8.2)). In addition to the computational efficiency and convergence guarantees, a third advantage of polynomial approximation methods is that they can be implemented in a distributed setting [33].