  This article presents a randomized matrix-free method for approximating the
trace of $f({\bf A})$, where ${\bf A}$ is a large symmetric matrix and $f$ is a
function analytic in a closed interval containing the eigenvalues of ${\bf A}$.
Our method uses a combination of stochastic trace estimation (i.e.,
Hutchinson's method), Chebyshev approximation, and multilevel Monte Carlo
techniques. We establish general bounds on the approximation error of this
method by extending an existing error bound for Hutchinson's method to
multilevel trace estimators. Numerical experiments are conducted for common
applications such as estimating the log-determinant, nuclear norm, and Estrada
index, and triangle counting in graphs. We find that using multilevel
techniques can substantially reduce the variance of existing single-level
estimators.
