  We present a new sublinear time algorithm for approximating the spectral
density (eigenvalue distribution) of an $n\times n$ normalized graph adjacency
or Laplacian matrix. The algorithm recovers the spectrum up to $\epsilon$
accuracy in the Wasserstein-1 distance in $O(n\cdot \text{poly}(1/\epsilon))$
time given sample access to the graph. This result compliments recent work by
David Cohen-Steiner, Weihao Kong, Christian Sohler, and Gregory Valiant (2018),
which obtains a solution with runtime independent of $n$, but exponential in
$1/\epsilon$. We conjecture that the trade-off between dimension dependence and
accuracy is inherent.
  Our method is simple and works well experimentally. It is based on a
Chebyshev polynomial moment matching method that employees randomized
estimators for the matrix trace. We prove that, for any Hermitian $A$, this
moment matching method returns an $\epsilon$ approximation to the spectral
density using just $O({1}/{\epsilon})$ matrix-vector products with $A$. By
leveraging stability properties of the Chebyshev polynomial three-term
recurrence, we then prove that the method is amenable to the use of coarse
approximate matrix-vector products. Our sublinear time algorithm follows from
combining this result with a novel sampling algorithm for approximating
matrix-vector products with a normalized graph adjacency matrix.
  Of independent interest, we show a similar result for the widely used
\emph{kernel polynomial method} (KPM), proving that this practical algorithm
nearly matches the theoretical guarantees of our moment matching method. Our
analysis uses tools from Jackson's seminal work on approximation with positive
polynomial kernels.
