A new variational variance reduction (VVR) technique is developed for improving the efficiency of Monte Carlo multigroup nuclear reactor eigenvalue and eigenfunction calculations. The VVR method employs a variational functional, which requires detailed estimates of both the forward and adjoint fluxes. The direct functional, employed in standard Monte Carlo calculations, requires only limited information concerning the forward flux. The variational functional requires global information about the forward and adjoint fluxes and hence is more expensive to evaluate but is more accurate than the direct functional. In calculations, this increased accuracy outweighs the extra expense, resulting in a more efficient Monte Carlo simulation. In our work, we evaluate the variational functional using Monte Carlo-calculated forward flux estimates and deterministically calculated adjoint flux estimates. Also, we represent the adjoint flux as a low-order polynomial in space and angle, which is accurate for diffusive systems. (In such systems, which are common in reactor analysis problems, the angular flux is locally nearly linear in space and angle.) Using this adjoint representation, we develop specific VVR methods for eigenvalue problems, in which an estimate of the eigenvalue k in a criticality calculation is desired, and eigenfunction problems, in which an estimate of a detector response due to a fission neutron source during a criticality calculation is desired. The resulting VVR method is very efficient for the problems of interest. With a set of example problems, we demonstrate the increased efficiency of the VVR method over standard Monte Carlo.