AbstractThe Hutchinson method can be used to estimate the trace of a function of a matrix A, that is, , stochastically. With the variance of the Hutchinson estimator decreasing as , techniques to reduce the variance of this estimator become necessary due to the high accuracy required by some calculations in scientific computing. In lattice QCD, where A is the Dirac operator D stemming from a discretization of the Dirac equation, the accuracy of some calculations involving renders them practically unfeasible if one uses the Hutchinson estimator to approximate the trace. One such variance reduction technique, known as multigrid multilevel Monte Carlo, has been introduced recently to deal with for the Laplace 2D, Gauge Laplacian, and Schwinger problems. We use it here in the context of lattice QCD, and moreover, we analyze its impact when dealing with , where is the matrix that emerges when a displacement in one of the four spacetime dimensions of the lattice is introduced. This latter trace is of importance when computing generalized parton distributions in lattice QCD.