The distribution of centered Gaussian random fields (GRFs) indexed by compacta such as smooth, bounded Euclidean domains or smooth, compact and orientable manifolds is determined by their covariance operators. We consider centered GRFs given as variational solutions to coloring operator equations driven by spatial white noise, with an elliptic self-adjoint pseudodifferential coloring operator from the Hörmander class. This includes the Matérn class of GRFs as a special case. Using biorthogonal multiresolution analyses on the manifold, we prove that the precision and covariance operators, respectively, may be identified with bi-infinite matrices and finite sections may be diagonally preconditioned rendering the condition number independent of the dimension p of this section. We prove that a tapering strategy by thresholding applied on finite sections of the bi-infinite precision and covariance matrices results in optimally numerically sparse approximations. That is, asymptotically only linearly many nonzero matrix entries are sufficient to approximate the original section of the bi-infinite covariance or precision matrix using this tapering strategy to arbitrary precision. The locations of these nonzero matrix entries can be determined a priori. The tapered covariance or precision matrices may also be optimally diagonally preconditioned. Analysis of the relative size of the entries of the tapered covariance matrices motivates novel, multilevel Monte Carlo (MLMC) oracles for covariance estimation, in sample complexity that scales log-linearly with respect to the number p of parameters. In addition, we propose and analyze novel compressive algorithms for simulating and kriging of GRFs. The complexity (work and memory vs. accuracy) of these three algorithms scales near-optimally in terms of the number of parameters p of the sample-wise approximation of the GRF in Sobolev scales.

We formulate standard and multilevel Monte Carlo methods for the kth moment M_ε^k[ξ] of a Banach space valued random variable ξ:Ω→E, interpreted as an element of the k-fold injective tensor product space ⊗_ε^kE. For the standard Monte Carlo estimator of M_ε^k[ξ], we prove the k-independent convergence rate [Formula presented] in the L_q(Ω;⊗_ε^kE)-norm, provided that (i) ξ∈L_kq(Ω;E) and (ii) q∈[p,∞), where p∈[1,2] is the Rademacher type of E. By using the fact that Rademacher averages are dominated by Gaussian sums combined with a version of Slepian's inequality for Gaussian processes due to Fernique, we moreover derive corresponding results for multilevel Monte Carlo methods, including a rigorous error estimate in the L_q(Ω;⊗_ε^kE)-norm and the optimization of the computational cost for a given accuracy. Whenever the type of the Banach space E is p=2, our findings coincide with known results for Hilbert space valued random variables. We illustrate the abstract results by three model problems: second-order elliptic PDEs with random forcing or random coefficient, and stochastic evolution equations. In these cases, the solution processes naturally take values in non-Hilbertian Banach spaces. Further applications, where physical modeling constraints impose a setting in Banach spaces of type p<2, are indicated.