Robust computation of linear models by convex relaxation

Gilad Lerman; Michael McCoy; Joel A. Tropp; Teng Zhang

Source-linked AI summary

Robust computation of linear models by convex relaxation

Gilad Lerman, Michael McCoy, Joel A. Tropp, Teng Zhang

arXiv:1202.4044v2Feb 18, 2012cs.ITstat.COstat.ML

TL;DR

The paper addresses robust low-dimensional modeling when inliers are noisy and outliers may be numerous or adversarial. It introduces REAPER, a convex relaxation based on orthoprojectors, together with an efficient solver and a theory of approximation. The analysis and experiments support reliable recovery or approximation of linear structure under the paper’s data models.

Problem
Noisy inliers and substantial outliers can make PCA unreliable, creating a need for robust low-dimensional linear modeling.
Method
REAPER parameterizes subspaces with orthoprojectors, replaces their nonconvex set by its convex hull, and solves the resulting convex program with an iterative algorithm.
Results
When the stability statistic is large, REAPER can effectively approximate the inlier subspace despite noisy inliers; under the Haystack Model, positive stability gives perfect recovery.
Takeaways & Limitations
REAPER provides a theoretically analyzed and computationally implementable approach for finding low-dimensional linear structure in data with many outliers.
Takeaways & Limitations
The theoretical analysis assumes an In & Out Model with inliers near a fixed low-dimensional subspace and arbitrarily arranged outliers.

Abstract

from arXiv · show

Consider a dataset of vector-valued observations that consists of noisy inliers, which are explained well by a low-dimensional subspace, along with some number of outliers. This work describes a convex optimization problem, called REAPER, that can reliably fit a low-dimensional model to this type of data. This approach parameterizes linear subspaces using orthogonal projectors, and it uses a relaxation of the set of orthogonal projectors to reach the convex formulation. The paper provides an efficient algorithm for solving the REAPER problem, and it documents numerical experiments which confirm that REAPER can dependably find linear structure in synthetic and natural data. In addition, when the inliers lie near a low-dimensional subspace, there is a rigorous theory that describes when REAPER can approximate this subspace.

1 Introduction

The paper develops REAPER, a convex method for robustly fitting low-dimensional linear models when data contain noise and outliers. It parameterizes subspaces with orthoprojectors, relaxes their feasible set to a convex hull, and supports the approach with theory, algorithms, and experiments.

Motivation: Noisy data with substantial outliers can undermine PCA, motivating robust low-dimensional modeling.The paper situates this problem across computer vision, machine learning, bioinformatics, and related applications.
Validation: The paper combines rigorous performance analysis with an efficient numerical algorithm and experiments on synthetic and natural data.The experiments indicate that REAPER can be more effective than competitors for certain robust linear-modeling problems.
Convex formulation: REAPER replaces the nonconvex orthoprojector constraint with a larger convex set, enabling efficient convex optimization for robust linear modeling.Its relaxation uses symmetric matrices with eigenvalues in [0,1], equivalently the convex hull of rank-d orthoprojectors.
Theory: REAPER’s analysis is based on the geometry of the convex hull and the objective’s limited descent directions near a good inlier projector.These observations support the claim that optimization cannot move far from a suitable inlier subspace while decreasing the objective.
Model representation: Subspaces are represented by orthoprojectors, and the nearest rank-d orthoprojector to a REAPER solution yields the fitted linear model.The range of the resulting projector often fits inliers well even in the presence of many outliers.

2 Theoretical Analysis of the REAPER Problem

The theoretical analysis models noisy inliers near a target subspace alongside arbitrary outliers, then identifies geometric conditions under which REAPER approximates that subspace. Its guarantees show exact recovery without inlier noise, controlled sensitivity to noise, and dependence on outlier alignment.

2.1 A Deterministic Data Model: The In & Out Model places inliers near a fixed low-dimensional subspace while allowing outliers anywhere in the ambient space.The target is to approximate the subspace containing the inliers.
2.2 Summary Parameters for the In & Out Model: REAPER’s effectiveness depends on inlier evidence for the target subspace and limited linear structure among the outliers.The analysis summarizes these properties using permeance, total inlier residual, and alignment statistics.
2.2 Summary Parameters for the In & Out Model: A large stability statistic indicates that REAPER can approximate the target subspace even when the inliers are noisy.Stability increases when inliers support the subspace strongly and outliers provide relatively little distracting structure.
2.3 Performance of REAPER with Deterministic Data: When all inliers lie in the target subspace and stability is positive, REAPER recovers that subspace without error.The guarantee states that the recovered projector equals the target projector.
2.3 Performance of REAPER with Deterministic Data: As inliers move away from the target subspace, approximation error increases linearly with noise and proportionally to S(L)^−1.Thus, large stability makes inlier noise have a small effect on the approximation error.
2.3 Performance of REAPER with Deterministic Data: Adversarial outliers require very large permeance when they cluster along a one-dimensional subspace orthogonal to the target.The theory also controls all principal angles between the computed and target subspaces.
2.3 Performance of REAPER with Deterministic Data: REAPER is nearly capable of finding the oracle subspace that minimizes total unsquared inlier residuals, a criterion less sensitive to large errors than PCA’s squared residuals.This explains the method’s intended robustness advantage over PCA.
2.4 Proof of Theorem 2.1: The proof handles noisy inliers by replacing the original objective with a nearby one and combining perturbation bounds with uniqueness for the noiseless problem.This yields a bound relating the recovered projector to the target projector through the residual and stability terms.

3 Theoretical Example: The Haystack Model

The Haystack Model studies low-dimensional inliers with unstructured outliers, using stability analysis to characterize when REAPER recovers or approximates the target subspace. Theory and experiments show that recovery improves with inlier sampling and remains effective with noise and many outliers.

The Haystack Model represents low-dimensional inliers and totally unstructured outliers in a generative random setting.
The stability statistic grows linearly with the inlier sampling ratio and decreases linearly with the outlier sampling ratio.
A sufficient condition for exact recovery is that the inlier sampling ratio be at least a constant multiple of the outlier sampling ratio.Because Haystack inliers lie in the target subspace, positive stability yields perfect recovery.
Exact recovery experiment: Figure 3.1 shows a linear trend between the 50% empirical success outlier ratio and the inlier sampling ratio for d = 1 and d = 10.The experiments use ambient dimension D = 100.
Noisy inliers: Under noisy inliers, the stability relationship remains unchanged, while the inlier residual becomes nonzero and contributes to approximation error.
Noisy inliers: REAPER occupies the low-error region over more of the tested parameter regime than PCA under the Noisy Haystack Model.The blue region denotes mean error below 10% of the maximum possible error.

4 An Iterative Reweighted Least-Squares Algorithm for REAPER

REAPER is solved with an IRLS method that repeatedly performs weighted least-squares updates, each obtainable through spectral computation. The method has a convergence guarantee, practical linear convergence, and per-iteration cost O(ND^2).

IRLS replaces the REAPER semidefinite program with a sequence of weighted least-squares problems.Each weighted subproblem has a closed-form solution obtained from a single SVD of the weighted data.
The weighted least-squares subproblem is solved using spectral computation and a water-filling step enforcing 0 ≼ P ≼ I.
The weights emphasize observations explained well by the current model, while regularization prevents any single point from gaining undue influence.
Under the stated nondegeneracy assumption, IRLS converges to a feasible point whose objective value is bounded relative to an optimal REAPER point.
Computational cost: Each IRLS iteration requires O(ND^2) arithmetic operations, and practical linear convergence requires O(log(1/η)) iterations for error η.
Empirical convergence: Empirical experiments indicate linear convergence, with challenging synthetic instances requiring as many as 200 iterations and natural-data experiments usually requiring about 20.

5 Numerical Example: REAPER Applied to an Image Database

The experiment evaluates REAPER and related robust methods on face images mixed with random images, emphasizing preprocessing, computational tractability, and out-of-sample generalization. S-REAPER produces the clearest projections and the strongest measured generalization among the compared methods.

Algorithms: S-REAPER usually outperforms REAPER after spherization, while the method is computed using IRLS with weighted least-squares subproblems solved by a single SVD.The implementation uses fixed regularization and stopping tolerances before postprocessing the optimal point.
Practical procedure: The recommended pipeline optionally centers and spherizes observations, solves REAPER, and extracts a d-dimensional invariant subspace from the optimal point.Spherization is motivated by sensitivity to large-magnitude outliers, while the optimal matrix may have rank exceeding d.
Experimental setup: The experiment tests robust linear models on 32 in-sample face images mixed with 467 grayscale Caltech101 images, reserving 32 additional faces for testing.Images are centered using the Euclidean median and modeled in nine dimensions.
Image projections: All methods produce recognizable in-sample face projections, but robust methods capture out-of-sample facial features better than PCA, with S-REAPER producing the clearest images.The figure compares projections of in-sample faces, random images, and two out-of-sample faces.
Generalization: S-REAPER is the dominant method under ordered out-of-sample distances, demonstrating better generalization than the other robust subspaces.Points below the 1:1 line indicate closer robust-model projections than PCA projections.

6 Related Work

The related-work discussion contrasts classical robust modeling, convex demixing, and projector-based relaxations in terms of computational tractability, guarantees, and problem fit. It positions REAPER as a tractable alternative with theoretical and numerical advantages under the considered inlier–outlier model.

Classical methods: Classical robust linear-modeling methods often have strong robustness properties but involve intractable optimization or algorithms designed for different estimation problems.RANSAC can eventually identify a model under assumptions but lacks guarantees on the required number of iterations.
Classical methods: Spherical PCA is computationally practical and recommended as a reliable classical robust method, but it has limited theoretical guarantees.It robustly centers and rescales observations before applying standard PCA.
Theoretical comparison: For the In & Out Model, REAPER has an exact-recovery guarantee stronger than the cited demixing result in the d ≪ D regime, with analogous advantages for stability.The comparison is stated under the model assumptions described in the passage.
REAPER’s position: REAPER tends to be significantly faster than the cited ADMM implementations in limited numerical experiments, while its analysis removes an oracle condition and strengthens exact-recovery guarantees.The present work also extends earlier formulations by treating d as a known target dimension and tightening the projector relaxation.

A Theorem 2.1: Supporting lemmas

The appendix establishes the technical results underlying Theorem 2.1 while retaining the assumptions and notation of the In & Out Model.

Proof organization: The appendix proves the two technical results that support Theorem 2.1.These proofs use the notation and assumptions introduced for the In & Out Model.
Proof organization: The appendix maintains the assumptions and notation from the In & Out Model and Theorem 2.1 throughout.
Proof organization: The supporting lemmas provide the technical foundation for Theorem 2.1.

A.1 Controlling the Size of the Perturbation

The appendix controls perturbations around a target projector by comparing inlier and outlier contributions to directional derivatives and bounding the resulting matrix norms. These bounds establish the stability-related ingredients used in the theorem’s proof.

Perturbation control: The perturbation analysis compares an objective with a perturbed objective whose difference depends only on the inliers.The comparison uses triangle inequalities and identifies the inlier residual term through the Schatten 1-norm.
Directional derivative: The directional-derivative lower bound separates the contributions of outliers and inliers.The outlier term is controlled by the alignment statistic, while the inlier term is controlled by the permeance statistic.
Directional derivative: Combining the outlier and inlier bounds yields the stability statistic S(L), which controls the directional derivative near the target subspace.The proof identifies the relevant bracket using d = dim(L).
Inlier bound: The inlier lower bound is established by normalizing the perturbation, using a thin SVD, and minimizing a concave expression over a convex set.The argument reduces the minimization to an extreme point and uses unit vectors contained in L.
Norm comparison: The proof relates ∥ΔΠ_L∥S1 to ∥Δ∥S1 by block decomposition and the feasible-set constraints 0 ≼ P ≼ I and tr(P) = d.Positive-semidefinite block relations and norm inequalities complete the comparison.

B Analysis of the Haystack Model

The appendix derives exact-recovery conditions for the Haystack Model by combining probabilistic bounds with the deterministic recovery theorem. The simplified Haystack result follows from the more detailed theorem through numerical inequalities.

Theorem B.1 establishes exact-recovery conditions for datasets drawn from the Haystack Model.It applies to an arbitrary d-dimensional subspace L of R^D and uses the stability statistic.
The analysis studies the probabilistic behavior of the permeance and alignment statistics.These bounds are the probabilistic ingredients needed for the deterministic recovery result.
The stability bound in Theorem B.1 holds except with probability 3.5e−c2d/2.
Theorem 3.1 is obtained by applying numerical inequalities to weaken the bound from Theorem B.1.The proof invokes convexity and standard inequalities for sums of squared terms.

B.1.1 Tools for the Permeance Statistic

The permeance analysis separates a centered fluctuation from a mean term and controls the fluctuation using Gaussian concentration. The resulting bounds rely on symmetrization, comparison, and closed-form Gaussian moments.

The permeance statistic is estimated for data drawn from the Haystack model using probability inequalities.
The Gaussian calculation uses E|⟨u, gi⟩| for a half-Gaussian random variable and bounds the expectation by a root-mean-square expression.
The mean calculation uses Rademacher symmetrization and a comparison theorem with the absolute-value function.
The centered fluctuation is dominated by the centering term when n ≫d.
The fluctuation term is controlled with a concentration inequality for Lipschitz functions of independent Gaussian variables.The relevant function has Lipschitz constant √n under the Frobenius norm.

B.1.2 Tools for the Alignment Statistic

The alignment analysis combines Gaussian-matrix and sphere-valued random-vector bounds to control the outlier contribution. These estimates are then combined with the permeance bound to obtain a stability guarantee.

The alignment analysis requires a tail bound for the maximum singular value of a Gaussian matrix.The bound is described as a consequence of Slepian’s lemma.
Sphere-valued random matrices are analyzed by comparing their norms with those of Gaussian matrices through χ-distributed variables.
In the Haystack model, the inlier and outlier sampling ratios are ρin := Nin/d and ρout := Nout/D.
The inliers are drawn from a centered Gaussian distribution on L, while the outliers are independent centered Gaussian vectors in R^D.
The permeance and alignment bounds are combined to establish a stability-statistic bound except with probability 3.5e−c2d/2.

C Analysis of the IRLS Algorithm

The appendix analyzes the IRLS algorithm by verifying its weighted least-squares subproblem and relating its limit to the REAPER optimum.

The IRLS analysis first verifies that Algorithm 4.1 reliably solves the weighted least-squares subproblem.
It then argues that IRLS converges to a point near the true optimum of the REAPER problem.

C.1 Solving the Weighted Least-Squares Problem

This section gives an explicit solution to the weighted least-squares problem by diagonalizing the weighted sample covariance matrix and constructing an optimal feasible matrix. The proof verifies optimality separately according to whether the covariance rank is at most or greater than the target dimension.

Scope and conventions: The lemma’s construction is valid throughout the range 0 < d < D, with the convention 0/0 := 0.The result applies to the full stated dimension range rather than only the algorithm’s presented parameter regime.
Solution construction: The weighted least-squares solution is constructed from the eigenvalue decomposition of the weighted sample covariance matrix.The formulation assumes 0 < d < D and assigns nonnegative weights to the observations.
Rank-dependent cases: When rank(C) ≤ d, the constructed solution attains objective value zero and is therefore optimal.The nonzero eigenvalues of C lie among the first d eigenvalues, whose corresponding coefficients make the residual vanish.
Rank-dependent cases: When rank(C) > d, optimality follows by showing that the objective’s directional derivative is nonnegative in every feasible perturbation.Feasible perturbations satisfy 0 ≼ P⋆+∆ ≼ I and tr∆ = 0.
Optimality proof: A threshold θ is selected from the eigenvalues of C, and the associated eigenspace decomposition establishes the required derivative inequality.The proof uses the subspace spanned by eigenvectors whose eigenvalues do not exceed θ, together with positive- and negative-semidefinite compressions.

C.2 Convergence of IRLS

This section analyzes IRLS by relating its iterates to a regularized convex objective. It shows that the iterates converge to the regularized optimizer and that the resulting point is nearly optimal for the original REAPER problem.

Convergence target: IRLS iterates converge to a matrix Pδ associated with the optimizer Qδ of a regularized REAPER objective.The analysis uses Q = I−P and characterizes the limit through a constrained regularized problem.
Regularized objective: The regularized objective is continuously differentiable, and under the stated observation assumption it is strictly convex, making Qδ well defined.The assumption excludes observations lying in the union of two strict subspaces of R^D.
Proof mechanism: IRLS is analyzed as a generalized Weiszfeld method whose iterates are governed by auxiliary quadratic and potential functions.The proof relates the smooth quadratic surrogate G to the regularized objective F and uses the resulting monotonicity properties.
Algorithmic implication: The objective values decrease along the iterations, providing the stopping criterion used by the algorithm.The monotonicity relation gives α(k+1) = F(Q(k+1)) ≤ F(Q(k)) = α(k).
Near-optimality: The limiting point is nearly optimal for REAPER, with F0(Qδ)−F0(Q⋆) ≤ 1/2 δ|X|.This bound follows by comparing the regularized and unregularized objectives at their respective optimizers.