  The efficient simulation of the mean value of a non-linear functional of the
solution to a linear stochastic partial differential equation (SPDE) with
additive Gaussian noise is considered. A Galerkin finite element method is
employed along with an implicit Euler scheme to arrive at a fully discrete
approximation of the mild solution to the equation. A scheme is presented to
compute the covariance of this approximation, which allows for rapid sampling
in a Monte Carlo method. This is then extended to a multilevel Monte Carlo
method, for which a scheme to compute the cross-covariance between the
approximations at different levels is presented. In contrast to traditional
path-based methods it is not assumed that the Galerkin subspaces at these
levels are nested. The computational complexities of the presented schemes are
compared to traditional methods and simulations confirm that, under suitable
assumptions, the costs of the new schemes are significantly lower.
