Authors: Anya Fries, Markus Reichstein, David Blei, Jonas Peters
Organizations: Seminar for Statistics, ETH Zürich, Zürich, Switzerland · Max Planck Institute for Biogeochemistry, Jena, Germany · Department of Computer Science and Department of Statistics, Columbia University, New York, USA
Real-world data in health, economics, and environmental sciences are often collected across heterogeneous domains (such as hospitals, regions, or time periods). In such settings, distributional shifts can make standard PCA unreliable, in that, for example, the leading principal components may explain substantially less variance in unseen domains than in the training domains. Existing approaches (such as FairPCA) have proposed to consider worst-case (rather than average) performance across multiple domains. This work develops a unified framework, called wcPCA, applies it to other objectives (resulting in the novel estimators such as norm-minPCA and norm-maxregret, which are better suited for applications with heterogeneous total variance) and analyzes their relationship. We prove that for all objectives, the estimators are worst-case optimal not only over the observed source domains but also over all target domains whose covariance lies in the convex hull of the (possibly normalized) source covariances. We establish consistency and asymptotic worst-case guarantees of empirical estimators. We extend our methodology to matrix completion, another problem that makes use of low-rank approximations, and prove approximate worst-case optimality for inductive matrix completion. Simulations and two real-world applications on ecosystem-atmosphere fluxes demonstrate marked improvements in worst-case performance, with only minor losses in average performance.
Figures & tables
Figure 1: Proportion of explained variance by poolPCA and norm-maxRegret in the FLUXNET example on source (left) and target (right) regions (for one specific split). Unlike poolPCA , which optimizes a pooled explained variance, norm-maxRegret optimizes a worst-case criterion over the source domains, which, here, yields a worst-case improvement of 7.8%. Section 3.1 proves that this choice comes with a worst-case improvement over a larger class of target domains, too. Indeed, in this example, the worst-case explained variance over the target domains improves by 25.8%. (Over 20 different splits, the median relative increase in worst-case improvements equal 25.6% and 16.6% for source and target domains, respectively; the median decrease in average performance for source domains equals 7.5%; see Section 6.1 for more details.)
Figure 2: Visualization of possible solutions to the rank-1 approximation problem in a specific example. The figure shows contour plots of two distributions N(0,Σ1) (green) and N(0,Σ2) (orange), together with their supports (shaded planes). The methods poolPCA , sepPCA , minPCA , and maxRegret describe different ways to aggregate over different source domains in the objective function (see Section 2.2 for details); their solutions are denoted by Vpool , Vsep , VminPCA , and VmaxRegret , respectively (the values in parentheses indicate the pooled and worst-case explained variance). For example, Vpool and Vsep explain 0% of variance in the worst-case domain, as they are orthogonal to the support of that domain. In contrast, VminPCA maximizes the worst case explained variance (resulting in 36%) and VmaxRegret explains at least 24% of variance in each domain. Details are provided in Example 1 .
Objective function
Normalized
Explained variance
Reconstruction error
Regret
(Def. 2 )
(Def. 3 )
(Def. 4 )
No
minPCA
=
maxRCS
=
maxRegret
=
=
=
=
Yes
norm-minPCA
=
norm-maxRCS
=
norm-maxRegret
=
Table 1: The six worst-case PCA formulations considered in this work: three objectives (explained variance, reconstruction error, and regret) combined with two normalization choices (normalized vs. unnormalized). The symbols = and = indicate whether the solutions coincide (see Theorem 6 ). The equality holds if the methods are applied to the full (and usually unknown) distribution or to its empirical counterpart.
Figure 4
Figure 5: Finite-sample performance of maxRCS ; see Section 5.1.3 . Left : Difference in worst-case loss (reconstruction error) over the convex hull of the population covariance matrices between empirical and population maxRCS . Right : Analogous difference between between empirical maxRCS and empirical poolPCA . Negative values indicate that maxRCS attains lower worst-case error than poolPCA . Across sample sizes, maxRCS outperforms poolPCA .
Figure 6: Regret may be preferable under heterogeneous noise; see Section 5.1.4 . Left : As suggested by the population theory, solutions to maxRCS and maxRegret are almost identical in the noiseless case. With heterogeneous noise, maxRegret behaves similarly as if the data were noiseless. Right : When k is strictly smaller than the rank of the covariance matrices, solutions to maxRCS and maxRegret differ in the noiseless case. Nevertheless, under heterogeneous noise, maxRegret can still be preferable over maxRCS even if the evaluation criteria is reconstruction error.
Figure 7: Worst-case inductive matrix completion with partially observed source domains (90% missing entries; regime (ii)). Left: Difference in average and worst-case test reconstruction error between maxMC and poolMC . In line with Theorem 14 (which covers the population case and allows for some ϵ deviation), maxMC generally improves worst-case performance; in terms of average error, it incurs only minor changes. Right: Worst-case difference as the proportion of missing target entries varies, for different heterogeneity levels (α,β) . maxMC achieves lower median worst-case error than poolMC , with larger gains under stronger domain heterogeneity.
Figure 8: Left : Map of FLUXNET sites, colored by TransCom region. Right : Relative difference in worst-case proportion of explained variance between poolPCA and maxRegret , norm-maxRegret , norm-maxRCS , and avgcovPCA , evaluated on held-out target regions over 20 random train–test splits. Positive values indicate improved worst-case performance relative to pooled PCA. The worst-case objectives generally improve robustness, with norm-maxRegret achieving the strongest gains in median worst-case proportion of explained variance.
Figure 9: Comparison of dimensionality reduction methods on ecosystem function data across continents. Left : Proportion of variance explained on the pooled dataset. Centre : Worst-case proportion of variance explained across continents. Right : Loadings of the top three principal components for poolPCA and norm-maxRCS . norm-maxRCS improves worst-case proportion of explained variance with negligible loss in average performance, while largely retaining the physical interpretability of the principal axes.
Appendix figures & tables7 assets
Supplementary material from the paper’s appendix.
Appendix
Figure B.i: Comparison of projected gradient descent (PGD; used in our experiments) with the optimization schemes of Wang et al. (2025) (StablePCA) and Tantipongpipat et al. (2019) (MW and SDP) for solving minPCA . Results are shown for (p,E)∈{(10,5),(10,50),(50,5)} where p is the number of features and E is the number of domains. Top row : median difference in minimum explained variance relative to that of PGD ( negative values indicate that PGD explains more variance in the worst case) as a function of the number of components k . Bottom row : runtime in seconds for each algorithm. Shaded regions represent the 25%-75% quantiles. PGD attains similar objective values to both baselines, with acceptable runtime.
Figure B.ii: Comparison of PGD, StablePCA, MW, and SDP for solving maxRegret for (p,E)=(10,5) . PGD achieves comparable (or better) maximum regret values (positive values indicate that PGD performs better, i.e., has a lower worst-case regret).
Figure B.iii: Worst-case inductive matrix completion with fully observed source domains; regime (i). Left: Difference in average and worst-case test reconstruction error between maxMC and poolMC . In line with Theorem 14 (which covers the population case and allows for some ϵ deviation), maxMC generally improves worst-case performance; in terms of average error, it incurs only minor changes. Right: Worst-case difference as the proportion of missing target entries varies, for different heterogeneity levels (α,β) . maxMC achieves lower median worst-case error than poolMC , with larger gains under stronger domain heterogeneity.
Figure C.i: Scree plots for each TransCom region (shown in different colours). In Section 6.1 , we use rank-2 approximations.
Figure C.ii: Relative difference in worst-case proportion of explained variance between poolPCA and, from top to bottom, maxRegret , norm-maxRegret , norm-maxRCS , and avgcovPCA . Positive values indicate improved worst-case performance relative to poolPCA . Results are shown for 20 random region splits, evaluated on the source and target regions across varying numbers of principal components. Apart from the rank-one case, the worst-case objectives exhibit improved robustness or similar behavior to poolPCA .
Variable
Description
GPP sat
Maximum gross CO 2 uptake at light saturation
NEP max
Maximum net ecosystem productivity
ET max
Maximum evapotranspiration
EF
Evaporative fraction
EF ampl
EF amplitude
GS max
Maximum dry canopy surface conductance
Appendix
Table 2: Ecosystem functional properties used in the reanalysis of Migliavacca et al. (2021) .
Figure C.iii: Comparison of dimensionality reduction methods on the ecosystem function dataset. Left : Proportion of variance explained on the pooled dataset. Right : Worst-case proportion of variance explained across continents. Methods that account for domain heterogeneity outperform poolPCA in the worst-case continent, with norm-maxRCS most closely matching pooled average performance.
Low-rank approximation represents matrices, tensors, and model computations with fewer components to reduce storage and computation while preserving useful structure.
We address two open questions in streaming PCA via Oja's algorithm: sharp operator-norm convergence for general rank under sub-Gaussian data, and distributional inference for the resulting subspace estimator. Existing convergence analyses, even in the rank-one case, either assume bounded data or leave non-vanishing remainder terms that prevent adaptation to a polynomially vanishing tail spectrum, while existing distributional results are confined to the rank-one case. Our convergence theory removes these remainder terms and yields a sharp rate. In the dense-tail spiked covariance regime, this rate matches the minimax rate up to logarithmic factors. More generally, we prove a matching lower bound, up to logarithmic factors, across both dense-tail and sparse-tail regimes under a mild nondegeneracy condition. The analysis yields a linearization of Oja's iterates, which in turn enables a high-dimensional Gaussian approximation for the general-rank subspace estimation error with an explicit limiting covariance. We also establish a row-wise Gaussian approximation over convex sets for the aligned difference, recovering prior rank-one results as special cases. For practical inference, we develop an online multiplier bootstrap algorithm and prove its consistency. Beyond streaming PCA, our techniques contribute to Gaussian approximation and bootstrap inference for nonconvex stochastic approximation.
Principal Component Analysis (PCA) is a fundamental tool for statistical dimension reduction, yet no canonical extension to distribution-valued data has emerged within the framework of optimal transport. While Wasserstein barycenters and Wasserstein least-squares regression provide canonical liftings of Euclidean averaging and linear regression in the Wasserstein space, an analogous principle for PCA has remained unexplored. In this paper, we introduce Lifted PCA, a canonical extension of Euclidean PCA to Wasserstein space, obtained through the framework of lifting functionals in the space of probability measures by convex extension. Starting from the classical variational formulations of PCA on the Stiefel and Grassmann manifolds, we derive two equivalent minimization problems over probability measures supported on these manifolds that correspond to Wasserstein barycenters with generalized costs. We establish an exact analogue of the equivalence between the Stiefel and Grassmann formulations of Euclidean PCA, preserving PCA's original variational structure under lifting. We further derive dual formulations and optimality conditions, and establish equivalences with multimarginal optimal transport and matching-for-teams problems. We also establish statistical consistency for empirical estimators and develop computational algorithms based on the Wasserstein barycenter formulation over the Stiefel manifold. Whereas Euclidean PCA summarizes vectors through projections onto an optimal low-dimensional linear subspace, Lifted PCA summarizes probability distributions through an optimal probability measure supported on multiple low-dimensional subspaces. Numerical experiments on distributional datasets illustrate the practical performance of the proposed methodology and its benefits for distributional data analysis.
Jérémie Bigot, Clément Bonet
Institut de Mathématiques de Bordeaux, Université de Bordeaux, CNRS (UMR 5251) · CMAP, Ecole Polytechnique, Institut Polytechnique de Paris, CNRS (UMR 7641)
Distributed principal component analysis (PCA) produces node-level estimates of both a mean vector and a principal subspace. Robustly aggregating these heterogeneous objects requires a relative scale between mean error and subspace error. We study a scale-calibrated median-of-means estimator for this problem using the product geometry of Euclidean space and the Grassmann manifold. A node-level PCA expansion shows that the mean component has the usual linear influence, whereas the subspace component is an eigengap-weighted covariance perturbation. We prove a local reduction showing that the proposed product-manifold median-of-means estimator is asymptotically equivalent to a scaled spatial median of node influence errors. This yields fixed-node non-Gaussian limits, growing-node Gaussian limits with finite-block bias, and an explicit scale-dependent covariance formula. We propose robust block-scale and inference-optimal calibration rules, establish high-probability median-of-means bounds, characterize factorwise bad-node influence, and prove node-bootstrap validity. Simulations and large-scale single-cell RNA-seq data show that scale calibration adapts to eigengap-driven subspace uncertainty and provides a robust distributed PCA summary.
Kisung You
Department of Mathematics, Baruch College · Department of Mathematics, The Graduate Center, City University of New York