IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XXX, NO. XXX, MONTH YEAR

1

Intrinsic Grassmann Averages for Online Linear,
Robust and Nonlinear Subspace Learning
Rudrasis Chakraborty, Liu Yang, Søren Hauberg and Baba C. Vemuri∗ , Fellow, IEEE
Abstract—Principal Component Analysis (PCA) and Kernel Principal Component Analysis (KPCA) are fundamental methods in machinelearning for dimensionality reduction. The former is a technique for finding this approximation in finite dimensions and the latter is often in
an infinite dimensional Reproducing Kernel Hilbert-space (RKHS). In this paper, we present a geometric framework for computing theprincipal linear subspaces in both (finite and infinite) situations as well as for the robust PCA case, that amounts to computing the
intrinsic average on the space of all subspaces: the Grassmann manifold. Points on this manifold are defined as the subspaces spanned
by K -tuples of observations. The intrinsic Grassmann average of these subspaces are shown to coincide with the principal componentsof the observations when they are drawn from a Gaussian distribution. We show similar results in the RKHS case and provide an efficientalgorithm for computing the projection onto the this average subspace. The result is a method akin to KPCA which is substantially faster.
Further, we present a novel online version of the KPCA using our geometric framework. Competitive performance of all our algorithms aredemonstrated on a variety of real and synthetic data sets.
Index Terms—Online Subspace Learning, Robust, PCA, Kernel PCA, Grassmann Manifold, Fréchet Mean, Fréchet Median.

F

1

I NTRODUCTIONP

RINCIPAL component analysis (PCA), a key work-horseof machine learning, can be derived in many ways:
Pearson [1] proposed to find the subspace that minimizes the
projection error of the observed data; Hotelling [2] insteadsought the subspace in which the projected data has maximal
variance; and Tipping & Bishop [3] considered a probabilistic
formulation where the covariance of normally distributed
data is predominantly given by a low-rank matrix. All thesederivations led to the same algorithm. There are many
generalizations of PCA including but not limited to sparse
PCA [4], generalized PCA [5] and we refer the interested
reader to a recent survey article [6]. We will not be addressing
any of these in our current work. Recently, Hauberg et al.
[7] noted that the average of all one-dimensional subspaces
spanned by normally distributed data coincides with theleading principal component. They computed the averageover the Grassmann manifold of one-dimensional subspaces(cf. Sec. 2). This average was computed very efficiently, but
unfortunately their formulation does not generalize to higherdimensional subspaces.
In this paper, we provide a formulation for estimating the
average K -dimensional subspace spanned by the observed
data, and present a very simple online algorithm for computing this average. This algorithm is parameter free, andhighly reliable, which is in contrast to current state-of-theart online techniques, that are sensitive to hyper-parametertuning. When the data is normally distributed, we show that
our estimated average subspace coincides with that spanned•••

R. Chakraborty and B. C. Vemuri are with University of Florida,
Gainesville, FL, USA. ∗ indicates corresponding author
E-mail: rudrasischa@gmail.com, vemuri@ufl.edu
S. Hauberg is at the Technical University of Denmark, Kgs. Lyngby,
Denmark.
E-mail: sohau@dtu.dk
L. Yang is at the University of California, Berkeley.
E-mail: liu-yang@berkeley.edu

by the leading K principal components. Moreover, since ouralgorithm is online, it has a linear complexity in terms ofthe number of samples. Furthermore, we propose an online
robust subspace averaging algorithm which can be used to
get the leading K robust principal components. Analogous
to its non-robust counterpart, it has a linear time complexity
in terms of the number of samples. In this article, we extend
an early conference paper [8] and perform online non-linear
subspace learning which is an online version of the kernel
PCA. In comparison to our preliminary work in [8], thispaper contains a more detailed analysis in addition to thegeneralization akin to kernel PCA [9].
1.1

Related WorkIn this paper we consider a simple linear dimensionality
reduction algorithm that works in an online setting, i.e.
only uses each data point once. There are several existing
approaches in the literature that tackle the online PCA andthe online Robust PCA problems and we discuss some ofthese approaches here before discussing some related works
on kernel PCA:
Online PCA and online robust PCA: Oja’s rule [10] is aclassic online estimator for the leading principal components
of a dataset. Given a basis Vt−1 ∈ RD×K this is updated
recursively via Vt = Vt−1 + γt Xt (XtT Vt−1 ) upon receiving
the observation Xt . Here γt is the step-size (learning rate)
parameter that must be set manually; small values yields
slow-but-sure convergence, while larger values may leadto fast-but-unstable convergence. Several variants of Oja’s
method have been proposed in literature and analyzed
theoretically. Most recently, Allen-Zhu and Li [11] presented
global convergence analysis of Oja’s algorithm. Their algorithm convergence is however dependent on the gap between
k th and (k + 1)th eigen values (which in practice is unknown
without already knowing the eigen spectrum) and a gap freeIEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XXX, NO. XXX, MONTH YEAR

result which further depends on another parameter called‘virtual gap’. In contrast, our proposed online algorithm forcomputing the principal components is parameter free.
EM PCA [12] is usually derived for probabilistic PCA,
and is easily be adapted to the online setting [13]. Here, theE- and M-steps are given by:T
T(E-step) Yt =(Vt−1
Vt−1 )−1 (Vt−1
Xt )

(1)

(M-step) Ṽt =(Xt YtT )(Yt YtT )−1 .

(2)

The basis is updated recursively via the recursion, Vt =(1 − γt )Vt−1 + γt Ṽt , where γt is a step-size.
GROUSE and GRASTA [14], [15] are online PCA andmatrix completion algorithms. GRASTA can be applied toestimate principal subspaces incrementally on subsampled
data. Both of these methods are online and use rank-one
updating of the principal subspace at each iteration. GRASTA
is an online robust subspace tracking algorithm and canbe applied to subsampled data and specifically matrix
completion problems. He et al. [15] proposed an `1 -norm
based fidelity term that measures the error between thesubspace estimate and the outlier corrupted observations.
The robustness of GRASTA is attributed to this `1 -norm
based cost. Their formulation of the subspace estimation
involves the minimization of a non-convex function in an
augmented Lagrangian framework. This optimization iscarried out in an alternating fashion using the well-known
ADMM [16] for estimating a set of parameters involving
the weights, the sparse outlier vector and the dual vector in
the augmented Lagrangian framework. For fixed estimatedvalues of these parameters, they employed an incremental
gradient descent to solve for the low dimensional subspace.
Note that the solution obtained is not the optimum of the
combined non-convex function of GRASTA.
Recursive covariance estimation [17] is straight-forward, andthe principal components can be extracted via standard eigen
decompositions. Boutsidis et al. [17] considered efficientvariants of this idea, and provided elegant performancebounds. The approach does not however scale to highdimensional data as the covariance cannot practically be
stored in memory for situations involving very large data
sets as those considered in our work.
Candes et al. [18] formulated Robust PCA (RPCA) as
separating a matrix into a low rank (L) and a sparse matrix(S ), i.e., data matrix X ≈ L + S . They proposed Principal
Component Pursuit (PCP) method to robustly find theprincipal subspace by decomposing into L and S . They
showed that both L and S can be computed by optimizing an objective function which is a linear combinationof nuclear norm on L and `1 -norm on S . Recently, Loiset al. [19] proposed an online RPCA algorithm to solve
two interrelated problems, matrix completion and online
robust subspace estimation. Candes et al. [18] had someassumptions including a good estimate of the initial subspace
and that the basis of the subspace is dense. Though theauthors have shown correctness of their algorithm underthese assumptions, they are often not practical. In [20],authors proposed an alternative approach to solve RPCA
by developing an augmented Lagrangian based approach,
IALM.
A practical application to do robust PCA is where oneneeds to deal with missing data. Though incremental SVD

2(ISVD) [21] can be used to develop online PCA, the extension
to missing data is not straightforward. In literature, there areseveral methods to do online robust PCA for missing datacases including MDISVD [22], Brand’s [23] and PIMC [24].
All these approaches start with learning the projection matrix
to the coordinates which are present for a data sample. Thesemethods differ in the way they compute singular value
matrix. If Sn is the nth singular value matrix, then thesemethods try to find λS where λ = 1 for a natural extensionof ISVD for missing data, denoted by MDISVD [22]. Brand[23] used λ ∈ (0, 1) while in PIMC [24], the authors chose
λ based on the norm of the data observed so far. PETRELS[25] is an extension of PAST [26] on missing data. PAST usedsubspace tracking idea with an approximation of projection
operator. Several other researchers have proposed subspacetracking algorithm for missing data including REPROCS[27].
In another recent work, Ha and Barber [28] proposed
an online RPCA algorithm when X = (L + S)C where C
is a data compression matrix. They proposed an algorithm
to extract L and S when the data X are compress-sensed.
This problem is quite interesting in its own right but not
something pursued in our work presented here. Feng et al.
[29] solved RPCA using a stochastic optimization approach.
They showed that if each observation is bounded, then theirsolution converges to the batch mode RPCA solution, i.e.,their sequence of robust subspaces converges to the “true”
subspace. Hence, they claimed that as the “true subspace”(subspace recovered by RPCA) is robust, so is their online
estimate. Though their algorithm is online, the optimization
steps ( Algorithm 1 in [29]) are computationally expensive forhigh-dimensional data. In an earlier paper, Feng et al. [30]proposed a deterministic approach to solve RPCA (dubbed
DHR-PCA) for high-dimensional data. They also showedthat they can achieve maximal robustness, i.e., a breakdown
point of 50%. They proposed a robust computation of thevariance matrix and then performed PCA on this matrix
to get robust PCs. This algorithm is suitable for very highdimensional data. As most of our real applications in thispaper are in very high dimensions, we find DHR-PCA to
be well suited to carry out comparisons with. For further
literature study on this rich topic, we refer the reader to [31],[32].

Fig. 1. The average of two subspaces.

Online kernel PCA: In this paper, we propose an online
subspace averaging algorithm that we also extend to analgorithm akin to kernel PCA [9], i.e., to compute nonlinear subspaces. We show that with the popular kernel

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XXX, NO. XXX, MONTH YEAR

3

trick, we can extend our subspace averaging algorithm to is denoted by Gr(K, D), where K ∈ Z+ , D ∈ Z+ , D ≥ K .
compute a non-linear average subspace. But, because of A special case of the Grassmannian is when K = 1, i.e.,
the infinite dimensionality of the reproducing kernel Hilbert the space of one-dimensional subspaces of RD , which isspace (RKHS), it is not computationally feasible to make known as the real projective space (denoted by RP D ). A point
this an online algorithm. We however resolve this problem X ∈ Gr(K, D) can be specified by a basis, X , i.e., a set of K
by a finite approximation of the kernel using the method linearly independent vectors in RD (the columns of X ) that
proposed in [33] leading to an online non-linear subspace spans X . We have X = Col(X) if X is a basis of X , where
averaging algorithm. The past work in this context includes Col(·) is the column span operator. Given X , Y ∈ Gr(K, D),extension of Oja’s rule to perform kernel PCA [9]. Honeine with their respective orthonormal basis X and Y , the unique[34] proposed an online kernel PCA algorithm. They pointed geodesic ΓY
X : [0, 1] → Gr(K, D) between X and Y is givenout that as the principal vector is a linear combination of by:
the kernel functions associated with the available training
ΓY(3)
X (t) = Col X V̂ cos(Θt) + Û sin(Θt)
data, it becomes a bottleneck in making the kernel PCAY
T
online. They overcame this problem by controlling the order with ΓY=
X (0) = X and ΓX (1) = Y , where, Û Σ̂V̂TT
−1of the model. The algorithm starts with a set of pre-selected (I − XX )Y (X Y ) is the “thin” singular value decomkernel functions. Upon the arrival of a new observation, position (SVD), (i.e., Û is D × K and V̂ is K × K column
the algorithm decides whether to include or discard the orthonormal matrix, and Σ̂ is K × K diagonal matrix), andobservation from the set of kernel functions. Thus, by Θ = arctan Σ̂. The length of the geodesic constitutes the
restricting the number of kernel functions, they made their geodesic distance on Gr(K, D), d : Gr(K, D) × Gr(K, D)
kernel PCA algorithm an online algorithm. Ghashami et al. → R+ ∪ {0} which is as follows: Given X , Y with respective[35] proposed a streaming kernel PCA algorithm (SKPCA) orthonormal basis X and Y ,which uses the finite approximation of kernel and stores av
uK
small set of basis in an online fashion. Some other online
uX
t
d(X
,
Y)
,
arccos(σi )2 ,(4)
KPCA algorithms include [36], [37].
i=1Motivation and contribution: Our work is motivatedby the work presented by Hauberg et al. [7], who showed where Ū ΣV̄ T = X T Y is the SVD of X T Y , and,that for a data set drawn from a zero-mean multivariate [σ , . . . , σ ] = diag(Σ). Here arccos(σ ) is known as the1Ki
Gaussian distribution, the average subspace spanned by ith principal angle between subspace X and Y . Hence, thethe data coincides with the leading principal component. geodesic distance is defined as the ` -norm of principal
2
NThis idea is sketched in Fig. 1. Given, {xi }i=1 ⊂ RD , the angles.
1-dimensional subspace spanned by each xi is a point on the
Grassmann manifold (Sec. 2). Hauberg et al. then computed
the average of these subspaces on the Grassmannian using 2.2 The Intrinsic Grassmann Average (IGA)
1an “extrinsic” metric, i.e. the Euclidean distance. Besides We now consider intrinsic averages (IGA) on the Grassmannian.
Toexaminetheexistence
anduniqueness of IGA, we
the theoretical insight, this formulation gave rise to highly
needtodefine
an
open
ball
on
the
Grassmannian.
efficient algorithms. Unfortunately, the extrinsic approach
is limited to one-dimensional subspaces, and Hauberg et al. Definition 1 (Open ball). An open ball B(X , r) ⊂ Gr(K, D)
resorted to deflation methods to estimate higher dimensional of radius r > 0 centered at X ∈ Gr(K, D) is defined as
subspaces. We overcome this limitation with an intrinsic
B(X , r) = {Y ∈ Gr(K, D)|d(X , Y) < r} .
(5)metric, extend the theoretical analysis of Hauberg et al.,and provide an efficient online algorithm for subspace
Let κ be the supremum of the sectional curvaturein the√estimation. We further propose an online non-linear andball. Then, we call this ball “regular” [40] if 2r κ < π . Using
robust subspace averaging algorithm akin to online KPCA
the results in [41], we know that, for RP D with D ≥ 2,
and RPCA respectively. We also present a proof that in theκ = 1, while for general Gr(K, D) with min(K, D) ≥ 2,
limit, our proposed online robust intrinsic averaging method
0 ≤ κ ≤ 2. So, on Gr(K, D) the radius of a “regular geodesic
returns the leading robust principal components.
√ball” is < π/2 2, for min(K, D) ≥ 2 and on RP D , D ≥ 2,
the radius is < π/2.
2 A N O NLINE L INEAR S UBSPACE L EARNING A L Let X1 , . . . , XN be independent samples on Gr(K, D)
GORITHMdrawn from a distribution P (X ), then we can define an∗
In this section, we present an efficient online linear subspace intrinsic average (FM), [38], [39], M as:learning algorithm for finding the principal componentsNXof a data set. We first briefly discuss the geometry of theM∗ = argmin
d2 M, Xi .
(6)
Riemannian manifold of K -dimensional linear subspaces
M∈Gr(K,D) i=1
in RD . Then, we present an online algorithm using the
As mentioned before, on Gr(K, D), IGA exists and isgeometry of these subspaces to get the first K principal
unique if the support of P (X ) is within a “regular geodesic
components of the D-dimensional data vectors.
√ball” of radius < π/2 2 [42]. Note that for RP D , we canchoose this bound to be π/2. For the rest of the paper, we
2.1 The Geometry of Subspaces
have assumed that data points on Gr(K, D) are within a
The Grassmann manifold (or the Grassmannian) is definedas the set of all K -dimensional linear subspaces in RD and
1. These are also known as Fréchet means [38], [39].

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XXX, NO. XXX, MONTH YEAR√
“regular geodesic ball” of radius < π/2 2 unless otherwise

specified. With this assumption, the IGA is unique. Note that
this assumption is only needed for proving the theorem presented
below.
Observe, when n → ∞, FM in Eq. 6 turns out to be
Fréchet expectation (FE), [38], [39]:h iM∗ = argmin E d2 M, X .
(7)
M∈Gr(K,D)

Note that, here the expectation is with respect to the underlying density from which the samples {X } are drawn.
The IGA may be computed using a Riemannian steepest
descent, but this is computationally expensive and requires
selecting a suitable step-size [43]. Recently Chakraborty etal. [44] proposed a simple and efficient recursive/inductive
Fréchet mean estimator given by:

M1 = X1 ,1 , ∀k ≥ 1.
(8)
k+1
This approach only needs a single pass over the data set toestimate the IGA. Consequently, Eq. 8 has linear complexity
in the number of observations. Furthermore, it is truly an
online algorithm since each iteration only needs one newobservation.
Eq. 8 merely performs repeated geodesic interpolation,
which is analogous to standard recursive estimators of
Euclidean averages: Consider observations xk ∈ RD , k =1, . . . , N . Then the Euclidean average can be computed
recursively by moving an appropriate distance away from thek th estimator mk towards xk+1 on the straight line joining
xk+1 and mk . The inductive algorithm (8) for computing
the IGA works in the same way and is entirely based on
traversing geodesics in Gr(K, D) and without requiring any
optimization.
X

Mk+1 = ΓMk+1k

Theorem 1. (Weak Consistency [44]) Let X1 , . . . , XN be
samples on Gr(K, D) drawn from a distribution P (X ). Then,N
MN (Eq. 8) converges to the IGA of {Xi }i=1 (in terms of
intrinsic distance defined in Eq. 4) in probability as N → ∞.
2.3

Principal Components as Grassmann Averages

Following Hauberg et al. [7] we cast the linear dimensionality
reduction problem as an averaging problem on the Grassmannian. We consider an intrinsic Grassmann average (IGA),
which allows us to reckon with K > 1 dimensional subspaces. We then propose an online linear subspace learning andshow that for the zero-mean Gaussian data, the expected IGA
on Gr(K, D), i.e., expected K -dimensional linear subspace,coincides with the first K principal components.
Given {xi }N
i=1 , the algorithm to compute the IGA that
produces the leading K -dimensional principal subspace issketched in Algorithm 1.
Let {Xi } be the set of K -dimensional subspaces asconstructed by IGA in Algorithm 1. Moreover, assume that√the maximum principal angle between Xi and Xj is < π/2 2,∀i 6= j . The above condition is assumed to ensure that the
IGA exists and is unique on Gr(K, D). This condition canbe ensured if the angle between xl and xk is smaller than√π/2 2, ∀x , x belonging to different blocks. For x , x inlklka same block, the angle must be below π/2. Note that, this4

Algorithm 1: The IGA algorithm to compute PCs
D
Input: {xi }N
i=1 ⊂ R , K > 0
Output: {v1 , . . . , vK } ⊂ RD
N
1 Partition the data {xj }j=1 into blocks of size D × K ;th
2 Let the iblock be denoted by, Xi = [xi1 , . . . , xiK ] ;3 Orthogonalize each block and let the orthogonalized
block be denoted by Xi ;4 Let the subspace spanned by each Xi be denoted by
Xi ∈ Gr(K, D) ;∗5 Compute IGA, M , of {Xi } ;∗6 Return the K columns of an orthogonal basis of M ;these span the principal K -subspace.

assumption is needed to prove Theorem 2. In practice, even if IGA
is not unique, we find a local minimizer of Eq. 6 [38], which serves
as the principal subspace.
Theorem 2. (Relation between IGA and PCA) Let x ∼N (0, Σ), with a positive definite matrix Σ. The eigen vectorsof Σ span the Fréchet expectation as defined in Eq. 7.
Proof. Let X be the corresponding orthonormal basis of X .
Let X = [x1 , · · · , xK ], where each xi are samples drawn
from N (0, Σ). Let M = [M1 , . . . , MK ] be an orthonormal basis of an arbitrary M ∈ Gr(K, D). The squared
distance between X and M is defined as d2 (X , M) =PK
2T
= M T X be the SVD,
j=1 (arccos (Sjj )) , where, Ū S V̄and Sjj ≥ 0. Substituting into the expression for FE fromequation 7 we obtain:


K
X
2
∗
M = argmin E (arccos (Sjj ))  .
(9)
M

j=1

Now recall that the Taylor expansion of arccos(x) is given
P∞
2n+1
(2n)!
by arccos(x) = π2 − n=0 22n (n!)2 x2n+1 Substituting thisexpansion into the equation 9 for FE, we get:


K
K
X
X
2
M∗ = argmax E 
π f (Sjj ) −(f (Sjj ))  . (10)
M

j=1

j=1

P∞
(2n)! (Sjj )2n+1
where, f (Sjj ) ==n=0 22n (n!)2
2n+1 . Let sn

2
(2n)!
(2n)!
− (2n+1)2
if n is even and sn = π (2n+1)22n (n!)2 if
2n (n!)2
n is odd. Then we can rewrite the objective function to be
optimized on the RHS of equation 10 denoted by F as:


∞
K 
n/2
X
X
2

F =
sn E 
(Sjj )

=

n=1
∞
X

j=1



sn E trace

M T XX T M

n/2 

(11)

n=1

PK
2
In the last equality
above, we use=
j=1 (Sjj )
TT
trace M XX M . Note that since M is a column orthonormal matrix and hence lies on a Stiefel manifold, St(K, D),
K×K
M T Mh = I . Hence, forall,
i n > 0 and for all Λ ∈ Rn/2
trace Λ M T M − I= 0. Since (M T XX T M ) is a

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XXX, NO. XXX, MONTH YEAR

5

symmetric positiveh semi-definite matrix,
for n ≥ 2, maxi
n/2
TT
imization of trace (M XX M )with the constraint

RKHS, we will need to cope with an infinite-dimensional
Grassmannian. Fortunately, we observe that the distance
between two subspaces in RKHS is same as the distances
M T M = I is achieved by optimizing the following objective:
between span of the coefficient matrices with respect to an∞n/2  orthogonal basis. Hence, instead of performing IGA on theXF =sn E trace M T XX T M − Λ(M T M − I), subspaces in RKHS, we will perform IGA of the span ofn=1
the coefficients, which are finite dimensional. The IGA isfor some diagonal matrix Λ ∈ RK×K . To optimize F , we take then computable using the kernel trick. A key advantage ofits derivative (on the manifold) with respect to M , ∇M F := KRIGA is that it does not require an eigen decomposition of the
∂F T
∂F
−
M
M
∈ TM St(K, D), and equate it to 0. Now, Gram matrix. Furthermore, we extend this formulation to∂M∂M
propose an online KRIGA algorithm by approximating the∂F= 0. Here,observe that ∇M F = 0 if and only if ∂Mkernel function.
∞
h
X
∂F
=
sn E n(XX T M − M Λ)
3.1 Deriving the Kernel Recursive Intrinsic Grassmann∂Mn=1n/2−1 
Average (KRIGA)
M T XX T M − Λ(M T M − I)(12) Let X = {x , . . . , x }, where x ∈ RD , for all i. We seek
1

In the above equation, the interchange of expectation andderivative uses the well known dominated convergence
theorem. Rewriting the above equation 12 using Z =M T XX T M − Λ(M T M − I) , taking the norm on bothsides and then applying the Cauchy-Schwartz’s inequality,
we get
∞
X

h 

i
∂F
n/2−1
k=k
sn E n XX T M − M Λ (Z)
k
∂M
n=1
∞
h

i
X
n/2−1
=k
nsn E (Z)XX T M − M Λ k

k

n=1

≤

∞
X

 i 12h
i1 h 
n/2−1 2 2
E k XX T M − M Λ k2
n|sn |E k (Z)kn=1

Now, if M is formed from the top K eigen vectors of∂F∂FΣ = 1/kE[XX T ], then, k ∂M
k = 0, which implies, ∂M= 0,which in turn makes ∇M F = 0. Thus, M is a solution of F
in Eq. 11. This completes the proof.
Now, using Theorem 1 and Theorem 2, we replace the line
5 of the IGA Algorithm 1 by Eq. 8 to get an online subspacelearning algorithm that we call, Recursive IGA (RIGA), tocompute leading K principal components, K ≥ 1.
The key advantages of our proposed RIGA algorithm tocompute PCs are as follows:
1) In contrast to work in [7], “IGA” will return the first K
PCs, K ≥ 1.
2) Using the recursive computation of the FM in IGA leads
to RIGA, an online PC computation algorithm. Moreover,
Theorem 1 ensures the convergence of “RIGA” to “IGA”.
Hence, our proposed RIGA is an online PCA algorithm.
3) Unlike previous online PCA algorithms, RIGA is parameter free.

3

A K ERNEL E XTENSIONIn this section, we extend RIGA to perform the principal
component analysis in a Reproducing Kernel Hilbert Space(RKHS). We extend RIGA to obtain an efficient nonlinear
subspace estimator in RKHS akin to kernel PCA [9] and dub
our algorithm Kernel RIGA (KRIGA). The key issue to beaddressed in KRIGA is that, in order to perform IGA in theNi

K principal components, K ≤ D. Let K(·, ·) be the kernel
associated with RKHS H and let φ(X) = [φ(x1 ), . . . , φ(xN )],ee 1 ), . . . , φ(xe N )] be the= [φ(xwhere φ : X → H . Let φ(X)eorthogonalization of φ(X). Since each φ(xi ) ∈ Col(φ(X))
eewe have φ(xi ) = φ(X)hφ(X), φ(xi )iH , where,ehφ(X),φ(xi )iHe 1 ), φ(xi )iH , . . . , hφ(xe N ), φ(xi )iH ]t .
= [hφ(x(13)

Here, h·, ·iH is the inner productin the RKHS.


Let
Yl
=
φ(xK(l−1)+1 ), . . . , φ(xKl )=eee
Yl iH . Here,
φ(X)h
φ(X),
Yl iH , where Cl = hφ(X),e i ), φ(xK(l−1)+j )iH . Note that, Cl is a
(Cl )ij = hφ(xmatrix of dimension N × K .
Weobservethatd(Col(Yi ), Col(Yj ))d(Col(Ci ), Col(Cj )) as proved in the following lemma.

=

Lemma 1. Using the above notations, d(Col(Yi ), Col(Yj )) =d(Col(Ci ), Col(Cj )).
Proof. Yi = φ̃ (X) Ci and Yj = φ̃ (X) Cj . Now, using (4), we
can see that, d(Col(Yi ), Col(Yj )) = d(Col(Ci ), Col(Cj )) iff
the SVD of YiT Yj is same as that of CiT Cj . Now, as φ̃ (X) is
column orthogonal, hence the result follows.
So, instead of IGA on {Yl }, we can perform IGA on {Cl },where Yl = Col(Yl ) and Cl = Col(Cl ), ∀l.
e
The orthogonalization φ(X)is achieved using the GramSchimdt orthogonalization process as follows:e i ) = φ(xi ) −
φ(x

i−1
X

e j )iH φ(xe j)hφ(xi ), φ(x

(14)

j=1

e i ) = φ(xe i )/kφ(xe i )k.
φ(x(15)

e i ), φ(xK(l−1)+j )iH can be
The elements of Cl , i.e., hφ(xcomputed using the kernel K(·, ·) as given in the Lemma 2
below.
Let the basis matrix of the IGA of {Cl } be M . Then theee = φ(X)M
basis matrix of the IGA of {Yl } is denoted by U.
e
The columns of U will give the PCs in RKHS. Note that thistallies with the Representer theorem [45], which tells us that
the PC in RKHS is a linear combination of the mapped data
vectors, φ(X).
The following corollary holds by virtue of Theorem 2.

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XXX, NO. XXX, MONTH YEAR

Corollary 1. The expected IGA of {Yl } is the same as the PC of
X = {x1 , . . . , xN } in the Hilbert space H .
e is given by Proj(xi ) =The projection of xi onto Uehφ(xi ), φ(X)iM,
where,Hehφ(xi ),φ(X)iHe 1 )iH , . . . , hφ(xi ), φ(xe N )iH ].
= [hφ(xi ), φ(x(16)

The following Lemma, gives the analytic form ofe m ), φ(xi )iH .
hφ(xe m ), φ(xi )iH =Lemma 2. hφ(x

Pm−1
e
e
i )−
 K(xm ,x
j=1 hφ(xm ),φ(xj )iH hφ(xj ),φ(xi )iHq,K(xm ,xm )−

Pm−1
j=1

i≥m

e j )iHhφ(xm ),φ(x0,otherwise.

e j )iH φ(xe j)e m ) = φ(xm ) − Pm−1 hφ(xm ), φ(xProof. φ(x
j=1
e
e
e
e
and
φ(xm ) = φ(xm )/kφ(xm )k where kφ(xm )k =qPeK(xm , xm ) − m−1
Now,
since
j=1 hφ(xm ), φ(xj )iH .
Ne i )}{φ(xserves
as
a
basis
to
represent
each
φ(x
)
,
clearly,ii=1e m ), φ(xi )iH = 0 when, i < m, m = 1, · · · , N .
hφ(xSo, consider m ∈ {1, · · · , N }, i1e m ), φ(xi )iH =hφ(xe m )k (K(xm , xi )−
kφ(x

Pm−1
e j )iH hφ(xe j ), φ(xi )iH .
hφ(x
),
φ(x
m
j=1

≥

m, then,

P Note that, thus far, we have implicitly assumed thati φ(xi ) = 0, i.e., the mapped data in RKHS is centered. For
non-centered data, we have to first center the data. Given the
non-centered data, {φ(xi )}N
data be dei=1 , let the centeredP
N
1
φ̄(x
)
=
φ(x
)
−
noted by {φ̄(xi )}N
,
where
i
i
i=1
j=1 φ(xj ).
N
Then, using Lemma 2, the coefficient matrix Cl is computedusinge m ), φ̄(xi )iHhφ(xe m ), φ(xi )iH −= hφ(xN
1 X ehφ(xm ), φ(xj )iH .
N j=1

(17)

Thus, the coefficient matrices {Cl } can be computed using
only the Gram matrix K as can be seen from Eq. 17. In termsof computational complexity, KRIGA takes O(N 3 − N 2 )computations while KPCA takes O(N 3 ) computations.
Now, observe that the above KRIGA algorithm (analog
to KPCA) is not online because of the following reasons:(i) centering step of the data; (ii) choice of basis, i.e.,N{φ(xi )}i=1 . Consistent with the online KPCA algorithm, we
assume data to be centered. Then, in order to make theabove algorithm online, we need to find a predefined basis in
RKHS. We will use the idea proposed by Rahimi et al. [33], toapproximate the shift-invariant kernel K. They observed that
infinite kernel expansions can be well-approximated using
randomly drawn features. For shift-invariant K, this relates
to Bochner’s lemma [46] as stated below.
Lemma 3. K is positive definite iff K is the Fourier transform ofa non-negative measure, µ (w).

6This in turn implies the existence of a probability density
p (w) := µ (w) /C , where C is the normalizing constant.
Hence,ZK (x, y) = C exp −jwt (x − y) p (w) dw
= CEw cos wt (x − y) .
The above expectation can be approximated using
Monte Carlo methods, more specifically, we will draw
M i.i.d. RD vectors from p(w) and form matrix W
of size M × D. Then, we can approximate K(x, y)tby,p K(x, y) = ψ (W x) ψt(W y) , where ψ (W x) =C/M (cos (W x) , sin (W x)) . Now, depending on the
choice of K, p(w) will change. Forexample, for the Gaussiankx−yk2
RBF kernel, i.e., K (x, y) = exp − 2σ2. w should be



2 −1
sampled from N 0, diag σ.
Rahimi et al. [33] provided a bound on the error in
the approximation of the kernel. Now, because of thisapproximation, in our KRIGA algorithm, we can replace φ by
ψ . As ψ is finite dimensional, we will choose the canonicale in the above derivation bybasis in R2M , i.e., replacing φ2M{ei }i=1 . This gives us an online KRIGA algorithm using therecursive IGA.

4A R OBUST O NLINE L INEAR S UBSPACE L EARN ING A LGORITHM
In this section, we will propose an online robust
PCA algorithm using intrinsic Grassmann averages. Let{X1 , X2 , · · · , XN } ⊂ Gr(K, D),√K < D be inside a regular
geodesic ball of radius < π/2 2 s.t., the Fréchet Median(FMe) [47] exists and is unique. Let X1 , X2 , · · · , XN be thecorresponding orthonormal basis, i.e., Xi spans Xi , for all i.
The FMe can be computed via the following minimization:M∗ = arg min
MNXd(Xi , M).

(18)

i=1

With a slight abuse of notation, we use the notation M∗ (M )to denote both the FM and the FMe (and their orthonormal
basis). The FMe is robust as was shown by Fletcher et al.
[47], hence we call our estimator Robust IGA (RoIGA). In thefollowing theorem, we will prove that RoIGA leads to therobust PCA in the limit as the number of the data samples
goes to infinity. An algorithm to compute RoIGA is obtainedby simply replacing Step 5 of Algorithm 1 by computationof RoIGA via minimization of Eq. 18 instead of Eq. 6. This
minimization can be achieved using the Riemannian steepest
descent, but instead, here we use the stochastic gradientdescent of batch size 5 to compute RoIGA. As at each
iteration, we need to store only 5 samples, the algorithm
is online. The update step for each iteration of the online
algorithm to compute RoIGA (we refer to our online RoIGA
algorithm as Recursive RoIGA (RRIGA)) is as follows:

M 1 = X1 ,Mk+1 = ExpMk

Exp−1
Mk (Xk+1 )(k + 1)d(Mk , Xk+1 )

!

.

(19)

where, k ≥ 1, Exp and Exp−1 are Riemannian Exponential
and inverse Exponential functions as defined below.

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XXX, NO. XXX, MONTH YEAR

Definition 2 (Exponential map). Let X ∈ Gr(K, D). Let
B (0, r) ⊂ TX Gr(K, D) be an open ball centered at the origin
in the tangent space at X , where r is the injectivity radius[38] of Gr(K, D). Then, the Riemannian Exponential map is
a diffeomorphism ExpX : B (0, r) → Gr(K, D).
Definition 3 (Inverse Exponential map). Since, inside B (0, r),
Exp is a diffeomorphism, hence the inverse Exponential map−1
is defined and is a map
 ExpX : U → B (0, r), whereU = ExpX (B (0, r)) := ExpX (U ) |U ∈ B (0, r) .
We refer the readers to [48] for the consistency proof ofthe estimator. Notice that, Eq. 19 can be rewritten as,

M1 = X1 ,X

Mk+1 = ΓMk+1(γk ) .
k(20)

1
Where, γk = (k+1)d(M. This is because, ΓYX (t) cank ,Xk+1 )−1−1
be written as ExpX tExpX (Y) , where, ExpX (Y) is thevelocity vector when moving from X to Y and ExpX (U ) isa point on Gr(K, D) which can be reached by going fromX along the shortest geodesic given by the velocity vector
U . Let Mk and Xk+1 be orthonormal basis of Mk and Xk+1
respectively. Then, the shortest geodesic can be expressed as:

70, then we can claim that our objective function in Eq. 18 isrobust [50].
This can be justified by noting that in the presence ofoutliers, i.e., when m → ∞ as variance becomes larger, wewant no change in the objective function, i.e., the gradientwith respect to m should be zero. This is due to the fact that,the other terms in the Taylor expansion
bounded√PK are upper
qP
2
∂E[
K
j=1 ((Si )jj ) ]2bygoes to 0,
j=1 ((Si )jj ) , hence, as∂mso do the gradients of the other terms.
φ(m+0.5)−φ(m)Now, ψ(m) = Γ(m)Γ(m + 0.5)
, where φ isΓ(m)2
the polygamma function [51] of order 0. After some simple
calculations, we get,lim (φ(m + 0.5) − φ(m)) = lim log(1 + 1/(2m))
m→∞


∞ 
X
1
1
−
+ lim
Bk
m→∞
k mk
k(m + 0.5)k
k=1

m→∞

= lim log(1 + 1/(2m)) + 0 = 0.
m→∞

Here, {Bk } are the Bernoulli numbers of the second kind[52]. Hence, limm→∞ ψ(m) = 0.


We would like to point out that the outlier corrupted datacanbe modeled using a mixture of independent randomXΓMk+1(t) = Col (Mk V cos (tS) + U sin (tS)) ,(21) variables, Y1 , Y2 , where Y1 ∼ N (0, Σ1 ) (to model non-outlierk−1data samples) and Y2 ∼ N (µ, Σ2 ) (to model outliers), i.e.,
where, U ΣV T = Xk+1 MkT Xk+1
− Mk is the thin(∀i), xi = w1 Y1 + (1 − w1 )Y2 , w1 > 0 is generally large,
singular value decomposition and S = atan (Σ).
so that the probability of drawing outliers is low. Then, as
Theorem 3. (Robustness of RoIGA) Assuming the above the mixture components are independent, (∀i), xi ∼ N ((1 −hypotheses and notations, as N → ∞, the columns of M are w1 )µ, w12 Σ1 +(1−w1 )2 Σ2 ). A basic assumption in any online
robust to outliers, where M is the orthonormal basis of M∗ as PCA algorithm is that data is centered. So, in case the data
are not centered (similar to the model of xi ), the first step
defined in Eq. 18.
of PCA would be to center the data. But then the algorithm
Proof. Let, Xi = [xi1 , · · · , xiK ] and xij be i.i.d. samples
cannot be made online, hence our above assumption thatdrawn from N (0, Σ), with a positive definite matrix Σ.
xi ∼ N (0, Σ) is a common assumption in an online scenario.
Let, M = [M1 , · · · , MK ] be an orthonormal basis of M.
But, in a general case, after centering the data as the first step
Recall that theqdistance between Xi and M is defined asof PCA, the above theorem is valid.
PKT
2d(Xi , M) ==j=1 (arccos((Si )jj )) , where Ūi Si Vi
M T Xi is the SVD, and (Si )jj ≥ 0.
5 E XPERIMENTAL R ESULTS
Observethat Eq. 18 is equivalent to maximizingq
PN
PK
We evaluate the performance of the proposed recursive2
i=1
j=1 (arcsin((Si )jj )) . Here, we will only foestimators on both real and synthetic data. Our overallcus on the first term of the Taylor
expansion (whichfindings are that the RIGA estimator is more accurate thanPN qPK
bounds the other terms) of i=1(arcsin((Si )jj ))2 , other online linear subspace estimators since it is parameter
j=1
qP
PK
∞
free. The Kernel RIGA (KRIGA) is found to yield results
2ni.e.,
n=1 tnj=1 ((Si )jj ) , where {tn } denote thecoefficients in the Taylor expansion.
The first term in that are almost identical to Kernel PCA (KPCA) but at aqPsignificant reduction in run time. Below we consider RIGAK
the expansion is Ve =((Si )jj )2 . Note that, and KRIGA separately.
j=1
 P

PK
PK
K
1
2
2
2
j=1 ((Si )jj ) ∼ Γ 2
j=1 σMj , 2 , as,
j=1 ((Si )jj ) =trace(M T Xi XiT M ), and M T Xi ∼ N (0, ΣM ). Hence, Ve 5.1 Online Linear Subspace Estimation
PKP
2
2
follows Ng 21 j=1 σM, K
j=1 σMj , where Ng is the Nak- Baselines: We compare with Oja’s rule and and the onlinejagami distribution∞, the RHS of Eq.
q [49]. Now, asN →q 18 version of EM PCA (Sec. 1.1). For Oja’s rule we follow√
PKPK
common guidelines and consider step-sizes γt = α/D t
22becomes E
((S
)
)
.
E
((S
)
)
=
i
jj
i
jj
j=1
j=1
with α-values between 0.005 and 0.2. For EM PCA we
√
PK
PK
2
2
2Γ( j=1 σŪij Mj + 0.5)/Γ( j=1 σŪij Mj ), where Γ is the follow the recommendations in Cappé [13] and use stepq
sizes γt = 1/tα with α-values between 0.6 and 0.9 along withPK2
well-known gamma function. Thus, E= Polyak-Ruppert averaging.
j=1 ((Si )jj )√PK(Synthetic) Gaussian Data: Theorem 2 state that the2ρ(m) , 2Γ(m + 0.5)/Γ(m), where m = j=1 σMj . Hence,
RIGA estimates coincide in expectation with the leadingthe √influence function [50] of ρ is proportional to ψ(m) , principal subspace when the data is drawn from a zeroPK2∂E[
j=1 ((Si )jj ) ]and if we can show that limm→∞ ψ(m) = mean Gaussian distribution. We empirically verify this for∂mIEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XXX, NO. XXX, MONTH YEAR

8

0.9

0.8
0.6
0.4
RIGA
GROUSEEM PCA
Oja's rule

0.2

0.85

0.8
0.6
0.4
RIGA
GROUSEEM PCA
Oja's rule

0.21000 2000 3000 4000 5000 6000 7000 8000 9000 10000

0.81000 2000 3000 4000 5000 6000 7000 8000 9000 10000

Number of Observations

4000

5000

Fig. 2. Expressed variance as a function of number of observations. Left:
The mean and one standard deviation of the RIGA estimator computed
over 150 trials. In each trial data are generated in R50 and we estimate
a K = 2 dimensional subspace. Right: The performance of different
estimators for varying step-sizes. Here data are generated in R250 and
we estimate a K = 20 dimensional subspace.

an increasing number of observations drawn from randomly
generated zero-mean and 0.5I variance Gaussian. We measure the expressed variance which is the ratio of the variancecaptured by the estimated subspace to the variance capturedby the true principal subspace:
Expressed Variance =

K
X

PN

T (est)n=1 xn vk∈ [0, 1].
PN
T (true)
n=1 xn vk
k=1

(22)

An expressed variance of 1 implies that the estimated subspace captures as much variance as the principal subspace.
The right panel of Fig. 2 shows the mean (± one standard
deviation) expressed variance of RIGA over 150 trials. It isevident that for the Gaussian data, the RIGA estimator does
indeed converge to the true principal subspace.
A key aspect of any online estimator is that it should
be stable and converge fast to a good estimate. Here, we
compare RIGA to the above-mentioned baselines. Both Oja’s
rule and EM PCA require a step-size to be specified, so we
consider a larger selection of such step-sizes. The left panel ofFig. 2 shows the expressed variance as a function of numberof observations for different estimators and step-sizes. In Fig.
3, we have comparative performance analysis of EM PCA,GROUSE, Oja’s rule and RIGA. EM PCA was found to bequite stable with respect to the choice of step-size, though it
does not seem to converge to a good estimate. Oja’s rule, on
the other hand, seems to converge to a good estimate, but its
practical performance is critically dependent on the step-size(as evident from Fig. 2). GROUSE is seen to oscillate for small
data size however, with a large number of samples, it yields
a good estimate. On the other hand, RIGA is parameter free
and is observed to have good convergence properties. Thisbehavior of RIGA is consistent for D ≥ 100 and K ≥ 10 asobserved empirically and depicted in Fig. 4.
In the right panel of Fig. 5, we perform a stability analysis
of GROUSE and RIGA. Here, for a fixed value of N , wegenerate a data matrix and perform 200 independent runson the data matrix and report the mean (± one standard
deviation) expressed variance. On the left and middle panels,
we show the performance of GROUSE with varying D and
K for learning rates 0.0001 and 0.01 respectively. We cansee that with a larger learning rate, GROUSE can achieve
better expressed variance but with less stability. As can beseen from the figure, RIGA is very stable in comparison to
GROUSE.

D = 100

1

D = 150

1

0.98

0.98
K=2
K=7
K = 12
K = 17
K = 22

0.96

0.94
2000

4000

6000

8000

0.94
2000

10000

D = 200

1

K=2
K=7
K = 12
K = 17
K = 22

0.96

4000

6000

0.98

8000

10000

D = 250

1

0.98
K=2
K=7
K = 12
K = 17
K = 22

0.96

0.94
2000

4000

6000

8000

K=2
K=7
K = 12
K = 17
K = 22

0.96

0.94
2000

10000

4000

Number of Observations

6000

8000

10000

Number of Observations

Fig. 4. Performance of RIGA with varying D and K . We can see that formoderately large D and K ≥ 10, the performance of RIGA is very good.

In the rest of the paper, we will use average reconstruction
error (ARE) [9] to measure the “goodness” of the estimatedsubspace, which is defined as follows:N

Average Reconstruction Error =1 X
e n k2 , (23)kxn − x
N n=1

e is the reconstructed sample using the estimated
where, xK
principal subspace spanned by {vk }k=1 .
1

1
0.9

0.9

Expressed Variance

3000

0.8

0.7

0.6

0.5

0.4

0.3
1000

0.8
0.7
0.6
0.5
0.4
0.3
0.2

2000

3000

4000

5000

6000

7000

8000

9000

10000

0.1
1000

2000

Number of Observations

Expressed Variance

2000

Number of Observations

Expressed Variance

1000

Number of Observations

Fig. 3. Expressed variance as a function of number of observations. The
performance of different estimators. Left: Data are generated in R250
and we set K = 20. Right: Data are generated in R100 and we set
K = 10. We observe that our estimator is better than its competitors forother values of D > 100 and K ≥ 10.

Expressed Variance

0.75

Expressed Variance

Expressed Variance

0.95

Expressed Variance

Expressed Variance

1

3000

4000

5000

6000

7000

8000

9000

10000

Number of Observations

0.95
GROUSE
RIGA

0.9

0.85

2000

4000

6000

8000

10000

Number of Observations

Fig. 5. Left and Middle: Performance of GROUSE with varying D and K
and learning rates of 0.0001 and 0.01 respectively. D ranges between
50 and 500 in increments of 50 and K ranges between 2 and 22 inincrements of 5. The plots are color coded from hot to cold, with warmest
being D = 50, K = 2 and coolest being D = 500, K = 22. We cansee that though for relatively larger learning rate the performance of
GROUSE is better, it is not stable. Right: Stability analysis comparisonof GROUSE and RIGA (for a fixed N , we randomly generate a data
matrix, X , from a Gaussian distribution on R250 . We estimate K = 20
dimensional subspace and report the mean and standard deviation over
200 runs on X with learning rate 0.0001. This plot is a close-up look at
the leftmost plot with D = 250 and K = 20.)

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XXX, NO. XXX, MONTH YEAR

Human Body Shape: Online algorithms are generally
well-suited for solving large-scale problems as they, byconstruction, should have linear time-complexity in thenumber of observations. As an example we consider a large
collection of three-dimensional scans of human body shape[53]. This dataset contains N = 21862 meshes which eachconsist of 6890 vertices in R3 . Each mesh is, thus, viewed
as a D = 6890 × 3 = 20670 vector. We estimate a K = 10
dimensional principal subspace using Oja’s rule, EM PCA
and RIGA respectively. The average reconstruction error over
all meshes are 16.8 mm for Oja’s rule, 1.9 mm for EM PCA,
and 1.0 mm for RIGA. Note that both Oja’s rule and EM PCA
explicitly minimize the reconstruction error, while RIGA does not
but yet outperforms the baseline methods. We speculate that
this is due to RIGA’s excellent convergence properties and it
being a parameter free algorithm is not bogged down by the
hard problem of step-size tuning confronted in the baseline
algorithms used here.
Santa Claus Conquers the Martians: We now consider
an even larger scale experiment and consider all frames ofthe motion picture Santa Claus Conquers the Martians (1964)2 .
This consist of N = 145, 550 RGB frames of size 320 × 240,
corresponding to an image dimension of D = 230, 400. We
estimate a K = 10 dimensional subspace using Oja’s rule,
EM PCA and RIGA respectively. Again, we measure theaccuracy of the different estimators via the reconstruction
error. Pixel intensities are scaled to be between 0 and 1.
Oja’s rule gives an average reconstruction error of 0.054,
EM PCA gives 0.025, while RIGA gives 0.023. Here RIGA
and EM PCA gives roughly equally good results, with aslight advantage to RIGA. Oja’s rule does not fare as well.
As with the shape data, it is interesting to note that RIGA
outperforms the other baseline methods on the error measurethat they optimize even though RIGA optimizes a different
measure.
5.2

Nonlinear Subspace Estimation

We now analyze comparative performance of the proposed
online KRIGA with KPCA as the baseline. In our experiments,
we use a Gaussian kernel with σ = 1. The performanceis compared in terms of the time required and averagereconstruction error (ARE). In Fig. 6, we present a synthetic
experiment, where the data generated is in the form of three
concentric circles. We can see that both KPCA and online
KRIGA yield similar cluster separation. As expected, we
can see that with very few dimensions for approximation(i.e., with small M ), the performance of online KRIGA ispoor. Similar observation can be made from Fig. 7, wherewith M = 500, we get almost as good result as KPCA. We
also compared KRIGA’s performance with SKPCA [35] andthough in Fig. 6, the results using SKPCA are worse thanours, the results in Fig. 7 are however comparable.
Now, we assess the performance of KRIGA and KPCA interms of ARE and computation time, based on randomly
generated synthetic data. Here, we compare our online
KRIGA with KPCA. In order to make a fair comparison,
we have used the MATLAB ‘eigs’ function of KPCA which issignificantly faster than KPCA. The results for ARE andcomputation time are shown in Fig. 8. We can see that
2. https://archive.org/details/SantaClausConquerstheMartians1964(a) synthetic data(b) KPCA(d) online KRIGA (e) SKPCA M=500(with M=50)9(c) online KRIGA(with M=500)(f) SKPCA M=50

Fig. 6. Results from KRIGA and SKPCA on synthetic data.

(a) synthetic data(b) KPCA(d) online KRIGA (e) SKPCA M=500(with M=50)(c) online KRIGA(with M=500)(f) SKPCA M=50

Fig. 7. Results from KRIGA and SKPCA on 3D synthetic data.

our online KRIGA is faster than KPCA with ‘eigs’ without
sacrificing much ARE. For this experiment, we chose K = 5.
Though, the ARE of KPCA is better than that of KRIGA,
we can see from Fig. 8 that with increasing number of PCs,
performance of KRIGA is similar to KPCA.
Finally, we test the KPCA and the online KRIGA algorithms on the entire movie, Santa Claus Conquers the Martians(1964) samples at 10 FPS, and present the time comparisonin Fig. 8. We observe a cubic time growth for KPCA while
for the online KRIGA, the time is almost a constant. This
demonstrates the scalability of our proposed method.

5.3

Robust Subspace Estimation

We now present a comparative experimental evaluation ofrobust extension (RRIGA). Here we use several baseline
methods and measure the performance using the reconstruction error (RE). We use UCSD anomaly detection database[54] and the Extended YaleB database [55]. Before presenting
these experiments, we present a synthetic experiment to show
comparisons of several robust and non-robust algorithms in asimulated setting. In this simulated setting, we demonstrate
the necessity of robust PCA algorithm in the case of anincreased amount of noise present in the data.

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XXX, NO. XXX, MONTH YEAR
8

600

10

KPCA(eigs)
KRIGA

Time(s)500

KPCA(eigs)
KRIGA

KPCA(eigs)
KRIGA

Time(s)Avg. reconstr. error in RKHS

Avg. reconstr. error in RKHS

are walking towards the camera) and “Peds2” (people arewalking
parallel to the camera plane). We will only consider
400
6“Peds1”
inthis experiment. In “Peds1” there are 36 training
300
5
200
4
and 34 testing videos where each video contains 200 frames
100
3
of dimension 158 × 238 (D = 37604). The test frames do not
0
2
0
20
40
60
80
100
0
20
40
60
80
100
have any anomalous activities. Some sample frames (with# of PCs# of PCs
and without outliers) are shown in Fig. 10. We first extract(a) (Left) ARE and (Right) time with number of PCs
K principal components on the training data (including
30
0.6
anomalies) and then compute reconstruction error on the test
25
0.55
frames (without anomalies) using the computed principal
20
0.5
15
components. It is expected that if the PC computation
0.45
10technique is robust, the reconstruction error will be good
0.4
5
since PCs should not be affected by the anomalies in training
0.35
0
0
5000
10000
15000
0
5000
10000
15000
samples. In Figs. 11,12, we compare the performance of# of samples# of samples
RRIGA with GRASTA, DHR-PCA, IALM, REPROCS in terms(b) (Left) ARE and (Right) time with number of samplesof RE and the (computation) time required. In terms of
100time it is evident that RRIGA is very fast compared to the80
competitors. RRIGA also outperforms the state-of-the-art in60
terms of the RE.
40
Yale ExtendedB database: This data contains 2414 face
20
imagesof 38 subjects. Each image was cropped to a 32 × 32
0
0
0.5
1
1.5
2
# of frames
10
image (D = 1024). Due to varying lighting conditions, some(c) Time with number of framesof the face images are shaded/ dark and can be treated as
outliers (this experimental setup is similar to the one in [57]).
Fig. 8. Comparison of KPCA and KRIGA in terms of ARE and running
In Fig. 10 some sample face images (outliers and non-outliers)
time.
are shown. One can see that due to poor lighting condition,
though the face in the middle of the top row is a face image,it appears completely dark and as an outlier. For testing, we
5.3.1 Synthetic experiment
We follow the exact same setup for the synthetic experi- use 142 non-outlier face images of 38 subjects and the rest
ment as in [56]. We choose D = 200 and K = 10 and are used to extract the PCs. We report RE (with varying K )select the ground truth subspace as U ∗ = Col (U ∗ ), where and computation time required for RRIGA, GRASTA, DHRU ∗ is composed of i.i.d. samples from standard normal PCA, IALM, REPROCS methods in Figs. 13, 12 respectively.
distribution. The coefficients are drawn from a standard As evident, RRIGA is faster than the competitors while
normal distribution. The samples of the data are generated outperforming all of the state-of-the-art methods except
REPROCS in terms of reconstruction error.
by adding a Gaussian noise with zero mean and σ standard
deviation. σ is set to lie in the range 10−2 , 10−5 , 0 . Given ane of the estimated subspace, we computeorthonormal basis, UeUe T U ∗ k. We compare 6 C ONCLUSIONS
the projected error using, kU ∗ − U
our results with an extensive sets of algorithms including In this paper, we present a geometric framework to compute
MDISVD [22], Brand’s [23], PIMC [24], PETRELS [25], IALM principal linear subspaces in finite and infinite dimensional[20], GROUSE and Oja’s algorithm. We should mention that reproducing kernel Hilbert spaces (RKHS). We compute an
for IALM requires the number of observations to be ≥ D, intrinsic Grassmann average as a proxy for the principal
which is the reason for using a constant value for the first linear subspace and show that if the samples are drawn fromD − 1 observations. The performance is depicted in Fig. 9, a Gaussian distribution, the intrinsic Grassmann average
which clearly indicates that RRIGA outperforms others in coincides with the principal subspace in expectation. We
the comparisons.
further show that the approach extends to the RKHS setting.
A robust version of the online PCA is also presented along
5.3.2 Real experimentwith several experiments demonstrating its performance inNow, we present experiments on real datasets. In these set comparison to the state-of-the-art. The approach has severalof experiments we compare performance of RRIGA with advantages. Unlike the work by Hauberg et al. in [7], ourGRASTA, DHR-PCA, IALM and REPROCS. As mentioned estimator returns the first K ≥ 1 components. The proposed
earlier, we present results of experiments performed on algorithm is inherently online, which also makes it scalable(a) the UCSD anomaly detection database, (b) the Extended to large datasets. We have demonstrated this by performing
YaleB database and the (c) Wallflower database.
principal component analysis of an entire Hollywood movie.
UCSD anomaly detection database: This data contains Unlike most other online algorithms there are no step-sizes orimages of pedestrian movement on walkways captured other parameters to tune; a very useful property in practicalby a stationary mounted camera. The crowd density on settings. We extend the approach to RKHS and thereby
the walkway varies from sparse to very crowded. The provided an algorithm that serves the same purpose as kernel
anomaly includes bikers, skaters, carts, people in wheelchair PCA. A benefit of our formulation is that, unlike KPCA, ouretc. This database is divided in two sets: “Peds1” (people estimator does not require an eigen decomposition of the
KPCA(eigs)
KRIGA

7

Time(s)KPCA(eigs)
KRIGA

4

Reconstruction Error

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XXX, NO. XXX, MONTH YEAR

noise var=1e-02

101

noise var=1e-05

101

noise var=0e+00
10 1
10 6

100

10 1

10 1

10 3

10 2

10 5

10 11
10 16

10 3

10 7

10 21

10 4

10 9

10 26

0

200

400

600

800 1000

0

200

400

600

800 1000

Number of Observations

0

200

400

600

800 1000

11

GROUSEMDISVD
Brand
PIMC
PETRELS
Oja
IALM
RRIGA

Fig. 9. Comparative analysis for different noise levels. For all the competing methods, we use the suggested parameters in [56]Fig. 10. Top and bottom row contain outliers (identified in a rectangularbox) and non-outliers frames of UCSD and YaleExtendedB data respectively.

Gram matrix. Empirically, we observe that our algorithm issignificantly faster than KPCA while giving similar results.

ACKNOWLEDGEMENTS
We thank Chun-Hao Yang for several helpful suggestions.
This research was in part funded in part by the NSF grants
IIS-1525431 and IIS-1724174 to BCV. SH was supported by aresearch grant (15334) from VILLUM FONDEN. This project
has received funding from the European Research Council(ERC) under the European Union’s Horizon 2020 researchand innovation programme (grant agreement no 757360).

R EFERENCES[1]

K. Pearson, “On lines and planes of closest fit to system of points
in space,” Philosophical Magazine, vol. 2, no. 11, pp. 559–572, 1901.
[2] H. Hotelling, “Analysis of a complex of statistical variables into
principal components.” Journal of educational psychology, vol. 24,
no. 6, p. 417, 1933.
[3] M. E. Tipping and C. M. Bishop, “Probabilistic principal component
analysis,” Journal of the Royal Statistical Society: Series B (Statistical
Methodology), vol. 61, no. 3, pp. 611–622, 1999.
[4] H. Zou, T. Hastie, and R. Tibshirani, “Sparse principal component
analysis,” Journal of computational and graphical statistics, vol. 15,
no. 2, pp. 265–286, 2006.
[5] R. Vidal, Y. Ma, and S. Sastry, “Generalized principal component
analysis (gpca),” IEEE transactions on pattern analysis and machine
intelligence, vol. 27, no. 12, pp. 1945–1959, 2005.
[6] I. T. Jolliffe and J. Cadima, “Principal component analysis: a review
and recent developments,” Phil. Trans. R. Soc. A, vol. 374, no. 2065,
p. 20150202, 2016.
[7] S. Hauberg, A. Feragen, R. Enficiaud, and M. J. Black, “Scalable
robust principal component analysis using grassmann averages,”
IEEE Transactions on Pattern Analysis and Machine Intelligence(TPAMI), 2015.
[8] R. Chakraborty, S. Hauberg, and B. C. Vemuri, “Intrinsic grassmann
averages for online linear and robust subspace learning,” inProceedings IEEE Conf. on Computer Vision and Pattern Recognition(CVPR), July 2017.
[9] B. Schölkopf, A. Smola, and K.-R. Müller, “Kernel principal
component analysis,” in Artificial Neural Networks. Springer, 1997,
pp. 583–588.
[10] E. Oja, “Simplified neuron model as a principal component
analyzer,” Journal of mathematical biology, vol. 15, no. 3, pp. 267–
273, 1982.

[11] Z. Allen-Zhu and Y. Li, “First efficient convergence for streaming
k-PCA: a global, gap-free, and near-optimal rate,” in Proceedingsof the IEEE Symposium on Foundations of Computer Science (FOCS),
Berkely, CA, october 2017.
[12] S. Roweis, “EM algorithms for pca and spca,” Advances in neural
information processing systems, pp. 626–632, 1998.
[13] O. Cappé, “Online expectation-maximisation,” Mixtures: Estimation
and Applications, pp. 31–53, 2011.
[14] L. Balzano, R. Nowak, and B. Recht, “Online identification andtracking of subspaces from highly incomplete information,” inCommunication, Control, and Computing (Allerton), 2010 48th Annual
Allerton Conference on. IEEE, 2010, pp. 704–711.
[15] J. He, L. Balzano, and A. Szlam, “Incremental gradient on thegrassmannian for online foreground and background separation in
subsampled video,” in CVPR, 2012, pp. 1568–1575.
[16] S. Boyd and L. Vandenberghe, Convex optimization. Cambridge
university press, 2004.
[17] C. Boutsidis, D. Garber, Z. Karnin, and E. Liberty, “Online principal
components analysis,” in Proceedings of the Twenty-Sixth Annual
ACM-SIAM Symposium on Discrete Algorithms. SIAM, 2015, pp.
887–901.
[18] E. J. Candès, X. Li, Y. Ma, and J. Wright, “Robust principal
component analysis?” Journal of the ACM (JACM), vol. 58, no. 3,
p. 11, 2011.
[19] B. Lois and N. Vaswani, “Online robust pca and online matrix
completion,” arXiv preprint arXiv:1503.03525, 2015.
[20] Z. Lin, M. Chen, and Y. Ma, “The augmented lagrange multipliermethod for exact recovery of corrupted low-rank matrices,” arXiv
preprint arXiv:1009.5055, 2010.
[21] J. R. Bunch and C. P. Nielsen, “Updating the singular value
decomposition,” Numerische Mathematik, vol. 31, no. 2, pp. 111–
129, 1978.
[22] R. Kennedy, L. Balzano, S. J. Wright, and C. J. Taylor, “Online algorithms for factorization-based structure from motion,” Computer
Vision and Image Understanding, vol. 150, pp. 139–152, 2016.
[23] M. Brand, “Incremental singular value decomposition of uncertain
data with missing values,” in European Conference on Computer
Vision. Springer, 2002, pp. 707–720.
[24] R. Kennedy, C. J. Taylor, and L. Balzano, “Online completion of
ill-conditioned low-rank matrices,” in 2014 IEEE Global Conferenceon Signal and Information Processing (GlobalSIP). IEEE, 2014, pp.
507–511.
[25] Y. Chi, Y. C. Eldar, and R. Calderbank, “PETRELS: Parallel subspaceestimation and tracking by recursive least squares from partial
observations,” IEEE Transactions on Signal Processing, vol. 61, no. 23,
pp. 5947–5959, 2013.
[26] B. Yang, “Projection approximation subspace tracking,” IEEE
Transactions on Signal processing, vol. 43, no. 1, pp. 95–107, 1995.
[27] H. Guo, C. Qiu, and N. Vaswani, “An online algorithm forseparating sparse and low-dimensional signal sequences from theirsum,” IEEE Transactions on Signal Processing, vol. 62, no. 16, pp.
4284–4297, 2014.
[28] W. Ha and R. F. Barber, “Robust pca with compressed data,” inAdvances in Neural Information Processing Systems, 2015, pp. 1936–
1944.
[29] J. Feng, H. Xu, and S. Yan, “Online robust pca via stochastic
optimization,” in Advances in Neural Information Processing Systems,
2013, pp. 404–412.
[30] ——, “Robust pca in high-dimension: A deterministic approach,”
arXiv preprint arXiv:1206.4628, 2012.

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XXX, NO. XXX, MONTH YEAR

12

Fig. 11. Comparison of reconstructions of a random test sample. (Left-to-Right) Original sample, RRIGA, GRASTA, DHR-PCA, IALM, REPROCS
outputs respectively. For GRASTA, default parameters suggested in [15] are used. An average reconstruction error of 9.51, 18.02, 21.46, 30.67, 21.01
is obtained for RRIGA, GRASTA, DHR-PCA, IALM, REPROCS respectively.

6000

Time(s)

4000

50
DHR-PCAGRASTAREPROCS
IALM
RRIGA

40

Time(s)

5000

3000

DHR-PCAGRASTAREPROCS
IALM
RRIGA

30

20

2000
10

1000
0

10

20

30

40

50

60

Number of PCs(K)

70

80

90

100

0

10

20

30

40

50

60

70

80

90

100

Number of PCs(K)Fig. 12. Computation time required for “Peds1” (Left) and “YaleB” (Right)
respectively.

Fig. 13. Comparison of reconstruction of random test samples. (Leftto-Right) Original sample, RRIGA, GRASTA, DHR-PCA, IALM and
REPROCS outputs respectively. Default parameter settings suggested in[15] are used for GRASTA. An average reconstruction error of 2.91, 4.85,
3.80, 3.03, 1.65 is obtained for RRIGA, GRASTA, DHR-PCA, IALM and
REPROCS respectively.
[31] X. Xu, “Online robust principal component analysis for backgroundsubtraction: A system evaluation on toyota car data,” Master’s
thesis, University of Illinois at Urbana-Champaign, 2014.
[32] C. Qiu, N. Vaswani, B. Lois, and L. Hogben, “Recursive robust pca
or recursive sparse recovery in large but structured noise,” IEEE
Transactions on Information Theory, vol. 60, no. 8, pp. 5007–5039,
2014.
[33] A. Rahimi and B. Recht, “Random features for large-scale kernel
machines,” in Advances in neural information processing systems, 2008,
pp. 1177–1184.
[34] P. Honeine, “Online kernel principal component analysis: A
reduced-order model,” IEEE transactions on pattern analysis andmachine intelligence, vol. 34, no. 9, pp. 1814–1826, 2012.
[35] M. Ghashami, D. J. Perry, and J. Phillips, “Streaming kernel
principal component analysis,” in Artificial Intelligence and Statistics,
2016, pp. 1365–1374.
[36] K. I. Kim, M. O. Franz, and B. Scholkopf, “Iterative kernel principal
component analysis for image modeling,” IEEE transactions on
pattern analysis and machine intelligence, vol. 27, no. 9, pp. 1351–1366,
2005.
[37] S. Günter, N. N. Schraudolph, and S. Vishwanathan, “Fast iterative
kernel principal component analysis,” Journal of Machine Learning
Research, vol. 8, no. Aug, pp. 1893–1918, 2007.
[38] H. Karcher, “Riemannian center of mass and mollifier smoothing,”
Communications on pure and applied mathematics, vol. 30, no. 5, pp.
509–541, 1977.
[39] M. Fréchet, “Les éléments aléatoires de nature quelconque dans
un espace distancié,” in Annales de l’institut Henri Poincaré, vol. 10.
Presses universitaires de France, 1948, pp. 215–310.
[40] W. S. Kendall, “Probability, convexity, and harmonic maps withsmall image i: uniqueness and fine existence,” Proceedings of the
London Mathematical Society, vol. 3, no. 2, pp. 371–406, 1990.
[41] Y.-C. Wong, “Sectional curvatures of grassmann manifolds,” Proceedings of the National Academy of Sciences, vol. 60, no. 1, pp. 75–79,
1968.
[42] B. Afsari, “Riemannian Lp center of mass: existence, uniqueness,and convexity,” Proceedings of the American Mathematical Society, vol.
139, no. 2, pp. 655–673, 2011.
[43] X. Pennec, “Intrinsic statistics on riemannian manifolds: Basic tools
for geometric measurements,” Journal of Mathematical Imaging andVision, vol. 25, no. 1, pp. 127–154, 2006.
[44] R. Chakraborty and B. C. Vemuri, “Recursive frechet mean computation on the grassmannian and its applications to computer vision,”in Proceedings of the IEEE International Conference on Computer Vision,
2015, pp. 4229–4237.
[45] B. Schölkopf, R. Herbrich, and A. J. Smola, “A generalized
representer theorem,” in Computational learning theory. Springer,
2001, pp. 416–426.
[46] W. Rudin, Fourier analysis on groups. Courier Dover Publications,
2017.
[47] P. T. Fletcher, S. Venkatasubramanian, and S. Joshi, “The geometric
median on riemannian manifolds with application to robust atlasestimation,” NeuroImage, vol. 45, no. 1, pp. S143–S152, 2009.
[48] S. Bonnabel, “Stochastic gradient descent on riemannian manifolds,”
IEEE Transactions on Automatic Control, vol. 58, no. 9, pp. 2217–2229,
2013.
[49] M. Nakagami, “The m-distribution-a general formula of intensity
distribution of rapid fading,” Statistical Method of Radio Propagation,
1960.
[50] P. J. Huber, Robust statistics. Springer, 2011.
[51] M. Abramowitz, I. A. Stegun et al., “Handbook of mathematical
functions,” Applied mathematics series, vol. 55, p. 62, 1966.
[52] S. Roman, The umbral calculus. Springer, 2005.
[53] G. Pons-Moll, J. Romero, N. Mahmood, and M. J. Black, “Dyna: A
model of dynamic human shape in motion,” ACM Transactions on
Graphics, vol. 34, no. 4, pp. 120:1–120:14, Aug. 2015.
[54] V. Mahadevan, W. Li, V. Bhalodia, and N. Vasconcelos, “Anomaly
detection in crowded scenes.” in CVPR, vol. 249, 2010, p. 250.
[55] A. S. Georghiades, P. N. Belhumeur, and D. J. Kriegman, “Fromfew to many: Illumination cone models for face recognition undervariable lighting and pose,” IEEE transactions on pattern analysis andmachine intelligence, vol. 23, no. 6, pp. 643–660, 2001.
[56] L. Balzano, Y. Chi, and Y. M. Lu, “Streaming pca and subspacetracking: The missing data case,” Proceedings of the IEEE, vol. 106,
no. 8, pp. 1293–1310, 2018.
[57] W. Jiang, F. Nie, and H. Huang, “Robust dictionary learning withcapped l 1-norm,” in Proceedings of the 24th International Conferenceon Artificial Intelligence. AAAI Press, 2015, pp. 3590–3596.
Rudrasis Chakraborty received his Ph.D. in computer science from theUniv. of Florida in 2018. He is currently a post doctoral researcher at
UC Berkeley. His research interest lies in the intersection of Geometry,
Machine Learning and Computer Vision.
Liu Yang is an undergraduate student in Xi’an Jiaotong University. She
is currently an exchange student in UC Berkeley. Her research interest
lies in Machine Learning and Computer Vision.
Søren Hauberg received his Ph.D. in computer science from the Univ. of
Copenhagen in 2012. He has been a visiting scholar at UC Berkeley(2010), and a post doc at the Perceiving Systems department at the Max
Planck Institute for Intelligent Systems in Tübingen, Germany (2012–2014). He is currently an associate professor at the Section for Cognitive
Systems, at the technical Univ. of Denmark. His research is at theinterplay of geometry and statistics.
Baba C. Vemuri is the Wilson and Marie Collins professor of Engineering
at the Department of Computer Information Science and Engineering andthe Department of Statistics, University of Florida. His research interests
lie in Geometric Statistics, Machine Learning, Computer Vision andMedical Imaging. He is an associte editor of the IJCV and MedIA journals.
He received the IEEE Computer Society’s Technical Achievement Award(2017) and is a Fellow of the ACM (2009) and the IEEE (2001).