  In lattice QCD the calculation of disconnected quark loops from the trace of
the inverse quark matrix has large noise variance. A multilevel Monte Carlo
method is proposed for this problem that uses different degree polynomials on a
multilevel system. The polynomials are developed from the GMRES algorithm for
solving linear equations. To reduce orthogonalization expense, the highest
degree polynomial is a composite or double polynomial found with a polynomial
preconditioned GMRES iteration. Matrix deflation is used in three different
ways: in the Monte Carlo levels, in the main solves, and in the deflation of
the highest level double polynomial. A numerical comparison with optimized
Hutchinson is performed on a quenched \(24^4\) lattice. The results demonstrate
that the new Multipolynomial Monte Carlo method can significantly improve the
trace computation for matrices that have a difficult spectrum due to small
eigenvalues.}
