  We construct quasi-Monte Carlo methods to approximate the expected values of
linear functionals of Galerkin discretizations of parametric operator equations
which depend on a possibly infinite sequence of parameters. Such problems arise
in the numerical solution of differential and integral equations with random
field inputs. We analyze the regularity of the solutions with respect to the
parameters in terms of the rate of decay of the fluctuations of the input
field. If $p\in (0,1]$ denotes the "summability exponent" corresponding to the
fluctuations in affine-parametric families of operators, then we prove that
deterministic "interlaced polynomial lattice rules" of order $\alpha = \lfloor
1/p \rfloor+1$ in $s$ dimensions with $N$ points can be constructed using a
fast component-by-component algorithm, in $\mathcal{O}(\alpha\,s\, N\log N +
\alpha^2\,s^2 N)$ operations, to achieve a convergence rate of
$\mathcal{O}(N^{-1/p})$, with the implied constant independent of $s$. This
dimension-independent convergence rate is superior to the rate
$\mathcal{O}(N^{-1/p+1/2})$, for $2/3\leq p\leq 1$ recently established for
randomly shifted lattice rules under comparable assumptions. In our analysis we
use a non-standard Banach space setting and introduce "smoothness-driven
product and order dependent (SPOD)" weights for which we show fast CBC
construction.
