  Hutchinson's method estimates the trace of a matrix function $f(D)$
stochastically using samples $\tau^Hf(D)\tau$, where the components of the
random vectors $\tau$ obey an isotropic probability distribution. Estimating
the trace of the inverse of a discretized Dirac operator or variants thereof
have become a major challenge in lattice QCD simulations, as they represent the
disconnected contribution to certain observables. The Hutchinson Monte Carlo
sampling, however, suffers from the fact that its accuracy depends
quadratically on the sample size, making higher precision estimation very
expensive. Meyer, Musco, Musco and Woodruff recently proposed an enhancement of
Hutchinson's method, termed \texttt{Hutch++}, in which the sample space is
enriched by several vectors of the form $f(D)\zeta$, $\zeta$ a random vector as
in Hutchinson's method. Theoretical analyses show that under certain
circumstances the number of these added sample vectors can be chosen in a way
to reduce the dependence of the variance of the resulting estimator from the
number $N$ of samples from $\mathcal{O}(1/N)$ to $\mathcal{O}(1/N^2)$.
  In this study we combine \texttt{Hutch++} with our recently suggested
multigrid multilevel Monte Carlo approach. We present results for the Schwinger
discretization of the $2$-di\-men\-si\-onal Dirac operator, revealing that the
two approaches contribute additively to variance reduction.
