  We test a method to reduce unwanted sample variance when predicting
Lyman-$\alpha$ (ly$\alpha$) forest power spectra from cosmological
hydrodynamical simulations. Sample variance arises due to sparse sampling of
modes on large scales and propagates to small scales through non-linear
gravitational evolution. To tackle this, we generate initial conditions in
which the density perturbation amplitudes are {\it fixed} to the ensemble
average power spectrum -- and are generated in {\it pairs} with exactly
opposite phases. We run $50$ such simulations ($25$ pairs) and compare their
performance against $50$ standard simulations by measuring the ly$\alpha$ 1D
and 3D power spectra at redshifts $z=2$, 3, and 4. Both ensembles use periodic
boxes of $40$ Mpc/h containing $512^3$ particles each of dark matter and gas.
As a typical example of improvement, for wavenumbers $k=0.25$ h/Mpc at $z=3$,
we find estimates of the 1D and 3D power spectra converge $34$ and $12$ times
faster in a paired-fixed ensemble compared with a standard ensemble. We
conclude that, by reducing the computational time required to achieve fixed
accuracy on predicted power spectra, the method frees up resources for
exploration of varying thermal and cosmological parameters -- ultimately
allowing the improved precision and accuracy of statistical inference.
