Two Proposals for Robust PCA using Semidefinite Programming

Michael McCoy; Joel Tropp

Source-linked AI summary

Two Proposals for Robust PCA using Semidefinite Programming

Michael McCoy, Joel Tropp

arXiv:1012.1086v3Dec 6, 2010stat.CO

TL;DR

Outliers make variance-based PCA unreliable, motivating robust alternatives. The paper proposes MDR and LLD, both based on semidefinite programming, and reports numerical evidence for their effectiveness while identifying limitations in multi-component rounding and basis invariance.

Problem
Variance-based PCA is highly sensitive to outliers, while many robust-statistics approaches rely on intractable optimization or lack a principled foundation.
Method
The paper proposes MDR, which uses mean absolute deviation with SDP relaxation and randomized rounding, and LLD, which separates data into low-leverage and corrupted components.
Results
MDR achieves optimality ratios of 0.99999, 0.99992, and 0.97253 for its first three components in one experiment, while outperforming standard PCA on the final outlying observations.
Takeaways & Limitations
The experiments support MDR and LLD as useful robust-PCA techniques for data containing outlying or corrupted observations.
Takeaways & Limitations
The paper reports that multi-component MDR rounding is not competitive with projection pursuit, and that one LLD pruning procedure does not necessarily produce a low-leverage solution.

Abstract

from arXiv · show

The performance of principal component analysis (PCA) suffers badly in the presence of outliers. This paper proposes two novel approaches for robust PCA based on semidefinite programming. The first method, maximum mean absolute deviation rounding (MDR), seeks directions of large spread in the data while damping the effect of outliers. The second method produces a low-leverage decomposition (LLD) of the data that attempts to form a low-rank model for the data by separating out corrupted observations. This paper also presents efficient computational methods for solving these SDPs. Numerical experiments confirm the value of these new techniques.

1. Introduction

The paper develops two semidefinite-programming approaches to robust PCA: MDR replaces variance with a less outlier-sensitive spread measure, while LLD separates corrupted observations from a low-leverage component. It also provides efficient SDP algorithms and numerical experiments.

PCA identifies high-variation directions for dimensionality reduction, but variance-based methods are highly sensitive to outliers.
MDR replaces variance with mean absolute deviation and approximates the resulting non-convex optimization problem using an SDP and randomized rounding.
LLD decomposes the data into a low-leverage matrix and a matrix of corrupted observations using a semidefinite program.
The paper develops efficient algorithms for solving the SDPs and reports numerical experiments supporting the effectiveness of both techniques.

2. Maximum Mean Absolute Deviation Rounding

The paper motivates robust PCA by replacing variance with robust scales, focusing on mean absolute deviation and developing a provably approximate algorithm despite computational intractability.

Robust PCA replaces variance with a scale intended to reduce sensitivity to outliers, with mean absolute deviation as the paper’s chosen alternative.
The paper’s first method is a specific projection-pursuit PCA formulation using the mean absolute deviation scale.
The formulation is computationally intractable, but the paper develops an algorithm that provably approximates its solution.

2.2. Projection-Pursuit PCA.

Projection-pursuit PCA finds directions maximizing a chosen robust scale and adds orthogonality constraints to obtain further components. With standard deviation, it reduces to classical PCA.

A robust component is the direction maximizing a robust scale of the projected data.
Further robust components are constructed by imposing orthogonality constraints.
Projection-pursuit PCA reduces to classical PCA when the scale is standard deviation.
Implementing the method requires a procedure for finding the first component, with orthogonality enforced separately.

2.3. PP-PCA with the MD Scale is NP-Hard.

Using the mean absolute deviation scale makes the projection-pursuit PCA optimization NP-hard, limiting existing general-purpose algorithms that seek both accuracy and efficiency.

Replacing the ℓ2 norm with the ℓ1 norm makes the projection-pursuit optimization NP-hard.
Accurately computing the relevant operator norm in polynomial time would imply P = NP.
Because existing projection-pursuit algorithms claim generality across scales, they cannot provide both accurate and efficient solutions for this problem.

2.4. Approximating the ℓ2 →ℓ1 Norm using Randomized Rounding.

The section relaxes the NP-hard ℓ2→ℓ1 norm computation to a semidefinite program and uses randomized rounding to recover an approximate unit vector.

Relaxation and rounding: The ℓ2→ℓ1 norm is NP-hard to compute, but can be approximated efficiently using semidefinite-programming and randomized-rounding ideas.The approximation combines the little Grothendieck theorem with techniques from Nesterov, Burer and Monteiro, and a new factorization step.
Relaxation and rounding: Relaxing the rank-one constraint Z = yy∗ to Z ≽ 0 produces the semidefinite relaxation.The original formulation imposes Z = yy∗ and [Z]ii = 1; the relaxation retains the diagonal constraints while replacing rank one with positive semidefiniteness.
Approximation guarantee: Grothendieck’s inequality bounds the SDP approximation factor, and the resulting bound is asymptotically best possible.The cited inequality is described as asymptotically tight for the relevant class of examples.
Relaxation and rounding: The algorithm solves the SDP, samples randomized sign vectors, and rounds an optimal solution Z⋆ back to a unit-norm vector v⋆.The rounding procedure targets a vector satisfying ∥Xv⋆∥≈∥X∥2→1.

2.5. The MDR Algorithm.

The MDR algorithm uses SDP-based randomized rounding to approximate the robust principal component objective, with guarantees controlled by the number of trials and practical ratios often near 0.95 or higher.

Algorithm and guarantee: MDR rounds an SDP solution repeatedly and selects the best random outcome to reduce the probability of an inadequate approximation.The method’s novel Step 2(b) computes v from the sampled sign vector, while Step 3 chooses the best trial.
Algorithm and guarantee: Theorem 2.2 guarantees α⋆≥∥X∥2→1 and provides a probability-dependent approximation bound for the returned vector.The guarantee is parameterized by the repetition count K and a threshold θ < 1.
Algorithm and guarantee: K = 94 yields ∥Xv⋆∥1 > 0.75∥X∥2→1 with probability at least 0.999.The approximation ratio ρ = ∥Xv⋆∥1/α⋆ is used to measure solution quality.
Algorithm and guarantee: In practice, MDR typically achieves an approximation ratio of 0.95 or higher, while the theoretical limit 2/π > 0.79 is asymptotically tight.The practical observation does not imply that the theoretical analysis is loose.
Computation: Interior-point methods solve the SDP in polynomial time in principle, but can be too slow for large-scale problems.The paper therefore prefers a Burer–Monteiro approach for practical computation.

2.6. Implementation of Algorithm 1.

The implementation replaces the SDP matrix with a lower-dimensional factor, optimizes the resulting nonconvex objective, and iteratively applies the procedure to obtain multiple orthogonal robust components.

Factorized SDP implementation: The Burer–Monteiro implementation factors the semidefinite matrix as Z = RR∗ and optimizes directly over R.This avoids explicitly forming Z and converts the diagonal constraints into unit-norm row constraints on R.
Factorized SDP implementation: Row normalization is incorporated into an unconstrained nonconvex objective, which is maximized using a conjugate-gradient algorithm.The implementation uses the Hager–Zhang conjugate-gradient method.
Factorized SDP implementation: Restricting R to an n×k matrix with k = O(√n) suffices to solve the SDP exactly.Under a mild technical condition, any local minimum of the factorized problem yields a global minimum of the original SDP.
Multiple components: To obtain the kth component, the method restricts X to the subspace perpendicular to previously computed components and reapplies Algorithm 1.Householder reflections provide a more numerically stable implementation and reduce the restricted matrix dimension.

2.7. Extending the Rounding to Multiple Components.

The paper extends the rounding formulation to multiple orthogonal components through a joint SDP, but finds projection pursuit with sequential rounding more competitive in practice.

Joint formulation: The multiple-component formulation imposes orthogonality through ⟨vi, vj⟩ = δij, so its optimum consists of orthogonal unit vectors.When T = 1, the formulation reduces to the single-component problem.
Joint formulation: The joint nonconvex problem is approximated by a semidefinite relaxation whose optimal value is within a logarithmic factor of the original optimum.The relaxation’s rounding does not itself produce orthogonal vectors, requiring an additional orthogonalization step.
Joint formulation: Algorithm 2 outputs a p × T matrix V⋆ with orthogonal columns by computing an SVD and retaining its first T right-singular vectors.The extracted columns are taken from V in the SVD P⋆ = UΣV∗.
Limitation: The joint SDP does not appear competitive with projection pursuit, whose feasible vectors typically achieve a larger objective value than rounded vectors followed by orthogonalization.The paper identifies improved rounding for this relaxation as a direction for further research.

3. The Low-Leverage Decomposition

The low-leverage decomposition (LLD) uses a semidefinite program to separate corrupted observations from a low-leverage surrogate, whose singular vectors define robust components. Its guarantees and parameter behavior clarify both the method’s statistical interpretation and practical limitations.

3. The Low-Leverage Decomposition: LLD decomposes the data matrix into a low-leverage matrix and a matrix of corrupted observations.The decomposition is obtained with a semidefinite program and interprets the two components as surrogate data and identified corruption.
3. The Low-Leverage Decomposition: The optimization penalizes the nuclear norm of P to promote low rank and the row-group norm of C to promote group sparsity.The constraint P + C = X enforces the decomposition of the input data matrix.
3. The Low-Leverage Decomposition: Robust components are the right singular vectors of the surrogate P⋆, while C⋆ represents identified corruption.The optimal P⋆ additionally forms a low-leverage set of observations.
3. The Low-Leverage Decomposition: The method’s optimality conditions characterize feasible solutions through a matrix Q satisfying subgradient conditions for the two norms.These conditions support the leverage-score bound established for the optimal decomposition.
3.1.1. Leverage scores.: The optimal decomposition has leverage scores bounded above by γ^2, giving γ a statistical interpretation as a user-specified leverage bound.The same bound also implies an upper bound on the rank of P⋆ through the trace of its hat matrix.
3.2. The Choice of γ.: Forcing low rank can produce dense corruption when data are not truly low rank, and γ choices trade off low-rank modeling against outlier identification.The paper reports that γ somewhat less than p/n gives a very good low-rank model but performs poorly for outlier identification; large-scale instances also require methods beyond generic interior-point software.

4. Previous Work

Previous robust PCA methods include heuristic, nonconvex, and convex approaches with different guarantees and invariance properties. The paper positions its proposals against these methods, emphasizing tractability, theoretical support, and basis-invariance limitations.

4. Previous Work: Many robust PCA proposals lack rigorous algorithms, pairing mathematical formulations with heuristics or algorithms without clear formulations.Convex approaches are described as unusual in the robust PCA literature.
4.1. Projection-pursuit PCA: MDR is a projection-pursuit PCA method, but computing the underlying problem is NP-hard and existing algorithms generally rely on heuristics.The paper notes that finite-direction and Monte Carlo methods do not provide the same tractability guarantees.
4.2. A convex approach.: The N + L1 approach decomposes data into low-rank and sparse components using a tractable convex program with exact-recovery conditions under stated assumptions.Its low-rank and sparsity penalties are based on matrix norms.
4.2. A convex approach.: N + L1 is not invariant to changes in the observation basis because its matrix norm is not invariant under orthogonal transformations.The paper argues this limitation matters when PCA’s goal is to identify a more effective coordinate basis.
4. Previous Work: Spherical PCA normalizes each observation before applying standard PCA and is used as a baseline because of its practical performance and ease of implementation.Maronna’s study is cited as motivating its use as a default robust PCA choice.
4. Previous Work: Other approaches include robust covariance estimation, nonconvex iteratively reweighted methods, and random trimming with theoretical guarantees in a high-dimensional regime.These methods differ in their assumptions, optimization procedures, and available guarantees.

5. Numerical Experiments

Numerical experiments compare MDR and LLD with standard PCA and other robust methods on environmental, iris, regression-surface, and movie-rating data. The results show that robust methods generally increase projected spread, better explain contaminated data, and alter component-based rankings, with method-specific trade-offs.

Experimental design: The experiments compare the proposed MDR and LLD methods with standard PCA and other robust methods across several data sets and tasks.The study includes projections of Oslo NO2 and iris data, a multiple-component regression-surface experiment, and movie-rating component rankings.
Oslo NO2: Every robust method produces a larger IQR than PCA on the Oslo NO2 data, with MDR largest and LLD smallest among robust methods.Except for N + L1, the methods identify directions with relatively many outliers, indicating heavy tails in the data.
Oslo NO2: 0.978 is the reported approximation ratio for the top MDR component, described as near optimal.
Iris: On iris data, LLD, sphPCA, and N + L1 each achieve an IQR at least 3 times that of PCA, while LLD and sphPCA appear as effective as setosa-only PCA.These methods do not clearly distinguish the three populations, whereas MDR discriminates among them and has an IQR 1.95 times that of setosa along dominant PCA.
Regression surfaces: Each robust method explains the data better than PCA for more than 75% of points in the regression-surface experiment.MDR and LLD become worse than PCA only after 95% of the data, and LLD dominates through the latter part; MDR explains final outliers better than the other methods.
MovieLens: Robust components change movie rankings while preserving sign consistency: all methods rank GoodFellas first on component one, and robust methods rank Mommie Dearest first on component two.The paper reports that robust methods assign more importance to discriminating movies, while PCA assigns the smallest weight to nearly every movie.

Appendix A. Proof of Theorem 2.2

Appendix A proves Theorem 2.2 by establishing bounds for the semidefinite-programming relaxation and the randomized rounding used in Algorithm 1. The proof combines operator-norm relations, independent rounding trials, and a probability bound based on the Paley–Zygmund inequality.

SDP bound: The proof lower-bounds the optimization value using ∥XX∗∥∞→1 and invokes the SDP relaxation discussion to establish the corresponding guarantee.
Probability bound: The probability step applies a bounded-variable Paley–Zygmund variant to lower-bound the chance that ∥Xv⋆∥1 exceeds a fraction of its expectation.
Theorem statement: Theorem A.4 states that Algorithm 1 outputs α⋆ satisfying α⋆≥∥X∥2→1, together with an additional inequality for θ∈[0,1].The theorem assumes X is an n×p matrix and K rounding trials.
Randomized rounding: For each randomized sign vector y, the proof constructs v=X∗y/∥X∗y∥2 and uses this choice to derive the needed relation for the rounded vectors.
Trial aggregation: Independence of the rounding vectors v(k) lets the proof translate per-trial bounds on ∥Xv(k)∥1 into a bound for the selected vector v⋆.