We present a control-variate method for reducing the variance of power spectrum covariance matrix estimates from simulations of large-scale structure. The key idea is to pair each mock simulation with a cheap Zeldovich-approximation realization sharing the same initial conditions, and to use the known statistical properties of the Zeldovich field to remove correlated sample variance from the covariance estimator. Under a Gaussian disconnected approximation, we derive fully analytic expressions for both the optimal control-variate coefficient, $\beta (k,\ell ;k\prime ,\ell \prime )$, and the corresponding correlation, $\rho (k,\ell ;k\prime ,\ell \prime )$, in terms of the auto- and cross-power spectra of the target and control fields. In the monopole case, the correlation takes the particularly simple form $\rho (k,k\prime )=r^2(k),r^2(k\prime )$, where $r(k)$ is the standard cross-correlation coefficient between the target and Zeldovich fields, implying that covariance estimation remains highly efficient whenever the two fields are strongly correlated. For masked redshift-space lognormal mocks, resembling Luminous Red Galaxies from the Dark Energy Spectroscopic Instrument (DESI), we find that the control-variate estimator reduces the variance of the covariance matrix by approximately an order of magnitude on large scales, $k \lesssim 0.05\,h\,{\rm Mpc}^{-1}$, precisely where accurate covariance estimation is most challenging. The gains are smaller for higher $k$ but typically accelerate convergence by a factor of 2-3, substantially lowering the computational cost of covariance estimation for current and upcoming large-scale structure surveys. We demonstrate that after a computationally cheap parametric-bootstrap calibration of the precision-matrix bias, the control-variate covariance yields a ${\chi }^2$ distribution consistent with the expected number of degrees of freedom. Due to its simplicity, this method is readily implementable in current imaging and spectroscopic surveys (e.g., DESI, , LSST, PFS, ).