  A first-order, Monte Carlo ensemble method has been recently introduced for
solving parabolic equations with random coefficients in [26], which is a
natural synthesis of the ensemble-based, Monte Carlo sampling algorithm and the
ensemble-based, first-order time stepping scheme. With the introduction of an
ensemble average of the diffusion function, this algorithm leads to a single
discrete system with multiple right-hand sides for a group of realizations,
which could be solved more efficiently than a sequence of linear systems. In
this paper, we pursue in the same direction and develop a new multilevel Monte
Carlo ensemble method for solving random parabolic partial differential
equations. Comparing with the approach in [26], this method possesses a
high-order accuracy in time and further reduces the computational cost by using
the multilevel Monte Carlo method. Rigorous numerical analysis shows the method
achieves the optimal rate of convergence. Several numerical experiments are
presented to illustrate the theoretical results.
