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.