  We introduce the multivariate decomposition finite element method for
elliptic PDEs with lognormal diffusion coefficient $a=\exp(Z)$ where $Z$ is a
Gaussian random field defined by an infinite series expansion
$Z(\boldsymbol{y}) = \sum_{j\ge1} y_j\,\phi_j$ with $y_j\sim\mathcal{N}(0,1)$
and a given sequence of functions $\{\phi_j\}_{j\ge1}$. We use the MDFEM to
approximate the expected value of a linear functional of the solution of the
PDE which is an infinite-dimensional integral over the parameter space. The
proposed algorithm uses the multivariate decomposition method (MDM) to compute
the infinite-dimensional integral by a decomposition into finite-dimensional
integrals, which we resolve using quasi-Monte Carlo (QMC) methods, and for
which we use the finite element method (FEM) to solve different instances of
the PDE.
  We develop higher-order quasi-Monte Carlo rules for integration over the
finite-dimensional Euclidean space with respect to the Gaussian distribution by
use of a truncation strategy. By linear transformations of interlaced
polynomial lattice rules from the unit cube to a multivariate box of the
Euclidean space we achieve higher-order convergence rates for functions
belonging to a class of anchored Gaussian Sobolev spaces, taking into account
the truncation error.
  Under appropriate conditions, the MDFEM achieves higher-order convergence
rates in term of error versus cost, i.e., to achieve an accuracy of
$O(\epsilon)$ the computational cost is $O(\epsilon^{-1/\lambda-d'/\lambda}) =
O(\epsilon^{-(p^*+d'/\tau)/(1-p^*)})$ where $\epsilon^{-1/\lambda}$ and
$\epsilon^{-d'/\lambda}$ are respectively the cost of the quasi-Monte Carlo
cubature and the finite element approximations, with $d' = d \, (1+\delta')$
for some $\delta' \ge 0$ and $d$ the physical dimension, and $0 < p^* \le
(2+d'/\tau)^{-1}$ is a parameter representing the sparsity of
$\{\phi_j\}_{j\ge1}$.
