Particles imaged in cryo-electron microscopy (cryo-EM) often do not exist in a single molecular structure but exhibit significant structural variability. To determine the biological function of a particle, it is important to characterize this variability given noisy projection images, a task known as the heterogeneity problem in single-particle cryo-EM. We present an efficient and accurate method for estimating the low-rank covariance matrix of the three-dimensional voxel structure, which in turn allows for reconstruction of the molecular structure corresponding to each image. The method poses the covariance estimation task as a linear inverse problem in the covariance matrix elements and calculates its least-squares estimator. An efficient algorithm is obtained by exploiting the convolutional structure of the normal equations and reformulating them as a deconvolution problem which can be solved using the conjugate gradient method. The resulting method is the first computationally tractable algorithm for three-dimensional covariance estimation in cryo-EM that is also consistent. Its performance is evaluated on various datasets, demonstrating its effectiveness in solving the heterogeneity problem for both simulations and experimental data.