This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at http://dx.doi.org/10.1109/TPAMI.2015.2511743
SUBMITTED TO IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XXX, NO. XXX, MONTH YEAR

1

Scalable Robust Principal Component
Analysis using Grassmann Averages
Søren Hauberg, Aasa Feragen, Raffi Enficiaud, and Michael J. Black
Abstract—In large datasets, manual data verification is impossible, and we must expect the number of outliers to increase withdata size. While principal component analysis (PCA) can reduce data size, and scalable solutions exist, it is well-known thatoutliers can arbitrarily corrupt the results. Unfortunately, state-of-the-art approaches for robust PCA are not scalable. We notethat in a zero-mean dataset, each observation spans a one-dimensional subspace, giving a point on the Grassmann manifold.
We show that the average subspace corresponds to the leading principal component for Gaussian data. We provide a simple
algorithm for computing this Grassmann Average (GA), and show that the subspace estimate is less sensitive to outliers than
PCA for general distributions. Because averages can be efficiently computed, we immediately gain scalability. We exploit robust
averaging to formulate the Robust Grassmann Average (RGA) as a form of robust PCA. The resulting Trimmed Grassmann
Average (TGA) is appropriate for computer vision because it is robust to pixel outliers. The algorithm has linear computational
complexity and minimal memory requirements. We demonstrate TGA for background modeling, video restoration, and shadow
removal. We show scalability by performing robust PCA on the entire Star Wars IV movie; a task beyond any current method.
Source code is available online.
Index Terms—Dimensionality reduction, Subspace estimation, Robust Principal Component Analysis

F

1

I NTRODUCTIONACROSS many fields of science and in many application
domains, principal component analysis (PCA) is one of
the most widely used methods for dimensionality reduction,modeling, and analysis of data. It is a core technique usedthroughout computer vision. While methods exist for scaling
PCA to large datasets, one fundamental problem has yet
to be addressed. Large datasets are often collected either
fully or semi-automatically and are too large to be verified
by humans. As a result, we should expect a significant
number of outliers. In essence, “big data” implies “bigoutliers.” This requires a form of PCA that is robust to
outliers. While there are several solutions that make PCA
robust to outliers [1]–[6], these generally do not scale tolarge datasets. Typically, robust methods for PCA are notscalable and scalable methods for PCA are not robust. Ourcontribution is a novel formulation and a scalable algorithm
for robust PCA that outperforms previous methods, has low
computational complexity, and scales to large datasets.
Our first contribution is to formulate subspace estimation
as the computation of averages of subspaces. Given a zero• S. Hauberg is with the Section for Cognitive Systems at the Technical
University of Denmark.
E-mail: sohau@dtu.dk• A. Feragen is with the Department of Computer Science at the Universityof Copenhagen, Denmark.
E-mail: aasa@diku.dk• R. Enficiaud is with the Software Workshop at the Max Planck Institutefor Intelligent Systems in Tübingen, Germany.
E-mail: raffi.enficiaud@tue.mpg.de• M.J. Black is with the Perceiving Systems Department at the Max Planck
Institute for Intelligent Systems in Tübingen, Germany.
E-mail: black@tue.mpg.de

Fig. 1. Left: A zero-mean dataset represented as a setof points. Center: The same data represented as a setof one-dimensional subspaces. Right: The blue dottedsubspace is the average subspace.

mean dataset x1:N ⊂ RD , we observe that each observation
xn spans a one-dimensional subspace. Mathematically,each of these subspaces can be described as a point on
the Grassmann manifold (the space of subspaces of RD ;see Sec. 3). We then develop an averaging operator on
the Grassmann manifold that formalizes the notion of theaverage subspace spanned by the data; see Fig. 1 for an
illustration. We show how this Grassmann Average (GA)relates to standard PCA and prove that, for Gaussian data,
the subspace found by the GA corresponds to that of standard
PCA. We show both formally and experimentally that GA
is already more robust than PCA as it optimizes an L1 styleenergy. But robustness can be further improved.
Understanding how PCA can be reformulated as thecomputation of averages (of subspaces) leads to a key observation: Robust averages are theoretically well understood andcan be computed efficiently. We leverage this fact to define
the Robust Grassmann Average (RGA), which generalizes
PCA to robust PCA. There are a range of robust averaging
methods, and any of these are appropriate for RGA but some

Copyright (c) 2016 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at http://dx.doi.org/10.1109/TPAMI.2015.2511743
SUBMITTED TO IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XXX, NO. XXX, MONTH YEAR

scale better than others.
In many fields and applications, outliers exist within adata vector while the majority of the vector is perfectly fine.
This is the common situation in computer vision, where, e.g.
a few pixels in an image may be corrupted. We show thatthe robust averaging in RGA can be done at the “element”
level, and we adopt a trimmed average that is both robust
and can be computed efficiently. The resulting Trimmed
Grassmann Average (TGA) can cope with both isolated
outliers within a vector (image) or entire data outliers.
The resulting algorithm is straightforward, theoreticallysound, and computationally efficient. We provide a detailed
evaluation of the method’s robustness and efficiency incomparison with related methods. We focus on problems
involving images and image sequences although the method
is general. We compare with Candes et al. [5] on their dataand show similar or superior results. Furthermore, we showsuperior results to both Candes et al. and De la Torre andBlack [6] on a problem of archival film restoration; here
we are able to provide a quantitative comparison. Finally,
we show how the method scales by applying it to video
modeling problems where all previous robust PCA methods
fail; for example, we compute TGA for all 179,415 framesof Star Wars IV. Our implementation is available online1 ,which makes it easy to apply it to new problems.
This paper extends a previous conference paper [7] witha probabilistic interpretation of subspace averages, a moreelaborate review of the literature, more intuitive illustrations,and a more thorough experimental evaluation.

2

R ELATED W ORK

Consider the problem of estimating a one-dimensional
subspace. PCA finds the subspace in which the projected
observations are as similar as possible to the original
observations [8]qPCA = arg minNXv, kvk=1 n=1

2

kxn − (x|n v) vk ,(1)

where the resulting unit vector, v, spans the estimatedsubspace. This least-squares energy implies that outliers
have a large influence on the results. Indeed, a single outlier
can arbitrarily corrupt the estimated components. We saythat PCA has a break-down point of 0% [2]; i.e. if morethan 0% of the data are outliers we cannot trust the results.
2.1

Robust PCA

Many suggestions have been made for how to make PCA
robust to outliers. The most classic approaches attempt toestimate the covariance matrix in a robust way by assigningweights to each data point, depending on how likely it is tobe an outlier [1]–[4], [9], [10]. In this scenario, the principal
components are defined as the eigenvectors ofΣrobust =NX

αn xn x|n , with αn ≥ 0.

n=1

1. http://http://ps.is.tue.mpg.de/research projects/robust-pca(2)2

Xu and Yuille [4] restrict the weights α1:N to be binary andestimatethese by optimizing Eq. 1 with an added regularizerPν n (1−αn ) that avoids the trivial solution αn = 0. This is
optimized using a mean field approximation, which reduces
the problem to solving a non-linear program. This works
well and provides a classification of the data points into
inliers / outliers, but is computationally demanding.
Feng et al. [3] iteratively estimate the weights α1:N byperforming ordinary PCA and assigning lower weights topoints that contribute substantially to the total projectedvariance. The approach can theoretically handle up to 50%
outliers, which makes it very powerful. Unfortunately, therepeated computation of ordinary principal components
prevents the algorithm from scaling to large datasets.
Rather than estimating weights for each data point,
several authors have suggested replacing the squared errorin (1) with a measure that is less sensitive to outliers[5], [6], [11]. Ding et al. [11] suggest replacing k · k2
with k · k in (1), leading to the so-called R1 -PCA. The
resulting energy can be optimized by iteratively estimating
a robust covariance matrix, which makes the approachcomputationally demanding. Alternatively, Eq. 1 can be
rewritten as the maximization of explained varianceqPCA = arg max

NX(x|n v)2 .

(3)v, kvk=1 n=1

Kwak [12] suggests replacing the square in this expression
with an absolute value to provide a robust energy, though
it is not robust to pixel-level outliers. Kwak provides a
local optimizer, while McCoy & Tropp [13] show that thegeneral problem is NP-hard and provide a semi-definite
approximation algorithm, which is guaranteed to provide
a near-optimal solution. Naor et al. [14] extend this workto estimate multiple components. We show that Kwak’s
approach gives an average subspace, and use this observation
to derive a general pixel-wise robust subspace average.
2.1.1

Robust M–Estimators

In computer vision, outliers often occur at the pixel level
rather than at the image level [5], [6]. De la Torre and Blackphrase robust subspace estimation in the framework of M–
estimators [2] and develop a robust extension of Eq. 1 using
the Geman-McClure error function [15]. This energy is
optimized using gradient descent with an annealing scheme
over the parameters of the Geman-McClure error function.
Their approach uses an initial singular value decomposition(SVD) followed by a series of iterative updates. While each
iteration of this scheme only has linear complexity, thereliance on SVD gives an overall cubic complexity [16].
This makes the algorithm impractical for large datasets.
2.1.2

Convex Programs

The approach from Candes et al. [5] decomposes the data
matrix X ∈ RN ×D into
X = L + S,

(4)

Copyright (c) 2016 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at http://dx.doi.org/10.1109/TPAMI.2015.2511743
SUBMITTED TO IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XXX, NO. XXX, MONTH YEAR

where L is low rank and S is sparse; i.e. it contains the pixellevel outliers. They show that under very broad assumptions
this can be solved as a convex program, and provide aseries of algorithms, where the so-called Inexact ALM [17]method is most commonly used. This algorithm repeatedly
computes an SVD of the data giving it cubic complexity [16].
Again, this complexity renders the approach impractical foranything beyond small-to-medium scale data sizes. This canbe alleviated somewhat through parallel implementations of
SVD [18], but this requires the availability of large compute
farms and does not address the underlying complexity issues.
Netrapalli et al. [19] show that the decomposition can be
performed with a per-iteration rank-reduced SVD, which
only has quadratic. complexity. While an improvement, thisis still impractical for large datasets.

3

Fig. 2. An illustration of the Grassmann manifold
Gr(1, 2) of 1-dimensional subspaces of R2 . Any subspace is represented by a point on the unit sphere withthe further identification that opposing points on thesphere correspond to the same point on Gr(1, D).

2.1.3 Approximation Techniques
The non-linear complexity of the different robust formula- 2.3 Subspaces in Visiontions of PCA has led to the development of a multitude ofapproximation techniques. Both Mackey et al. [20] and Lui Subspace estimation appears throughout computer vision,
et al. [21] propose subsampling algorithms for solving the so it is not surprising that the Grassmann manifold has
decomposition problem of Candes et al. [5]. This improves seen many applications; for a survey see e.g. [25]. In
scalability as they need only solve smaller SVD problems. this paper we consider the problem of estimating a single
While theoretical analysis shows that subsampling provides subspace from data. The more general problem of estimating
good approximations, it is inherently disappointing that to multiple subspaces can also be considered. For example
analyse large quantities of data one should ignore most of Generalized PCA [26] provides an elegant algebraic solution
it. In a similar spirit, Mu et al. [22] propose to replace the to this subspace clustering problem. The numerical solution
data matrix with a low-rank counterpart found through a requires fitting a higher order polynomial to the data, which
random projection during optimization. This again allows tends to be unstable in poor signal-to-noise scenarios.
In the absence of noise, the subspace clustering problem
for computational gains as smaller SVD problems need tocanbe studied by assuming that any data point can be
be computed when parts of the data are not used.
expressed as a linear combination of a few other data points.
This leads to a sparse convex optimization problem, which
2.2 Scalable PCA
can be made robust to noisy data [27]–[29]. The algorithms,
Most PCA implementations either perform eigenvalue however, do not scale to large datasets as they rely on SVD.
decomposition on the data covariance matrix or SVD on thedata matrix [8]. The former only works in low-dimensional
scenarios, where the covariance can easily be stored in 3 T HE G RASSMANN AVERAGE
memory, while the latter is impractical when either the
Our key contribution is to formulate dimensionality reducnumber of observations or the dimensionality is large.
Due to the importance of PCA for data analysis, scalable tion as the estimation of the average subspace of those
implementations are readily available [23], [24]. A common spanned by the data. Informally, if we seek an averaging
choice is the EM PCA algorithm [24], which computes each operation that returns a subspace, then the input to thisoperation must be a collection of subspaces. To formalize
principal component by repeatedly iterating
this idea we need the notion of “a space of subspaces,”
N
X
ṽ
which
is provided by the classical geometric constructionṽ ←(x|n vi−1 )xn
and
vi ←
(5)
kṽk
known
as the Grassmann manifold [30, p. 24]:n=1
until convergence. Here vi denotes the estimate of the Definition 1. The Grassmann manifold Gr(k, D) is theprincipal component at iteration i of the algorithm. This space of all k-dimensional linear subspaces of RD .
can either be derived as an EM algorithm [24] under an
uninformative prior, as a gradient descent algorithm or as a
Assuming we are given a dataset x1:N with the mean
power iteration [16]. One iteration of the algorithm has linear subtracted, then each observation spans a one-dimensional
complexity and the low memory use makes it highly scalable. subspace of RD and, hence, is a point in Gr(1, D). The
Unfortunately, the approach converges slowly and “zig-zags” space of one-dimensional subspaces in RD takes a particuaround the optimum, implying low numerical precision. larly simple form since any such subspace can be represented
This imprecision accumulates when multiple components by a unit vector u or its antipode −u. Thus, Gr(1, D) canare estimated, which practically limits the approach to only be written as the quotient S D−1 /{±1} of the unit sphereestimate a few of the leading components.
with respect to the antipodal action of the group {±1} [30].
Copyright (c) 2016 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at http://dx.doi.org/10.1109/TPAMI.2015.2511743
SUBMITTED TO IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XXX, NO. XXX, MONTH YEAR

We denote points in Gr(1, D) as[u] = {u, −u},(6)

and note that [u] = [−u]; i.e., given an observation, xn , werepresent the spanned subspace as [un ] = ±xn/kxn k, where
the sign is unknown (see Fig. 2).
3.1

The Average SubspaceWith zero-mean data, we now define the average subspacecorresponding to the first principal component. Weighted
averages are commonly defined as the minimizer of theweighted sum-of-squared distances; i.e.
[q] = arg minNX

wn dist2Gr(1,D) ([un ], [v]),(7)[v]∈Gr(1,D) n=1

where w1:N are weights and distGr(1,D) is a distance onGr(1, D). We define this distance through the quotient spaceconstruction of Gr(1, D) as [31, p. 65](8)

dist2Gr(1,D) ([u1 ], [u2 ])= min{dist2S D−1 (v1 , v2 ) | v1 ∈ [u1 ], v2 ∈ [u2 ]},where distS D−1 is a yet-to-be-defined distance on the unit
sphere S D−1 . With this, the subspace distance distGr(1,D) is
merely the shortest spherical distance between any possible
pair of unit vectors representing the subspaces. The sphericaldistance will be chosen to give an efficient algorithm forcomputing averages on Gr(1, D), but first we prove a useful
relationship between averages on Gr(1, D) and those on
the unit sphere S D−1 :Lemma 1. For any weighted average [q] ∈ Gr(1, D)
satisfying Eq. 7, and any choice of q ∈ [q] ⊂ S D−1 , there
exist u1:N ⊂ [u1:N ] ⊂ S D−1 such that q is a weighted
average on S D−1 :q = arg min

NXv∈S D−1 n=1

1) Pick representations u1:N of [u1:N ] that are closest
to the current average estimate q under distS D−1 .
2) Update q to be the spherical average of u1:N .
For this algorithm to be practical, we need the ability tocompute fast spherical averages. With this in mind, we pick1
ku1 − u2 k2 = 1 − u|1 u2 . (10)2Under this distance, the weighted average of data u1:N ⊂S D−1 is given in closed-form by [32, §2.1]
dist2S D−1 (u1 , u2 ) =q=µ(w1:N , u1:N ),kµ(w1:N , u1:N )k

(11)

where
µ(w1:N , u1:N ) =

N
X

wn

!−1 N
X

n=1

wn un

(12)

n=1

denotes the weighted Euclidean average. This is also the
max-likelihood estimator of the mean of the spherical vonMises-Fisher distribution [32]. As we need to normalize thedata, x1:N , to unit length before viewing each observation
as a subspace it is natural to pick subspace weights
wn = kxn k.

(13)Noting that the optimal sign for un is αn = sign(u|n qi−1 ),we arrive at the following iterative scheme for computing
an average subspace:
Algorithm: Grassmann Average (GA)xn
Let un = kx
, and i = 1.
nk
repeatαn ← sign(u|n qi−1 ) ∀n,µ(w1:N , (αu)1:N ),qi ←kµ(w1:N , (αu)1:N )ki ← i + 1,

(14)

until kqi − qi−1 k < .
wn dist2S D−1 (un , v).

(9)

Proof: Select any q ∈ [q]. Note, from Eq. 8, that
distS D−1 (q, un ) ≥ distGr(1,D) ([q], [un ]) for all n andun ∈ [un ], thus, the function minimized in Eq. 9 will
never have a smaller minimum value than the one found in
Eq. 7. Note moreover that we can find un ∈ [un ] such that
distS D−1 (un , q) = distGr(1,D) ([un ], [q]) for all n. But thenq is a minimizer of Eq. 9 and hence a weighted average on
the sphere S D−1 .
Lemma 1 is useful because it reduces the problem ofcomputing averages on Gr(1, D) to a problem of computing
averages on S D−1 , optimized over the possible representatives ±u1:N ⊂ [u1:N ].
3.2

4

Computing the Average Subspace

To compute averages on Gr(1, D) we suggest a straightforward algorithm, which repeats the following two steps
until convergence

The algorithm is initialized with a uniform random sample
q0 . Each iteration of this algorithm has complexity O(N D),
which is the same as EM PCA (5). This algorithm was alsostudied by Kwak [12] who proved finite-time convergence.
3.2.1 Multiple Components
The above algorithm only computes the leading componentof the data. We compute additional components in a greedy
fashion by adding the further constraint that components
should be orthogonal. We compute this using deflation [16]by removing the estimated component from the data andcomputing the next component from the residual. Formally,
if X ∈ RN ×D denotes the data matrix, thenX̂ ← X − (Xq)q |(15)

is the updated data matrix where the component q isremoved. The next component can then be computed onX̂ as above. This approach is optimal for Gaussian data
(see Sec. 3.2.3). While we do not have results for other
distributions we show that it works well in practice. As

Copyright (c) 2016 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at http://dx.doi.org/10.1109/TPAMI.2015.2511743
SUBMITTED TO IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XXX, NO. XXX, MONTH YEAR

5

by the data. As PN (0, Σ) is mirror-symmetric we havep([un ]) = PN ({un , −un } | 0, Σ)

Fig. 3. Left: Zero-mean Gaussian data (orange) andthe unit circle (black). Right: The distribution of theGaussian data projected onto the unit circle. As thedata has zero mean, the projected distribution has two
equal-sized modes at opposing points. The Grassmann
average is such a mode.

(18)= PN (un | 0, Σ) + PN (−un | 0, Σ)(19)= 2PN (un | 0, Σ).

(20)

As PN (un | 0, Σ) = PN (−un | 0, Σ) we note that the
distribution of subspaces p([un ]) is unimodal with a mode
corresponding to the first principal component of the data;
see also Fig. 3.
Since p([un ]) is symmetric around the mode, it is naturalto think that the mode and the mean coincide. This, however,depends on the metric used to compare subspaces. The
argument holds for any rotationally invariant metric; in thefollowing we provide a proof for the metric in Eq. 16.

projection on a subspace does not remove outliers, thesame level of robustness is needed for every component. Theorem 2. Let x1:N ∼ N (0, Σ) be sampled from aDThis greedy approach is problematic for EM PCA as multivariate normal distribution on R . In the limit of
inaccuracy in one estimated component influences the infinite data, the average subspace (17) coincides with thefollowing components. This is not the case for the GA first principal component.
PN
algorithm as the found local minimum is exact to numerical
Proof: First note that N1 n=1 |x|n v| is the fiprecision due to the closed-form spherical average.
nite sample Rapproximation of the expectation value
E(|x| v|) = |x| v|p(x)dx. Therefore, in the limit of
PN
3.2.2 An Energy Minimization View
infinite data, arg maxv∈S D−1 N1 n=1 |x|n v| coincides with|The GA algorithm represented by (14) is derived as an arg max E(|x v|). Since x| is sampled from a normal
x v follows a univariate normal
algorithm for computing the average subspace spanned by distribution, the projection
2
distribution,
N
(0,
σ
),
with
σv2 = v | Σv [35, §2.3]. Thus,v
the data, but it proves insightful to look at which underlying|
energy is being optimized. From (8) and (10) we see that |x v| follows a half-normal distribution with expectedvalue proportional to σv [36]. The standard deviation, σv ,|2v is the principal eigenvector of Σ,
distGr(1,D) ([u1 ], [u2 ]) = 1 − |u1 u2 |,(16) is maximized whenthus arg max E(|x| v|) coincides with the first principal
component as defined by ordinary PCA.
and the average (7) reduces to[q] = arg maxv∈S D−1

NX|x|n v|.

(17)

3.2.4

A Probabilistic Interpretation

n=1

The idea of computing an average subspace implicitlyThis is an energy in the projection pursuit family [33]. assumes a unimodal subspace distribution p([un ]). As
beCompared with PCA, which maximizes explained variance [un ] = {un , −un } this implicit assumption can also
D−1
seen
as
a
mirror-symmetricbimodaldistributiononS;(3), we see that the GA can be expected to be more resilient
1/2p([u ]).
i.e.
p(u
)
=
p(−u
)
=
n
n
n
to outliers than ordinary PCA as the square has been replaced
An equivalent formulation of the Grassmann average can,with an absolute value. However, as the GA relies on an
ordinary average (12) it still has a break-down point of 0%. thus, be attained as a two-class clustering problem on the unit
In Sec. 3.3 we will improve upon this with robust averages. sphere with antipodal cluster centers. For clusters following
a von Mises-Fisher distribution [32], the Grassmann averagealgorithm (14) can be derived as a spherical two-class k3.2.3 Relationship to PCA (Gaussian Data)
means with antipodal means. This view has the advantage
To show how GA is related to ordinary PCA, we now that it provides a generative model for subspaces. Let theshow that they coincide when the observed data follows spherical clustering have cluster means ±µ, i.e.
a normal distribution. Before stating a formal result, we
1
1
provide some intuitions. Let x1:N ∼ N (0, Σ) be sampledp(un ) = vMF(un | µ, κ) + vMF(un | −µ, κ) (21)D
22from a multivariate zero-mean normal distribution on R ,Dκ /2−1
with a covariance with distinct eigenvalues. The subspaceexp (κµ| un ) , (22)vMF(un | µ, κ) =spanned by data point xn is represented by the natural(2π)D/2 ID/2−1 (κ)
xnbasis vector un = kx
. This follows a projected normal
nk
distribution [34], un ∼ PN (0, Σ). Our interest is in the where vMF is the von Mises-Fisher density; Iv the modified
distribution of the one-dimensional subspaces [un ] spanned Bessel function of the first kind and order v; and κ ≥ 0 a
Copyright (c) 2016 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at http://dx.doi.org/10.1109/TPAMI.2015.2511743
SUBMITTED TO IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XXX, NO. XXX, MONTH YEAR

Fig. 4. Gaussian data (orange) along with the averagesubspace (black line). The light gray subspace interval
corresponds to the 95% confidence interval of the mean
estimate under the probabilistic interpretation (27).
concentration parameter. The subspace distribution is thenp([un ]) = p(un ) + p(−un )1
1= vMF(un | µ, κ) + vMF(−un | µ, κ)2
2
1
1+ vMF(un | −µ, κ) + vMF(−un | −µ, κ)22= vMF(un | µ, κ) + vMF(−un | µ, κ)= vMF([un ] | µ, κ).

(23)
(24)
(25)
(26)
(27)

This model immediately allows us to sample new subspaces[32], estimate confidence intervals for the subspace [37,§5.3.1–5.3.2], and derive central limit theorems for the mean
estimate [38]. As an example, Fig. 4 shows 95% confidence
intervals for the estimation of the leading component fromGaussian data. This type of subspace reasoning is notavailable with conventional PCA models, but follows easily
from the subspace estimation formulation.

6Here, µrob can be chosen to have the robustness properties
relevant to the application. While we do not have a formal
general convergence proof, empirical experience indicatesthat the algorithm converges in finite time.
In computer vision applications, we often deal with image
data where individual pixels are corrupted. Consequently,
when computing average images robustly it is common
to do so on a per pixel basis; i.e. as a series of onedimensional robust averages. A standard approach to createa robust one-dimensional average is the trimmed average[2]. This removes the largest and smallest observations prior
to computing the mean. If we remove P % of the data fromboth sides, we get a break-down point of P %. For P = 50%
we get the maximally robust median [2].
To build a subspace estimator that is robust with respect topixel outliers, we pick µrob as the pixel-wise trimmed mean.
We call the resulting estimator the Trimmed Grassmann
Average (TGA) and let TGA(P, K) denote the estimator thatfinds K components and trims P percent of the smallest andlargest elements in the data. Note that P and K are the only
parameters of the estimator, where the choice P = 50%
gives the maximally robust median [2].
Each iteration of TGA(P, K) has computational complexity O(KN D) as the trimmed average can be computed inlinear time using a selection algorithm [39, §9]. We thus
have a robust subspace estimator that can deal with pixelwise outliers and has the same computational complexity as
both GA (14) and EM PCA (5). In the appendix we showthat a version of TGA(50%, 1) converges in finite time toa local minimum of the trimmed version of Eq. 7. We do
not have similar results for other versions of RGA, but
empirically we find that an optima is always found.
3.4

Practical Considerations

It has great practical value that a robust subspace estimatorcan be created using only a robust average, as it is often
easy to compute a robust average which takes the structure
3.3 Robust Grassmann Averagesof the problem into account. With this, the data can be madeThe core computational part of the GA algorithm (14) is“zero mean”, and a robust Grassmann average (RGA) can
the spherical average (11) which, in turn, merely computes
be computed by repeated averaging and deflation.
a normalized Euclidean average (12). From this we see that
Picking a suitable dimensionality for the estimatedeven a single outlier can arbitrarily corrupt the result.
subspace is a general problem that influences all current
A straight-forward solution to this issue is to replace themethods either directly or through a regularization parameter.
average with a robust average, giving the following scheme:
For PCA it is common to pick enough dimensions to explain
Algorithm: Robust Grassmann Average (RGA)e.g. 90% of the variance in the data. This approach is,
Let µrob denote any robust averaging operator.
however, troublesome in the presence of outliers as they willxn
Let un = kx
and
i
=
1.
dominate the variance estimate. Our probabilistic Grassmannnkrepeataverage (Sec. 3.2.4) provides an estimate of the confidence
interval of a given subspace. One option is to keep increasing|
αn ← sign(un qi−1 ) ∀n,
the dimensionality of the estimated subspace until thisµrob (w1:N , (αu)1:N )qi ←,(28) confidence interval is larger than a specified threshold. In
kµrob (w1:N , (αu)1:N )kpractice this is, however, also sensitive to outliers. In ouri←i+1experiments, we either rely on prior knowledge or use atrimmed variance estimator to determine dimensionality.
until kqi − qi−1 k <  or i > max_iter.
The RGA algorithm finds a local optimum, so the average
The algorithm is initialized with a uniform random sample subspaces are not guaranteed to be computed in order of
q0 . We call this the Robust Grassmann Average (RGA). decreasing variance. This is a general issue for methods
Copyright (c) 2016 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at http://dx.doi.org/10.1109/TPAMI.2015.2511743
SUBMITTED TO IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XXX, NO. XXX, MONTH YEAR

Groundhog Day

7

Pulp Fictionthat estimate one basis vector at a time; e.g. both EM PCA(85 frames)(85 frames)
and the method from De la Torre and Black [6] are not
TGA(50%, 80)
0.01570.04040.01680.0443Inexact ALM [5]guaranteed to provide ordered components. The standard
De la Torre and Black [6]0.03490.0599solution is a simple post hoc reordering according to theGA (80 comp)0.35510.3773amount of variance captured by each subspace.
0.35930.3789PCA (80 comp)
The recursive deflation approach (15) for estimating mulTABLE 1
tiple components relies on the specification of a projection The mean absolute reconstruction error of the different
operator. As linear projection is known to be sensitive toalgorithms on recent movies with added noise
outliers [40] it may be beneficial to use a robust projection
estimated from Nosferatu.
operator. We investigate this option empirically in Sec. 4.9.

4

0.4

R ESULTS

0.35

Groundhog Day
Pulp Fiction

Mean absolute deviation

In this section we compare our method with the approaches
0.3
of Candes et al. [5] and De la Torre and Black [6]. For
0.25
comparative purposes, all algorithms are implemented in
0.2
Matlab, and, unless otherwise stated, the experiments are
0.15
performed on an 6 core 3.3GHz Intel Xeon W3680 with 24
0.1
GB of memory. When using TGA we subtract the pixel-wise
0.05
median to get a robust estimation of a “zero mean” dataset,0
and we trim at 50% when computing subspaces. For [5] we
0
10
20
30
40
50
Trimming percentage (%)
use the default parameters of the published code2 . For [6]we also use the default parameters of the published code3 , Fig. 6. The mean absolute deviation from ground-truth
but with the same number of components as TGA.
measured at noisy pixels for two sequences.
Section 4.1 explores the application of robust PCA to
the restoration of archival film data, while in Sec. 4.2 and4.3 we repeat the experiments from [5]. We quantitatively material. TGA removes or downweights the outliers while
compare the robustness of different methods in Sec. 4.4 and mostly preserving the inliers. Both algorithms from Candes
investigate the empirical impact of trimming in Sec. 4.5. et al. [5] and De la Torre and Black [6] oversmooth the
Theorem 2 is empirically verified in Sec. 4.6, and the impact frames and drastically reduce the visual quality of the video.
of outliers is studied in Sec. 4.7. The random initialization This is visually evident in the supplementary video.
is studied in Sec. 4.8, different projections are compared in
Note that this is not a fully general algorithm for film
Sec. 4.9, and finally Sec. 4.10 measures scalability.
restoration but with modifications for camera stabilizationand dealing with fast-moving regions, TGA might form acomponent in a film restoration system, cf. [41].
4.1 Video Restoration
To quantify the performance of the algorithms, we
Archival films often contain scratches, gate hairs, splices,
generate artificially corrupted videos by taking the detected
dust, tears, brightness flicker and other forms of degradation.
outliers from Nosferatu, binarizing them to create an outlier
Restoring films often involves manual labor to identify andmask, translating these masks and combining them inrepair such damage. As these degradations vary in bothrandom order to produce new outlier masks and thenspace and time, we treat them collectively as pixel-level
applying them to clean frames from Hollywood movies.
outliers. Assuming a static camera, we perform robust PCAThis corrupts approximately 10% of the pixels. To measureon the frames of a movie clip and reconstruct the clip using
the accuracy of the different reconstructions, we measure
the robust subspace. We perform the reconstruction using
the mean absolute deviation between the original and theorthogonal projection; this choice is investigated in Sec. 4.9.
reconstructed images at the pixels where noise was added.
As an example, we consider a scene from the 1922 classic
This is shown in Table 1 where TGA is consistently better
movie Nosferatu. We restore the frames using TGA, [5],than the other algorithms, though Inexact ALM does almostand [6]. We expect the true rank of the inliers to correspondequally well. For completeness, Table 1 also shows theto the number of frames in the sequence. This consists ofaccuracy attained by ordinary PCA and GA (14).
80 frames, so for TGA and [6] we compute 80 components.
These results are obtained with a pixel-wise median; i.e.
This may intuitively seem like “too many” components,
50% trimming. Figure 6 shows the impact of this parameter.
but as robust methods disregard parts of the data (e.g.
When we trim more than 10–15% we see substantial
through trimming) they often allow for more components
improvements compared to no trimming.
than observations. Figure 5 shows two representative frames
and the restored movies are available in the supplementary
4.2 Background Modeling
2. http://perception.csl.illinois.edu/matrix-rank/sample code.html#
RPCA
3. http://users.salleurl.edu/∼ftorre/papers/rpca2.html

We now repeat an experiment from Candes et al. [5] on
background modeling under varying illumination conditions.

Copyright (c) 2016 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at http://dx.doi.org/10.1109/TPAMI.2015.2511743
SUBMITTED TO IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XXX, NO. XXX, MONTH YEAR

Original Data

TGA(50%, 80)

Inexact ALM [5]De la Torre and Black [6]8Fig. 5. Two representative frames from the 1922 film Nosferatu as well as their restoration using TGA and thealgorithms from [5], [6]. TGA removes many outliers and generally improves the visual quality of the film, while
Inexact ALM oversmoothes the results and De la Torre and Black oversmooth and introduce artifacts. This isalso seen in the absolute difference between the data and the reconstruction, which is also shown (inverted andmultiplied by 2 for contrast purposes). We also refer the reader to the supplementary video.
A static camera records a scene with walking people, while illumination conditions, and, hence, with different castthe lighting conditions change throughout the sequence. shadows [43]. Assuming the faces are convex Lambertian
Robust PCA should be able to capture the illumination objects, they should lie near a nine-dimensional linear spacedifference such that the background can be subtracted.
[44]. In this model, the cast shadows are outliers that can be
The first sequence, recorded at an airport [42], was used removed with robust PCA. We estimate images with smoothby Candes et al. [5] to compare with the algorithm from De shading, but without cast shadows using TGA(50%, 9).
We study the same people as Candes et al. [5]. Figure 8
la Torre and Black [6]. We only compare with [5], whichshows the original data, the reconstruction using different
was reported to give the best results.
Figure 7 (left) shows select frames from the airport algorithms as well as the absolute difference between images
sequence as well as the reconstruction of the frame with and reconstructions. While there is no “ground truth”, bothInexact ALM and TGA with 50% trimming and 5 compo- TGA and Inexact ALM [5] remove cast shadows equallynents. See also the supplementary video. In general, TGA well, though TGA keeps more of the shading variation andspecularity while Inexact ALM produces a more matte result.
gives noticeably fewer ghosting artifacts than [5].
We repeat this experiment on another video sequence Note that the theory says that a nine-dimensional subspace isfrom the CAVIAR4 dataset. The results, shown in Fig. 7 sufficient to capture most of the variance in the reflectance(right), confirm the results from the airport sequence as of a smooth Lambertian surface under distant isotropic
illumination [44]. It does not say that it has to contain only
TGA(50%, 5) produces noticeably less ghosting than [5].
Lambertian reflectance, particularly if the reflectance itself
is low dimensional.
4.3 Shadow Removal for Face Data
We repeat the shadow removal experiment from Candes etal. [5], which considers images of faces under different
4. http://homepages.inf.ed.ac.uk/rbf/CAVIAR/

4.4 Vector-level Robustness
So far, we have considered examples with pixel-wise outliers.
For completeness, we also consider a case where the entire

Copyright (c) 2016 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at http://dx.doi.org/10.1109/TPAMI.2015.2511743
SUBMITTED TO IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XXX, NO. XXX, MONTH YEAR

Data

TGA(50%, 5)Inexact ALM [5]

Data

TGA(50%, 5)9Inexact ALM [5]Fig. 7. Reconstruction of the background from video sequences with changing illumination; three representative
frames. Left: the airport sequence (3584 frames), which was also used in [5]. Right: a sequence from the CAVIAR
dataset (1500 frames).
Data

TGA(50%, 9)Inexact ALM

Data

TGA(50%, 9)Inexact ALM

Fig. 8. Shadow removal using robust PCA. We show the original image, the robust reconstructions as well as theirabsolute difference (inverted). TGA preserves more specularity, while Inexact ALM produce more matte results.
vector observation is corrupted. We take two 425-frame clips
from contemporary Hollywood movies (the same as in thequantitative experiment in Sec. 4.1) and treat one clip as theinliers. We then add an increasing number of frames from thesecond clip and measure how this influences the estimatedcomponents. We use the standard error measure, expressedvariance [3], which is the ratio between the amount ofvariance captured by the estimated principal component andthat captured by the ground truth component vgt ,PN(x| v)2
E(v) = PNn=1 |n.
(29)
2n=1 (xn vgt )level outliers and is a good basic alternative to PCA. Whenevery vector has some data-level outliers, GA performs like
PCA (Table 1), and TGA offers significant robustness.
4.5

Statistical Efficiency vs. Robustness

In Fig. 9a it is worth noting that TGA(25%, 1) performs
better than TGA(50%, 1) for a small number of outliers.
Medians are known [2] to require more observations togive accurate answers than a non-robust average, but on theother hand the median is maximally robust. The trimming
parameter, thus, provides a trade-off between statistical
efficiency and robustness. This is well-known for robust
Here we only consider one component to emphasize the estimators and it is not surprising that the TGA is affecteddifference between the methods; for [5] we remove outliers by this as well. Thus, one should use as low a trimmingand then reduce to a one-dimensional subspace using PCA. parameter as the data allows to get as efficient an estimator
Figure 9a shows the results. Both [5] and [6] have difficulties as possible; see e.g. [2] for further discussions on the matter.
with this problem, but this is not surprising as they areSome insight into the statistical efficiency of the TGA
designed for pixel-wise outliers. Both TGA and GA handle estimator can be gained by studying its behavior on data
the data more gracefully, in particular versus PCA.
with no outliers. Figure 9b shows the expressed variance as
In summary GA is quite robust when there are vector- a function of the number of observations N and the level of
Copyright (c) 2016 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at http://dx.doi.org/10.1109/TPAMI.2015.2511743
SUBMITTED TO IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XXX, NO. XXX, MONTH YEAR

a

10

1.1

b

c

Expressed Variance

Expressed Variance

Expressed Variance

1

0.9

0.8

0.7

0.6

Outlier Percentage (%)

0.50Number of observations (N)

d

1

0

0.96
0.94

25

0.92

0.9
0.9
0.88
50
0

0.2

0.4

0.6
rank / D

0.8

800

10

1

1000

0.96

20
0.95

30
40

0.94

50
0.93

60
70

37.5

Outlier Percentage (%)

400
600
Number of observations (N)fInitializations

Outlier Percentage (%)

Expressed Variance

0.95

200

e

0.98
12.5

EM PCA (5 comp)GA (5 comp)
TGA(50, 5)

0.92

80

0.91

90
100

0

12.5

25
37.5
Outlier Percentage (%)

50

0.9

Fig. 9. (a) The expressed variance for different methods as a function of the percentage of vector-level outliers.
Note that EM PCA and Inexact ALM provide near-identical solutions. (b) The expressed variance as a function of
the number of observations N . Here the dimensionality is fixed at D = 100. The different curves correspond todifferent levels of trimming. (c) Expressed variance as a function of number of samples from a known Gaussian
distribution. The results are for D = 30 dimensional data, and are averaged over 5 experiments with randomly
generated covariance matrices. (d) Expressed variance for Gaussian inliers with an increasing number of outliers.
The results are for D = 30 dimensional data, and are averaged over 50 experiments with randomly generated
covariance matrices. (e) Expressed variance (color coded from 0.86 (blue) to 1 (yellow)) as a function of the relative
rank of the inliers verses the percentage of outliers. (f) Impact of initialization is studied in TGA by measuring
expressed variance (see color code) for different random initializations.
trimming. Here the dimensionality is fixed at D = 100. In worse than GA due to lower statistical efficiency as described
low-sample scenarios the amount of trimming has a large in Sec. 4.5. Overall GA and TGA provide results that are
impact on the quality of the estimator. As the number of very similar to PCA for Gaussian data.
observations increase, the impact of trimming decreases.
Throughout this paper, we consistently trim 50%; i.e.
4.7 Performance with Vector-level Outliersa pixel-wise median. In some experiments slightly better
results can be attained by trimming less, but we opted for We further measure the methods’ performance when the data
is Gaussian with vector-level outliers. We sample inliersa median to avoid parameter tuning.
from a Gaussian distribution with a random covariance Σ.
Outliers are sampled from a Gaussian distribution with a
4.6 Performance on Gaussian Datamean placed along the smallest eigenvector of Σ, giving
Theorem 2 states that for Gaussian data, the component a consistent bias away from the principal components of
estimated by GA coincides with the principal component the inliers. Figure 9d shows the expressed variance fordefined by PCA. The result is, however, only in expectation an increasing number of outliers for EM PCA, GA andand may not hold for finite data samples. We investigate TGA with 50% trimming. As expected, the performance of
this empirically by sampling data from a Gaussian distri- EM PCA quickly drops, while the performance of GA drops
bution with known covariance, such that the ground truth more slowly. TGA, on the other hand, is very robust andprincipal component can be estimated through an eigen- attains a strong performance until the number of outliers
decomposition. Figure 9c shows the expressed variance as exceeds that of the inliers.
a function of the number of samples for both EM PCA, GA
We repeat the experiment for TGA with inliers of varying
and TGA with 50% trimming. EM PCA has slightly better rank. We consider N = 1000 inliers in D = 100 dimensions
performance compared to GA, which, in turn, is slightly with a rank between 1 and D, and add up to N outliers forbetter than TGA. This is to be expected as EM PCA directly a total of 2N observations. We estimate a robust subspace
optimizes to get the principal component. TGA does slightly using TGA(50%, 5). Figure 9e shows the expressed variance
Copyright (c) 2016 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at http://dx.doi.org/10.1109/TPAMI.2015.2511743
SUBMITTED TO IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XXX, NO. XXX, MONTH YEAR

as a function of both rank and number of outliers. In contrast
to Fig. 9d, the expressed variance is now color coded withyellow being 1 and dark blue being approximately 0.86.
TGA is robust over a broad range of ranks and outlier
percentages as can be seen by the contours showing robust
solutions capturing 0.9 and 0.95 of expressed variance.
4.8

Quality of Local Optima

11

difference between robust and orthogonal projection is dueto the outliers having a maximal value of 255.
It is worth noting that outliers generally do not influence
the deflation estimator (15) of multiple components. Afterthe leading component has been subtracted from the data,
the outliers will generally still be present in the data, if theleading component was estimated robustly. As long as arobust subspace estimator is applied, this procedure can be
repeated until the outliers dominate the residual data.

TGA provides a locally optimal solution, which depends ona random initialization. To understand the quality of these 4.10 Scalability
optima, we sample a single set of inliers (N = 1000, D =100), and add an increasing number of outliers as in Sec. 4.7. An important feature of the Grassmann average is that it
Figure 9f shows the expressed variance as a function of the scales gracefully to large datasets. To show this, we reportoutlier percentage for 100 different initializations. While the running time of the different algorithms for increasing
small fluctuations appear, it is evident that almost equally numbers of observations. We use video data recorded witha static camera looking at a busy parking lot over a twogood solutions are found for all initializations.
day period. Each frame has a resolution of 320 × 240 and
we consider up to 38, 000 frames, requiring a total of 22
4.9 Robust vs. Orthogonal Projections
GB of memory. Figure 10(right) shows the running times.
In the experiments, we estimate robust subspaces with a Inexact ALM [5] and the algorithm from De la Torre andbasis Q ∈ RD×K , where K is the number of components Black [6] quickly become impractical due to their memory
and D is the dimensionality of the data space. Data is then use; both algorithms run out of memory with more thanprojected into the subspace using orthogonal projection,
6000 observations. TGA, on the other hand, scales roughly|
P = XQQ ,(30) linearly with the data. In this experiment, we compute 100components, which is sufficient to capture the variation inwhere X ∈ RN ×D is the data matrix. This projection finds the background. For comparison, we also show the running
the subspace point P closest to the original point X; i.e. a time of EM PCA [24]. TGA is only slightly slower thanleast-squares optimization problem where (30) is the solution EM PCA, which is generally acknowledged as being among[40]. As least-squares problems are sensitive to outliers, it the most practical algorithms for ordinary PCA on large
is reasonable to have a robust projection operator.
datasets. It may be surprising that a robust PCA algorithm
One such robust projection operator is derived by Black can achieve a running time that is comparable to ordinary
and Jepson [40] using the Geman-McClure error function PCA. While each iteration of EM PCA is computationally[15]. This projection is then optimized using gradient cheaper than for TGA, we find that the required number ofdescent.
iterations is often lower for TGA than for EM PCA, which
We have experimented with the robust operator, and find explains the comparable running times. For 38, 000 images,that it provides minor improvements. We first consider a TGA(50%, 100) used 2–134 iterations per component, withsynthetic experiment: we select an increasing number of an average of 50 iterations. EM PCA, on the other hand,random pixels from the Groundhog Day sequence and make used 3–729 iterations with an average of 203 iterations.
them white, thereby adding biased outliers. We estimate
In the above, all algorithms are implemented in Matlab
a subspace from the noise-free data, and reconstruct the and have seen the same level of code optimization. We also
noisy data by projecting it into the subspace. We consider have a C++ implementation of TGA that is significantly
both orthogonal and robust projection. We measure the faster than the Matlab implementation (Fig. 10 (right)).
mean absolute difference between the projections from noise- Both implementations are available online1 .
free and noisy data. Figure 10(left) shows that the robust
To further emphasize the scalability of TGA, we computeprojection generally is better than orthogonal projection.
the 30 leading components of the entire Star Wars IV movie.
For non-synthetic examples we observe less of a dif- This consists of 179,415 RGB images in full DVD resolution
ference. To illustrate this, we reconstruct a frame from (704 × 306). This requires 431 GB of memory using single
the Nosferatu experiment using both orthogonal and robust precision. Computing 30 components takes 34 hours on an
projection. The data and results are given in Fig. 10(center), AMD Opteron 6378 (64 cores) with 1TB of memory, runningwhere the different projections appear to give nearly identical at 2.4GHz, using the C++ implementation of TGA. It takes
results. Minor differences between the methods do appear 24 hours to compute the leading 30 principal components
and the figure shows the pixels where they disagree. When using a C++ implementation of EM PCA with the samethe images are represented using 8 bits per pixel (i.e. multi-threading scheme.
Figure 11 shows the leading odd-numbered components
between 0 and 255), the largest observed absolute difference
as computed by TGA with 50% trimming. We note somebetween the reconstructions was 1 gray level.
As the use of a robust projection only resulted in small similarities between the estimated components and those
improvements, we opted for the much faster orthogonal attained from the analysis of natural images [45]. An appliprojection throughout the paper. We speculate that the small cation made possible by the scalability of the Grassmann
Copyright (c) 2016 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at http://dx.doi.org/10.1109/TPAMI.2015.2511743
SUBMITTED TO IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XXX, NO. XXX, MONTH YEAR

0.18
0.16

12

Orthogonal projection
Robust projection

Mean absolute deviation

0.14
0.12
0.1

Original data

Orthogonal projection

Robust projection

Thresholded difference

0.08
0.06
0.04
0.02
0
0

10

20
30
Outlier percentage (%)

40

50Fig. 10. Left: Mean absolute deviation of a reconstruction attained by different projection operators. Robust
projection generally improves results. Center: Comparison between orthogonal and robust projection. The panel
shows the data, its reconstruction using both orthogonal and robust projections, as well as the thresholded
difference between the two projections methods; when pixels are represented using 8 bits (between 0 and 255),
the maximal observed absolute pixel-difference is 1. Right: Running time of the algorithms as a function of thesize of the problem. The dimensionality is fixed at D = 320 × 240 while the number of observations N increases.
TGA is comparable to EM PCA, while [5] and [6] are unable to handle more than 6000 observations.
average is the recent work on optical flow from Wulff &Black [46]. They learn a robust basis for full-frame opticalflow using TGA on 180,000 flow fields computed frommovies. Robustness is important because the computed flow
fields contain outliers.

5

D ISCUSSION

However, this is not a proper metric as it does not satisfy
the triangle inequality. It is unclear if a proper metric exists
under which standard PCA is the average subspace.
A perhaps more natural distance measure between subspaces is the angle between them; this is known as theintrinsic metric. Unfortunately, even on a sphere, thecomputation of an average under the intrinsic metric requires
non-linear optimization [47]. This may limit the practicaluse of the intrinsic metric.

The principle idea of this work is to rephrase linear subspaceestimation directly in the space of subspaces. In particular,
we suggest to compute the average subspace spanned by
the data, i.e. the first moment on the Grassmann manifold. 5.2 Probabilistic Interpretation
As standard PCA is concerned with the second moment On the Grassmann manifold Gr(1, D), it is not hard toof the data (covariance) our approach simplifies subspace develop distributions of the subspaces spanned by the data.
estimation. This is also evident in the efficient algorithm We considered a von Mises-Fisher distribution as its modeand the ease with which robust extensions are developed.
coincides with the Grassmann average. Other choices are
The concept of studying the subspaces spanned by the equally possible, though the choice of metric will influence
data (rather than the data itself) opens several new directions which distributions are easy to work with.
which can be explored; here we discuss a few.
Having a distribution over subspaces has several advantages not found in classic subspace estimators. It becomes
5.1 Choice of Metriceasy to specify a prior over which subspaces are of interest;
The general formulation of the average subspace (7) allows confidence intervals for the estimated subspace can befor any metric on the Grassmann manifold. We picked the derived, allowing hypothesis testing; new random subspaces
Euclidean metric (10) as it both allows for closed-form can be drawn from the distribution; etc.
spherical averages and makes it straightforward to extend to
pixel-wise outliers. While these advantages are substantial, 5.3 Online Algorithms
other metrics may be worth investigating.
We have not explored an online extension of the GA
A first question is: does standard PCA compute analgorithm, which would only see one data point at a time.
average subspace under a particular metric? Standard
As we are only considering averages, we expect such online
PCA minimizes the reconstruction error (1), which canalgorithms exist. As central limit theorems are available on
be rewritten as
the Grassmann manifold [38] it should be further possibleNX
to bound the performance of such online estimators.
2
2
qPCA = arg min
wn 1 − cos (θn (v)) ,(31)
v∈S D−1 n=1

where wn = kxn k and θn (v) denotes the angle between xn
and v. This implies that PCA computes an average subspace
under the metricdist2PCA (xn , v) = 1 − cos2 (θn (v)) .

(32)

5.4

Non-orthogonal Components?

In our formulation, we require orthogonal subspaces when
we compute multiple components. This has computational
advantages, but it may not always be ideal. If orthogonality
is not a desired property, one can, instead of computing

Copyright (c) 2016 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at http://dx.doi.org/10.1109/TPAMI.2015.2511743
SUBMITTED TO IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XXX, NO. XXX, MONTH YEAR

13

Fig. 11. The leading odd-numbered principal components of Star Wars IV as computed by TGA(50%, 30)). The
top row contains components 1, 3, 5, 7, and 9 from left to right, and so forth.
averages, seek modes of the distribution of subspaces.
Finding such modes can be done, for example, with k-meansor mean-shift [48] directly on the Grassmann manifold.
These are straightforward exercises.

6

C ONCLUSION

Principal component analysis is a fundamental tool for data
analysis and dimensionality reduction. Previous work hasaddressed robustness at both the data vector and vectorelement level but fails to scale to large datasets. This istroublesome for big data applications where the likelihoodof outliers increases as data acquisition is automated.
We introduce the Grassmann average (GA), which is a
simple and highly scalable approach to subspace estimationthat coincides with PCA for Gaussian data. We have further
shown how this approach can be made robust by using
a robust average, yielding the robust Grassmann average(RGA). For a given application, we only need a robust
average to produce a suitable robust subspace estimator. We
develop the trimmed Grassmann average (TGA), which isa robust subspace estimator working at the vector-elementlevel. This has the same computational complexity as thewell-known scalable EM PCA [24], and empirical results
show that TGA is not much slower while being substantiallymore robust. Our implementation is available online1 .
The availability of a scalable robust PCA algorithm opens
up to many new applications. We have shown that TGA
performs well on different tasks in computer vision, where
alternative algorithms either produce poor results or failto run at all. We can even compute robust components ofentire movies at full resolution in a reasonable time. Further,
the ubiquity of PCA makes RGA relevant beyond computer
vision; e.g. for large datasets found in biology, physics andweather forecasting. Finally, as our approach is based onstandard building blocks like trimmed averages, it is very
easy to speedup via parallelization.

ACKNOWLEDGMENTS
Søren Hauberg is funded in part by the Villum Foundationand the Danish Council for Independent Research (DFF),
Natural Sciences; Aasa Feragen is funded in part by the
DFF, Technology and Production Sciences. The authors are
grateful to anonymous reviewers for helpful feedback.

R EFERENCES
N. Campbell, “Robust procedures in multivariate analysis I: Robust
covariance estimation,” Applied Statistics, vol. 29, no. 3, 1980.
[2] P. J. Huber, Robust Statistics. Wiley: New York, 1981.
[3] J. Feng, H. Xu, and S. Yan, “Robust PCA in high-dimension: A
deterministic approach,” in International Conference on Machine
Learning, 2012.
[4] L. Xu and A. L. Yuille, “Robust principal component analysis by
self-organizing rules based on statistical physics approach,” IEEE
Transactions on Neural Networks, vol. 6, no. 1, pp. 131–143, 1995.
[5] E. J. Candes, X. Li, Y. Ma, and J. Wright, “Robust principal
component analysis?” Journal of the Association for Computing
Machinery, vol. 58, no. 3, 2011.
[6] F. De la Torre and M. J. Black, “A framework for robust subspacelearning,” International Journal of Computer Vision, vol. 54, pp.
117–142, 2003.
[7] S. Hauberg, A. Feragen, and M. J. Black, “Grassmann averages for
scalable robust PCA,” in Computer Vision and Pattern Recognition(CVPR), 2014.
[8] I. Jolliffe, Principal component analysis. Springer, 2002.
[9] F. H. Ruymgaart, “A robust principal component analysis,” Journalof Multivariate Analysis, vol. 11, pp. 485–497, 1981.
[10] H. Aanæs, R. Fisker, K. Åström, and J. Carstensen, “Robust
factorization,” IEEE Transactions on Pattern Analysis and MachineIntelligence, vol. 24, no. 9, pp. 1215–1225, Sep 2002.
[11] C. Ding, D. Zhou, X. He, and H. Zha, “R1 -PCA: rotational
invariant l1-norm principal component analysis for robust subspacefactorization,” in International Conference on Machine Learning,
2006, pp. 281–288.
[12] N. Kwak, “Principal component analysis based on l1-norm maximization,” IEEE Transactions on Pattern Analysis and MachineIntelligence, vol. 30, no. 9, pp. 1672–1680, Sept 2008.
[13] M. McCoy and J. A. Tropp, “Two proposals for robust pca using
semidefinite programming,” Electronic Journal of Statistics, vol. 5,
pp. 1123–1160, 2011.
[14] A. Naor, O. Regev, and T. Vidick, “Efficient rounding for thenoncommutative grothendieck inequality,” in Proceedings of the Fortyfifth Annual ACM Symposium on Theory of Computing, ser. STOC
’13. ACM, 2013, pp. 71–80.
[15] S. Geman and D. McClure, “Statistical methods for tomographic image reconstruction,” Bulletin of the International Statistical Institute,vol. 52, no. 4, pp. 5–21, 1987.
[16] G. H. Golub and C. F. Van Loan, Matrix Computations, 3rd ed. The
Johns Hopkins University Press, Oct. 1996.
[17] L. W. Z. Lin, M. Chen and Y. Ma, “The augmented lagrange multipliermethod for exact recovery of corrupted low-rank matrices,” UILUENG-09-2215, Tech. Rep., 2009.
[18] R. Tron and R. Vidal, “Distributed computer vision algorithms throughdistributed averaging,” in IEEE Conference on Computer Vision and
Pattern Recognition (CVPR), June 2011, pp. 57–63.
[19] P. Netrapalli, N. U N, S. Sanghavi, A. Anandkumar, and P. Jain, “Nonconvex robust pca,” in Advances in Neural Information Processing
Systems 27, Z. Ghahramani, M. Welling, C. Cortes, N. Lawrence, and
K. Weinberger, Eds. Curran Associates, Inc., 2014, pp. 1107–1115.
[1]

Copyright (c) 2016 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at http://dx.doi.org/10.1109/TPAMI.2015.2511743
SUBMITTED TO IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XXX, NO. XXX, MONTH YEAR[20] L. W. Mackey, A. S. Talwalkar, and M. I. Jordan, “Divide-and-conquer
matrix factorization,” in Advances in Neural Information Processing
Systems 24, J. Shawe-Taylor, R. Zemel, P. Bartlett, F. Pereira, and
K. Weinberger, Eds., 2011, pp. 1134–1142.
[21] R. Liu, Z. Lin, Z. Su, and J. Gao, “Linear time principal component
pursuit and its extensions using l1 filtering,” Neurocomputing, vol.
142, pp. 529 – 541, 2014, SI Computational Intelligence Techniques
for New Product Development.
[22] Y. Mu, J. Dong, X. Yuan, and S. Yan, “Accelerated low-rank visual
recovery by random projection,” in IEEE Conference on Computer
Vision and Pattern Recognition (CVPR), June 2011, pp. 2609–2616.
[23] R. Arora, A. Cotter, K. Livescu, and N. Srebro, “Stochastic optimization for PCA and PLS,” in Communication, Control, and Computing(Allerton). IEEE, 2012, pp. 861–868.
[24] S. T. Roweis, “EM algorithms for PCA and SPCA,” in Advances in
Neural Information Processing Systems (NIPS). MIT Press, 1998.
[25] Y. M. Lui, “Advances in matrix manifolds for computer vision,”
Image and Vision Computing, vol. 30, pp. 380–388, 2012.
[26] R. Vidal, Y. Ma, , and S. Sastry, “Generalized principal component
analysis,” IEEE Transactions on Pattern Analysis and MachineIntelligence, vol. 27, no. 12, pp. 1945–1959, 2005.
[27] P. Favaro, R. Vidal, and A. Ravichandran, “A closed form solution
to robust subspace estimation and clustering,” in Computer Vision
and Pattern Recognition (CVPR), 2011 IEEE Conference on. IEEE,
2011, pp. 1801–1807.
[28] E. Elhamifar and R. Vidal, “Sparse subspace clustering: Algorithm,
theory, and applications,” Pattern Analysis and Machine Intelligence,
IEEE Transactions on, vol. 35, no. 11, pp. 2765–2781, 2013.
[29] G. Liu, Z. Lin, and Y. Yu, “Robust subspace segmentation by low-rank
representation,” in Proceedings of the 27th International Conference
on Machine Learning (ICML-10), 2010, pp. 663–670.
[30] J. Lee, Introduction to Smooth Manifolds. Springer, 2002.
[31] M. R. Bridson and A. Haefliger, Metric spaces of non-positive
curvature. Springer, 1999.
[32] K. V. Mardia and P. E. Jupp, Directional Statistics. Wiley, 1999.
[33] P. J. Huber, “Projection pursuit,” Annals of Statistics, vol. 13, no. 2,
pp. 435–475, 1985.
[34] F. Wang, “Space and space-time modeling of directional data,” Ph.D.
dissertation, Duke University, 2013.
[35] C. M. Bishop, Pattern Recognition and Machine Learning. Springer,
2006.
[36] F. C. Leone, L. S. Nelson, and R. B. Nottingham, “The folded normal
distribution,” Technometrics, vol. 3, no. 4, pp. 543–550, 1961.
[37] N. I. Embleton, T. Fisher, and B. J. J. Lewis, Statistical analysis ofspherical data. Cambridge University Press, 1993.
[38] R. Bhattacharya and V. Patrangenaru, “Large sample theory of
intrinsic and extrinsic sample means on manifolds,” The Annalsof Statistics, vol. 31, no. 1, pp. 1–29, 2003.
[39] T. H. Cormen, C. E. Leiserson, R. L. Rivest, and C. Stein, Introduction
to Algorithms, 3rd ed. MIT Press, 2009.
[40] M. J. Black and A. D. Jepson, “Eigentracking: Robust matching andtracking of articulated objects using a view-based representation,”International Journal of Computer Vision, vol. 26, pp. 63–84, 1998.
[41] H. Ji, S. Huang, Z. Shen, and Y. H. Xu, “Robust video restoration
by joint sparse and low rank matrix approximation,” SIAM Journal
on Imaging Sciences (SIIMS), vol. 4, no. 4, pp. 1122–1142, 2011.
[42] L. Li, W. Huang, I. Gu, and Q. Tian, “Statistical modeling of complex
backgrounds for foreground object detection,” IEEE Transactions on
Image Processing, vol. 13, no. 11, pp. 1459–1472, 2004.
[43] A. Georghiades, P. Belhumeur, and D. Kriegman, “From few to many:
Illumination cone models for face recognition under variable lighting
and pose,” IEEE Transactions on Pattern Analysis and MachineIntelligence, vol. 23, no. 6, pp. 643–660, 2001.
[44] R. Basri and D. Jacobs, “Lambertian reflectance and linear subspaces,”
IEEE Transactions on Pattern Analysis and Machine Intelligence,
vol. 25, no. 2, pp. 218–233, 2003.
[45] P. J. Hancock, R. J. Baddeley, and L. S. Smith, “The principal
components of natural images,” Network: computation in neuralsystems, vol. 3, no. 1, pp. 61–70, 1992.
[46] J. Wulff and M. J. Black, “Efficient sparse-to-dense optical flow
estimation using a learned basis and layers,” in IEEE Conf. on
Computer Vision and Pattern Recognition (CVPR) 2015, Jun. 2015.
[47] X. Pennec, “Probabilities and statistics on Riemannian manifolds:
Basic tools for geometric measurements,” in Proceedings of Nonlinear
Signal and Image Processing, 1999, pp. 194–198.

14[48] H. Çetingul and R. Vidal, “Intrinsic mean shift for clustering on
Stiefel and Grassmann manifolds,” in IEEE Conference on Computer
Vision and Pattern Recognition (CVPR), June 2009, pp. 1896–1902.

Søren Hauberg received his Ph.D. in computer science from the Univ. of Copenhagenin 2011. He has been a visiting scholar at
UC Berkeley (2010), and a post doc at thePerceiving Systems department at the Max
Planck Institute for Intelligent Systems inTübingen, Germany (2012–2014). He is currently a post doc at the Section for Cognitive
Systems, at the technical Univ. of Denmark.
His research is at the interplay of geometry
and statistics.

Aasa Feragen holds a Ph.D. in mathematics
from the Univ. of Helsinki (2010) and hassince been post doc and later Associate
Professor at DIKU at the Univ. of Copenhagen.
From 2012 to 2014 she was a post doc in the
Machine Learning and Computational Biology
Research Group at the Max Planck institutefor Developmental Biology and Max Planck
Institute for Intelligent Systems, Tübingen,Germany. Her research is on geometric modeling for imaging and vision.

Raffi Enficiaud received his M.Sc. fromthe National Institute of Telecommunications,
Paris-Evry, France, in 2000, and his Ph.D.
from Paris School of Mines, France, in 2007.
He worked as a research engineer for DxO
Labs, Paris and INRIA Paris-Rocquencourt(2007-2013). He is currently running the Software Workshop at the Max Planck Institutefor Intelligent Systems, Tübingen, Germany.
His research interests include computer visionand image processing.

Michael J. Black received his B.Sc. fromthe Univ. of British Columbia (1985), his M.S.
from Stanford (1989), and his Ph.D. from Yale
(1992). After research at NASA Ames and theUniv. of Toronto, he was at Xerox PARC as amember of research staff and area manager(1993-2000). From 2000 to 2010 he was on
the faculty of Brown Univ. (Assoc. Prof. 20002004, Prof. 2004-2010). Since 2011 his oneof the founding director at the Max Planck
Institute for Intelligent Systems in Tübingen,Germany, where he leads the Perceiving Systems department.
Significant awards include the 2010 Koenderink Prize (ECCV) andthe 2013 Helmholtz Prize (ICCV). He is a foreign member of theRoyal Swedish Academy of Sciences and is a co-founder and board
member of Body Labs Inc.

Copyright (c) 2016 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.