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.