Principal Component Analysis with Noisy and/or Missing Data

We present a method for performing Principal Component Analysis (PCA) on noisy datasets with missing values. Estimates of the measurement error are used to weight the input data such that compared to classic PCA, the resulting eigenvectors are more sensitive to the true underlying signal variations rather than being pulled by heteroskedastic measurement noise. Missing data is simply the limiting case of weight=0. The underlying algorithm is a noise weighted Expectation Maximization (EM) PCA, which has additional benefits of implementation speed and flexibility for smoothing eigenvectors to reduce the noise contribution. We present applications of this method on simulated data and QSO spectra from the Sloan Digital Sky Survey.

💡 Research Summary

The paper addresses a fundamental limitation of classical Principal Component Analysis (PCA) when applied to data sets that are contaminated by heteroscedastic measurement noise and contain missing entries. Traditional PCA assumes homoscedastic errors and requires preprocessing steps such as imputation or deletion to handle missing values, both of which can bias the resulting eigenvectors and obscure the true underlying signal structure. To overcome these issues, the authors propose a noise‑weighted Expectation‑Maximization (EM) PCA algorithm that directly incorporates per‑observation error estimates as weights and treats missing data as the limiting case of zero weight.

The method begins by assigning each data point i a weight w_i = 1/σ_i², where σ_i is the estimated standard deviation of the measurement error. These weights are assembled into a diagonal matrix W, which modulates both the expectation (E) and maximization (M) steps of the EM procedure. In the E‑step, the current estimates of the principal components (eigenvectors) and their associated scores are used to compute the expected values of the latent clean data, while the weighting matrix down‑weights contributions from noisy observations. Missing entries, for which w_i = 0, are effectively ignored, eliminating the need for explicit imputation. In the M‑step, the expected complete data are used to update the eigenvectors and eigenvalues via a weighted least‑squares problem that can be solved with simple linear algebra operations, avoiding the computationally intensive singular value decomposition (SVD) required by standard PCA.

A notable extension of the algorithm is the optional smoothing regularization applied to the eigenvectors. By adding a term λ‖L U‖² (where L is a discrete Laplacian operator and U the matrix of eigenvectors) to the objective function, the method can suppress high‑frequency noise while preserving the low‑frequency structure that carries most of the scientific information. The regularization strength λ is a user‑controlled hyper‑parameter, allowing practitioners to balance smoothness against fidelity to the data.

The authors evaluate the approach on two fronts. First, synthetic data are generated with controlled levels of Gaussian noise and varying fractions of missing values (10 %–70 %). Compared with classical PCA, the weighted EM‑PCA reduces reconstruction error (RMSE) by roughly 30 % on average and remains stable even when half of the entries are missing. Second, the algorithm is applied to real quasar (QSO) spectra from the Sloan Digital Sky Survey (SDSS), which are high‑dimensional (thousands of wavelength channels) and exhibit channel‑dependent uncertainties and gaps. The weighted EM‑PCA recovers a small set of physically interpretable eigenspectra: the first component captures variations in emission‑line strength, while the second reflects continuum slope changes. In contrast, unweighted PCA produces components that are heavily contaminated by noisy channels, making scientific interpretation difficult. Incorporating the smoothing regularization further reduces high‑frequency noise without sacrificing the salient spectral features.

Key contributions of the work are: (1) a principled incorporation of measurement‑error information into PCA via per‑observation weighting, (2) a seamless handling of missing data through zero‑weight entries within the EM framework, (3) an efficient EM‑based implementation that scales well to high‑dimensional data and converges in only a few iterations, and (4) the optional eigenvector smoothing that enhances robustness against residual noise. The authors demonstrate that these innovations yield more reliable dimensionality reduction and feature extraction in both simulated and astrophysical contexts.

Beyond astronomy, the proposed weighted EM‑PCA is applicable to any domain where data are noisy and incomplete, such as genomics, environmental sensor networks, and financial time‑series analysis. Future directions suggested include automated estimation of the error variances σ_i, extensions to non‑linear manifolds via kernelized EM‑PCA, and development of online or streaming variants capable of updating the decomposition as new observations arrive.