Sparse outlier-robust PCA for multi-source data

1. Problem Statement

Principal Component Analysis (PCA) is a standard tool for dimensionality reduction, but traditional PCA faces three significant limitations in modern data science contexts:

1. Lack of Sparsity: Standard PCA loadings often involve all variables, making interpretation difficult in high-dimensional settings. Sparse PCA (SPCA) addresses this but is typically applied to single datasets.
2. Multi-Source Data: Many scientific applications involve multiple related datasets (e.g., time series, spatial groups, or different demographic subgroups). Applying PCA globally ignores local group-specific structures, while applying it locally ignores global connections. There is a need for a method that simultaneously captures global patterns (shared across sources) and local patterns (specific to a source).
3. Sensitivity to Outliers: Classical PCA relies on non-robust covariance estimators (like the sample covariance), which are heavily distorted by outliers. While robust PCA methods exist (e.g., ROSPCA), they generally do not handle the structured sparsity required for multi-source data.

The authors aim to fill this gap by developing a method that performs joint sparse PCA on multiple data sources while being robust to outliers.

2. Methodology

The proposed method, Sparse Multi-Source Robust PCA, integrates three core components: a structured sparsity penalty, a robust covariance plug-in, and an efficient optimization algorithm.

A. Objective Function and Structured Sparsity

The method formulates the estimation of the $k$-th principal component loading matrix $V^k$ (where columns represent sources and rows represent variables) as an optimization problem. The objective minimizes the negative variance explained while penalizing complexity:

$$\min_{V} -\sum_{i=1}^N v_{\cdot i}' \hat{\Sigma}_i v_{\cdot i} + \eta \gamma \sum_{j=1}^p \sum_{i=1}^N |v_{ji}| + \eta(1-\gamma)\sqrt{N} \sum_{j=1}^p \|v_{j\cdot}\|_2$$

• $\hat{\Sigma}_i$: Robust covariance estimator for source $i$.
• $\eta$: Controls the overall degree of sparsity.
• $\gamma \in [0, 1]$: Balances between local sparsity ($\gamma=1$, element-wise L1 penalty) and global/structured sparsity ($\gamma=0$, group-wise L2 penalty).
  • Local Sparsity: Selects specific variables for specific sources.
  • Global/Structured Sparsity: Forces the same variable to be zero (or non-zero) across all sources simultaneously, leveraging the interconnection between datasets.
• Constraints: Orthogonality constraints are enforced per source ($v_{\cdot i} \perp v_{\cdot l}$) to ensure distinct principal components.

B. Robustness via ssMRCD Plug-in

To handle outliers, the method uses the spatially smoothed Minimum Regularized Covariance Determinant (ssMRCD) estimator as a plug-in for $\hat{\Sigma}_i$ and $\hat{\mu}_i$.

• Mechanism: The ssMRCD selects a subset of "clean" observations for each source but smooths the covariance estimates across sources using a weight matrix $W$ (representing similarity between sources, e.g., temporal or spatial proximity).
• Benefit: This allows the method to borrow strength from neighboring sources to improve estimation accuracy while maintaining robustness against outliers in individual sources.

C. Optimization Algorithm (ADMM)

The optimization problem is non-convex due to the orthogonality constraints and the specific penalty structure. The authors propose a solution using the Alternating Direction Method of Multipliers (ADMM):

• Decomposition: The problem is split into consensus subproblems involving four variables ($v^{(0)}, v^{(1)}, v^{(2)}, v^{(3)}$) to handle the quadratic term, L1 penalty, and group L2 penalty separately.
• Convergence Enhancements:
  • Starting Values: A heuristic starting value is constructed by averaging the solution of standard PCA ($\eta=0$) and an extreme sparse solution ($\eta \to \infty$).
  • Sign Constraint: To prevent the algorithm from oscillating between opposite vectors (a common issue in non-convex PCA), a constraint $z' B_i v^k \geq 0$ is added to fix the direction.
  • Inexact Minimization: The root-finding subproblem is stopped early to save computational time without sacrificing accuracy.
• Hyperparameter Tuning:
  • $\gamma$ and $\eta$ are selected by maximizing the Trade-off Product (TPO) of standardized sparsity and explained variance.
  • The smoothing parameter $\lambda$ for ssMRCD is selected by minimizing the mean residual norm of the smallest $\alpha$-fraction of data.

3. Key Contributions

1. Novel Framework: First extension of robust PCA to jointly analyze multiple related datasets with structured sparsity (simultaneously detecting global and local patterns).
2. Robustness Integration: Successfully integrates the ssMRCD estimator into a sparse PCA framework, providing protection against outliers in multi-source settings.
3. Algorithmic Efficiency: Development of a tailored ADMM algorithm with specific enhancements (warm starts, sign constraints, inexact minimization) to solve the non-convex, constrained optimization problem efficiently.
4. Theoretical & Practical Validation: Rigorous simulation studies and real-world applications demonstrating superiority over existing methods.

4. Results

Simulation Studies

• Sparsity Recovery: The method successfully recovered both entry-wise and group-wise sparse patterns. It outperformed standard approaches in identifying the correct zero/non-zero structure, especially in complex scenarios with opposing variance directions across sources.
• Outlier Robustness: In the presence of contamination (up to 20% outliers), the proposed ssMRCD-based method significantly outperformed:
  • ROSPCA (Robust Sparse PCA): ROSPCA failed to maintain accuracy when outliers were present, particularly in local contamination scenarios.
  • Non-robust variants: The proposed method maintained low angle errors (subspace estimation) and low orthogonal distances (reconstruction error) even with outliers, whereas non-robust methods degraded rapidly.
• Local Contamination: When only one source was contaminated, the proposed method remained stable for that specific source, whereas ROSPCA's performance on the contaminated source collapsed.

Real Data Applications

1. Weather Analysis (Hohe Warte Station):
   • Data: 64 years of daily weather data (13 variables).
   • Findings: The method identified a global pattern driven by precipitation (PC1) and a seasonal pattern driven by temperature/radiation (PC2). Crucially, PC3 revealed a long-term trend (likely climate change) where loadings for temperature and wind changed over the 64-year period. This evolving local structure was invisible to standard global PCA.
2. Plant Geochemistry (NEXT Project):
   • Data: Elemental concentrations in 3 plant species and 3 organs (19 variables).
   • Findings: The method identified distinct loading patterns for different plant-organ combinations. It successfully differentiated between calcsilicate and mafic rock formations based on plant tissue chemistry. Specifically, Scots Pine bark showed the highest discriminative power for mineral exploration, a finding that would be obscured by global averaging.

5. Significance

This paper addresses a critical gap in multivariate analysis where data is inherently grouped but requires joint modeling.

• Interpretability: By enforcing structured sparsity, the method produces loadings that are easier to interpret, highlighting which variables drive variance globally versus locally.
• Reliability: In an era of "big data" where outliers and data contamination are common, the robustness of the ssMRCD plug-in ensures that the extracted patterns reflect the true underlying structure rather than noise.
• Versatility: The framework is flexible regarding source definition (time, space, demographics) and can be extended to other multivariate techniques like discriminant analysis or graphical modeling.

The proposed method offers a powerful, computationally efficient tool for researchers dealing with complex, multi-source, high-dimensional, and contaminated datasets.