Harmonic Mean PCA

• HM-PCA is a robust dimension-reduction technique that aggregates local covariance inverses via the harmonic mean to mitigate eigenvalue ordering errors under contamination.
• It preserves classical PCA efficiency for clean data while enhancing robustness in distributed, heavy-tailed, or contaminated environments.
• The method computes the inverse of an averaged set of matrix inverses using spectral decomposition and ridge regularization for stable subspace estimation.

Harmonic Mean Principal Component Analysis (HM-PCA) designates a class of dimension-reduction and subspace-estimation techniques in which the harmonic mean is used to aggregate covariance or scatter matrices—usually in the context of distributed or robust principal component analysis. Within the recently formalized $\phi$-PCA framework (Hung et al., 15 Oct 2025), HM-PCA is identified by the transformation $\phi(u)=u^{-1}$ applied to the eigenvalues before averaging. HM-PCA methods are distinguished by their optimal robustness against outlier-driven eigenvalue ordering errors and preservation of classical PCA efficiency in the absence of contamination. The methodology and theory of HM-PCA unify and extend earlier proposals involving harmonic averaging of positive semidefinite matrices, incorporating developments in perturbation inequalities (Sababheh, 2018), spectral asymptotics (Lodhia, 2019), and generalized matrix mean aggregation for distributed PCA (Jou et al., 1 Oct 2024).

1. Mathematical Structure of HM-PCA

Let $X \in \mathbb{R}^{n \times p}$ denote the (partitioned) data matrix, and suppose $X$ is split into $m$ (possibly distributed) subsets. Each subset yields a local covariance matrix $S_k$ (for $k=1,\ldots,m$). The $\phi$-PCA framework aggregates these as

$$\widehat{\Sigma}_\phi^{(m)} = \phi^{-1}\left( \frac{1}{m} \sum_{k=1}^{m} \phi(\widehat{S}_k) \right),$$

where $\phi$ is a monotone function acting on the spectrum of each $S_k$. For HM-PCA, $\phi(u) = u^{-1}$, and thus

$$\widehat{\Sigma}_{\text{HM}}^{(m)} = \left( \frac{1}{m} \sum_{k=1}^{m} \widehat{S}_k^{-1} \right)^{-1}.$$

This operation is defined on the cone of symmetric positive definite matrices and can be regularized with a ridge parameter $\epsilon > 0$ as $\widehat{S}_k + \epsilon I$. The method admits a spectral (eigendecomposition) implementation: $$\widehat{S}_k = \sum_{j=1}^p \lambda_{k,j} \gamma_{k,j} \gamma_{k,j}^\top,\quad \phi(\widehat{S}_k) = \sum_{j=1}^p \phi(\lambda_{k,j}) \gamma_{k,j} \gamma_{k,j}^\top,$$ with aggregation performed eigenvalue-wise and reconstruction via inverse spectral mapping.

In distributed PCA, the matrix harmonic mean underlies $\beta$-DPCA with $\beta=-1$ (Jou et al., 1 Oct 2024), constructed as

$$\overline{M}_{-1} = \left(\frac{1}{m}\sum_{\ell=1}^m M^{(\ell),-1}\right)^{-1},$$

where each $M^{(\ell)}$ is a projection matrix or covariance estimate from a local node.

For classical random-matrix-theoretic formulations, given $n$ independent Wishart matrices $W_i$ (sample covariances), the harmonic mean is

$H = n \left( \sum_{i=1}^n W_i^{-1} \right)^{-1},$

which exhibits distinct limiting spectral behavior compared to the arithmetic mean (Lodhia, 2019).

2. Theoretical Properties: Robustness and Efficiency

Ordering-Robustness:

HM-PCA is specifically designed to provide optimal protection to the eigenvalue (and hence eigenvector) ordering under data contamination. The harmonic mean transformation down-weights large (outlying) eigenvalues relative to small ones, reducing the influence of contaminated partitions on the aggregated covariance. The overarching theoretical result [(Hung et al., 15 Oct 2025), Thm 3] quantifies the maximal gain in robustness for HM-PCA compared to arithmetic or geometric mean aggregations, especially when outliers reside in the noise subspace orthogonal to the true signal.

Efficiency Preservation:

Crucially, HM-PCA preserves the first-order (asymptotic) normality and efficiency properties of standard PCA in the absence of outliers. The distribution of eigenvalues and eigenvectors of $\widehat{\Sigma}_{\phi}^{(m)}$ are asymptotically identical to those of single-batch (arithmetic) PCA when all partitions are clean, regardless of the choice of $\phi$ [(Hung et al., 15 Oct 2025), Thm 1].

Perturbation Analysis:

For the influence function expansions, the leading sensitivity of both eigenvalues and eigenvectors is unchanged among $\phi$-PCA methods. However, second-order (quadratic) terms, which dictate ordering robustness, are strictly smaller for HM-PCA under partition-level contamination [(Hung et al., 15 Oct 2025), Thm 2].

Spectral Inequalities:

HM-PCA constructions leverage matrix harmonic-arithmetic mean inequalities (Sababheh, 2018), ensuring that the harmonic mean of SPD matrices is always less than the arithmetic mean in the Löwner order, tightly bounded by explicit quadratic corrections. The spectral decomposition facilitates precise control of the effect of harmonic aggregation on each eigenmode.

3. Algorithmic Implementation and Computation

HM-PCA proceeds via the following steps:

1. Partition the data into $m$ subsets (for distributed, robust, or computational considerations).
2. Compute the local sample covariance (or robust scatter) matrices $\widehat{S}_k$ for each subset.
3. (Optional) Regularize each $$\widehat{S}_k \leftarrow \widehat{S}_k + \epsilon I$$ for stability.
4. Compute $$\widehat{\Sigma}_{\text{HM}}^{(m)} = (\frac{1}{m} \sum_{k=1}^m \widehat{S}_k^{-1})^{-1}$$.
5. Extract the leading $r$ eigenvectors of $\widehat{\Sigma}_{\text{HM}}^{(m)}$ for subspace estimation.

For generalized mean aggregation (matrix $\beta$-means), the following formula is used (Jou et al., 1 Oct 2024): $$\overline{M}_\beta = \left( \frac{1}{m} \sum_{\ell=1}^m [M^{(\ell)}]^\beta \right)^{1/\beta},$$ with $\beta=-1$ yielding HM-PCA.

In numerical practice, all matrix means are computed spectrally: eigen-decompose each matrix, raise eigenvalues to the appropriate power (inverse for $\beta=-1$), sum, apply the mean, and reconstruct.

For high-dimensional random matrix aggregation, the limiting spectral distribution of the harmonic mean, derived using free probability, admits explicit closed-form Stieltjes transform equations (Lodhia, 2019).

4. Comparative Analysis: Arithmetic, Geometric, and Harmonic Aggregation

Method	Aggregation Function $\phi(u)$	Ordering Robustness	Efficiency (Clean Data)	Suitability under Contamination
AM-PCA	$u$	Baseline	Optimal	Sensitive to outliers, non-robust
GM-PCA	$\log u$ (limit $\beta\to0$)	Moderate	Optimal	Partially robust
HM-PCA	$u^{-1}$	Maximal	Optimal	Optimal, especially for extreme outliers

HM-PCA outperforms AM-PCA and GM-PCA in scenarios where eigenvalue ordering underlies selection of the principal subspace and where outlier partitions or local corruptions are present. It has been shown (Jou et al., 1 Oct 2024, Hung et al., 15 Oct 2025) that for heavy-tailed, contaminated, or distributed data, HM-PCA maintains accurate eigenspace estimation while minimizing the risk of misordering the leading components.

5. Practical Applications and Computational Considerations

Distributed and Federated PCA:

HM-PCA is naturally suited to distributed settings. Each computational node or data silo computes a local covariance; the central server aggregates via the harmonic mean. This approach enhances robustness not only to classical adversarial outliers, but also to systemic differences among heterogeneous data sources (Jou et al., 1 Oct 2024, Hung et al., 15 Oct 2025).

High-dimensional/Contaminated Data:

For modern "big data" settings—such as image analysis or genomics—HM-PCA can recover principal subspaces robust to blocks of contaminated data or heavy tails. Simulation and real-data studies, including partitioned MNIST reconstructions, demonstrate that standard PCA suffers from dramatic eigenspace reordering, whereas HM-PCA preserves feature structure (Hung et al., 15 Oct 2025).

Spectral Regularization and Conditioning:

The harmonic mean diminishes the influence of large-eigenvalue blocks, mitigating overfitting to high-variance local distortions. However, care must be taken to avoid amplifying conditioning issues from near-singular local covariances; ridge regularization is commonly employed.

Computational Complexity:

HM-PCA requires inversion and spectral composition of $m$ local covariance matrices. While more expensive than direct averaging, these steps scale well in distributed architectures and can leverage efficient inversion for structured/sparse matrices. The tradeoff is justified for improved robustness in adversarial settings.

6. Limitations, Open Problems, and Theoretical Frontiers

While HM-PCA offers strong guarantees under partition-level contamination, several limitations and considerations remain:

• The inversion and spectral steps can be computationally intensive for large $p$ and $m$.
• The robust gain in eigenvalue ordering is specific to scenarios where contamination is limited to a fraction of partitions. Arbitrarily structured adversarial noise may still present challenges.
• The matrix harmonic mean does not always admit closed-form expressions for more than two blocks; efficient algorithms rely on iterative or spectral approaches.
• Broad extension to other aggregation/statistical learning tasks is an open direction, with the partition-aggregation principle constituting a general strategy for robust distributed inference (Hung et al., 15 Oct 2025).

7. Connections to Asymmetric Norm PCA and Generalized Matrix Means

Methodological developments in PCA under asymmetric loss functions and tail-sensitive objectives (Tran et al., 2014) complement HM-PCA by targeting tail risk and extreme-event structure. While these methods are not directly couched in matrix harmonic mean aggregation, there is a conceptual alignment: both approaches aim to control influence from atypical variations (in tails or contamination) and achieve robust low-rank representations. The iterative reweighted least squares and expectile/quantile-based PCA routines may be interpreted, in a generalized sense, as sharing the robustness-design goals that motivate HM-PCA.

Furthermore, the β-mean framework (Jou et al., 1 Oct 2024) unifies arithmetic, geometric, and harmonic means, enabling flexible adaptation to data properties. The robustness ordering achieved by HM-PCA is distinguished by an infinite tolerance threshold for eigenvalue perturbations (order invariance under any local contamination), not matched by geometric mean (finite tolerance) or arithmetic mean (fragile tolerance).

---

Harmonic Mean PCA establishes a rigorous, theoretically substantiated approach to robust and distributed subspace estimation. By leveraging reciprocal aggregation of covariance spectra, it preserves principal structure under contamination and achieves high-efficiency estimation when data are clean. Table-driven comparison with geometric and arithmetic mean-based frameworks reveal its superiority for applications demanding robust ordering of principal components. Its generalization within the $\phi$-PCA and matrix β-mean paradigms further supports extensions beyond PCA across statistical and machine learning practice.