  We present a Multi-Index Quasi-Monte Carlo method for the solution of
elliptic partial differential equations with random coefficients. By combining
the multi-index sampling idea with randomly shifted rank-1 lattice rules, the
algorithm constructs an estimator for the expected value of some functional of
the solution. The efficiency of this new method is illustrated on a
three-dimensional subsurface flow problem with lognormal diffusion coefficient
with underlying Mat\'ern covariance function. This example is particularly
challenging because of the small correlation length considered, and thus the
large number of uncertainties that must be included. We show numerical evidence
that it is possible to achieve a cost inversely proportional to the requested
tolerance on the root-mean-square error, for problems with a smoothly varying
random field
