  We present and test a method that dramatically reduces variance arising from
the sparse sampling of wavemodes in cosmological simulations. The method uses
two simulations which are fixed (the initial Fourier mode amplitudes are fixed
to the ensemble average power spectrum) and paired (with initial modes exactly
out of phase). We measure the power spectrum, monopole and quadrupole
redshift-space correlation functions, halo mass function and reduced bispectrum
at $z=1$. By these measures, predictions from a fixed pair can be as precise on
non-linear scales as an average over 50 traditional simulations. The fixing
procedure introduces a non-Gaussian correction to the initial conditions; we
give an analytic argument showing why the simulations are still able to predict
the mean properties of the Gaussian ensemble. We anticipate that the method
will drive down the computational time requirements for accurate large-scale
explorations of galaxy bias and clustering statistics, enabling more precise
comparisons with theoretical models, and facilitating the use of numerical
simulations in cosmological data interpretation.
