Towards a Theoretical Analysis of PCA for Heteroscedastic Data

arXiv:1610.03595v1 [math.ST] 12 Oct 2016

David Hong, Laura Balzano, and Jeffrey A. Fessler

Abstract— Principal Component Analysis (PCA) is a method
for estimating a subspace given noisy samples. It is useful in
a variety of problems ranging from dimensionality reduction
to anomaly detection and the visualization of high dimensional
data. PCA performs well in the presence of moderate noise and
even with missing data, but is also sensitive to outliers. PCA is
also known to have a phase transition when noise is independent
and identically distributed; recovery of the subspace sharply
declines at a threshold noise variance. Effective use of PCA
requires a rigorous understanding of these behaviors. This
paper provides a step towards an analysis of PCA for samples
with heteroscedastic noise, that is, samples that have nonuniform noise variances and so are no longer identically
distributed. In particular, we provide a simple asymptotic
prediction of the recovery of a one-dimensional subspace from
noisy heteroscedastic samples. The prediction enables: a) easy
and efficient calculation of the asymptotic performance, and b)
qualitative reasoning to understand how PCA is impacted by
heteroscedasticity (such as outliers).

I. I NTRODUCTION
Given noisy measurements of points from a subspace, one
may estimate the subspace with Principal Component Analysis (PCA). Estimating a k-dimensional subspace from noisy
samples y1 , . . . , yn ∈ Rd by PCA is accomplished by solving
the non-convex problem
Û =

argmin
U ∈Rd×k :U T U =I

min

zi ∈Rk

n
X
i=1

kyi − U zi k22 ,

(1)

which can be done efficiently via the singular value decomposition. PCA performs well in the presence of low
to moderate noise and even performs well with missing
data [1], [2]. Furthermore, for mean zero data, representing
the samples in the basis produced by PCA gives coordinates
that are uncorrelated and provide a convenient representation
of the data where the relevant factors have been decoupled.
As a result of such nice properties, PCA has been applied
in myriad contexts to accomplish tasks such as dimensionality reduction, anomaly detection and the visualization of high
dimensional data. A small sample of these settings include
medical imaging [3], anomaly detection on computer networks [4] and dimensionality reduction for classification [5].
It has also been used to model images taken of a scene under
various illuminations [6] as well as measurements taken in
environmental monitoring [7], [8], to name just a few.
To use PCA effectively in all these settings, it is important
to rigorously understand its performance under a variety of
conditions. It is known, for example, that PCA is sensitive
Work by D. Hong was supported by the National Science Foundation
Graduate Research Fellowship under DGE #1256260. Work by L. Balzano
was supported by the ARO Grant W911NF-14-1-0634. Work by J. Fessler
was supported by the UM-SJTU data science seed fund.

to outliers (i.e., gross errors) [9]. Thus for problems such
as computer vision modeling [10] or foreground-background
separation [11] where outliers may be expected (or even
of interest), robust variants [1] are used instead. PCA with
independent identically distributed noise is also known to
exhibit a phase transition; recovery of the subspace sharply
declines after the noise variance exceeds a threshold [12].
This paper provides a step towards extending such analysis
to the case where noise is heteroscedastic, that is, the case
where samples have non-uniform noise variances and so are
no longer identically distributed. In particular, we provide
a simple asymptotic prediction of the recovery of a onedimensional subspace from noisy heteroscedastic samples.
Forming the prediction involves connecting several results
from random matrix theory to obtain an initial complicated
asymptotic prediction and then exploiting its structure to find
a much simpler algebraic description.
The simple form enables: a) easy and efficient calculation
of the asymptotic prediction, and b) reasoning qualitatively
about the expressions to understand the asymptotic behavior
of PCA with heteroscedastic noise. We demonstrate these
benefits through an example calculation and a qualitative
analysis that explains a surprising phenomenon: the largest
noise variance seems to most heavily influence performance.
We also perform numerical experiments to illustrate how the
asymptotic prediction applies for particular (finite) choices
of ambient dimension and number of samples.
The rest of the paper is organized as follows. Section II
describes the model we consider (a one-dimensional signal in
heteroscedastic noise) and states the main result: an asymptotic prediction for the recovery of the one-dimensional
subspace by PCA. It also includes an example calculation
of the asymptotic prediction for a particular set of model
parameters, illustrating how the main result enables easy and
efficient calculation of the prediction. Section III compares
the prediction with experimental results simulated according
to the model. The simulations demonstrate good agreement
as the ambient dimension and number of samples grow large;
when these values are small the prediction and experiment
differ but have the same behavior. Section IV provides a
proof of the main result. Section V uses the main result to
provide a qualitative analysis of the behavior of PCA under
heteroscedastic noise, revealing some interesting phenomena about the negative impact of heteroscedasticity. Finally,
Section VI discusses the findings and describes avenues for
future work.

We model n heteroscedastic samples y1 , . . . , yn ∈ R from
a one dimensional subspace ũ ∈ Rd as
yi = θũzi + ηi εi

(2)

where
• θ ∈ R+ is the subspace amplitude,
• ηi ∈ R+ are the noise standard deviations,
iid
d
• ũ ∈ R is the subspace and has entries ũj ∼ F1 (0, 1/d)
with mean zero and variance 1/d,
iid
• zi ∼ F2 (0, 1) are random subspace coefficients and
have mean zero and unit variance, and
d
• εi ∈ R are independent noise vectors that have entries
iid
εij ∼ F3 (0, 1) with mean zero and unit variance,
such that the distributions F1 and F2 satisfy the log-Sobolev
inequality [13] and the distribution F3 satisfies condition
(1.3) from [14]. Notably, these conditions are satisfied by
Gaussian distributions F1 = F2 = F3 = N .
We further suppose that L noise levels σ1 , . . . , σL occur
in proportions p1 , . . . , pL . Namely, p1 of the samples have
noise level ηi = σ1 , p2 have ηi = σ2 and so on, where the
p` values sum to unity.
The following theorem is our main result and describes
how well the subspace ũ is recovered by PCA as the problem
dimensions grow.
Theorem 1: For fixed samples-to-dimension ratio c > 1,
the PCA estimate û is such that


A (β)
2
a.s.
T
(3)
ũ û
−→ max 0,
n, d → ∞
βB 0 (β)

1.2

1.2

1.0

1.0

0.8

0.8

|ũT û|2

d

|ũT û|2

II. M AIN RESULT

0.6

0.6

0.4

0.4

0.2
0

0.2
0

(a) Results for d = 102 and
n = 103 (10000 trials).

B (x) = 1 −

A(β)
βB 0 (β) ,

Example calculation: Here we calculate the asymptotic
prediction in (3) for the case:
c=5

p = (0.2, 0.8)

θ=2

σ = (1, 2)
2

Namely, we determine the limit of ũT û when d, n → ∞
for the case where there are 5 times as many samples as the
ambient dimension, the signal amplitude is 2, and 20% of
the samples have low noise with variance 1 (signal-to-noise
ratio θ2 /σ12 = 4) and 80% of the samples have high noise
with variance 4 (signal-to-noise ratio θ2 /σ22 = 1). The steps
are as follows.

(b) Results for d = 103 and
n = 104 (1000 trials).

obtaining
A (β)
= 0.705.
βB 0 (β)

4) Take the maximum with zero, and conclude that
ũT û

2

a.s.

−→

0.705.

n, d → ∞
n/d = c

and β is the largest real root of B.
Section IV presents the proof of this theorem. We illustrate
Theorem 1 with the following example calculation.

1

β = 23.466.
3) Evaluate

A (x) = 1 −

1/2
p2

1) Substitute the values of c, θ, p and σ into the formulas
for A and B, obtaining
!
0.8 · 24
0.2 · 14
A (x) = 1 − 5 ·
2 +
2
(x − 12 )
(x − 22 )
64
1
=1−
2 −
2
(x − 1)
(x − 4)


0.2
0.8
B (x) = 1 − 5 · 22 ·
+
x − 12
x − 22
4
16
=1−
−
.
x−1 x−4
2) Find the largest root of B, obtaining

L
X

p` σ`4
c
2 2
`=1 (x − σ` )
L
X
p`
cθ2
x − σ`2
`=1

0

Fig. 1: Simulation results for c = 10, θ = 1, σ = (1.8, 0.2)
where p2 is swept from 0 to 1 with p1 = 1 − p2 . Simulation mean
(blue curve) and interquartile interval (light blue ribbon) shown with
asymptotic prediction (red curve).

n/d = c

where

0

1

1/2
p2

Note that the second step can be easily done by clearing the
denominator of B and finding the real roots of the resulting
degree L polynomial (using off-the-shelf tools). Hence the
asymptotic prediction can be efficiently computed.
III. E XPERIMENTAL VERIFICATION
To illustrate the main result (Theorem 1) we performed a
numerical experiment for the two noise level case (L = 2):
c = 10

θ=1

σ = (1.8, 0.2)

where p2 is swept from 0 to 1 with p1 = 1 − p2 . This allows
us to investigate the accuracy of the asymptotic prediction
for a variety of settings: at the extremes (p2 = 0, 1) the
setup matches the homoscedastic setting and in the middle
(p2 = 1/2) the samples are split evenly between the two
noise levels.

We first suppose d = 102 and n = 103 . We performed
10000 trials with data generated by the Gaussian distribution.
Namely, F1 = F2 = F3 = N . Figure 1a shows the
simulation results with the mean (blue curve) and interquartile interval (light blue ribbon) shown with the asymptotic
prediction (red curve).
There is generally good agreement between the mean and
the asymptotic prediction for p2 > 0.6 (they deviate from
each other by no more than 0.025). However, the smaller the
value of p2 , the greater the deviation of the prediction from
the mean, with the asymptotic prediction underestimating the
(non-asymptotic) simulation by at most 0.182.
Figure 1 illustrates a general phenomenon we observed:
small asymptotic predictions typically underestimate the
simulation results, with smaller predictions underestimating
more. Intuitively, a small prediction of this inner product corresponds to a subspace estimate that is becoming increasingly
like an isotropically random vector and so it has vanishing
square inner product with the true subspace as the dimension
grows. Asymptotically, below the phase transition, theory
predicts the subspace estimate to be practically isotropically
random. However, in finite dimension there is a better chance
of alignment, resulting in a positive square inner product.
We illustrate this phenomenon with a second experiment
with a higher dimension of d = 103 and a higher number
of samples n = 104 (chosen to have the same sample-todimension ratio c = 10). The data were again generated
using the Gaussian distribution and we performed 1000 trials.
Figure 1b shows the simulation results for this case. Again,
the mean and interquartile interval are in blue and light blue,
and the asymptotic prediction is in red. This experiment
demonstrates better agreement between the mean behavior
and the asymptotic prediction. In particular, for p2 > 0.4
they deviate from each other by no more than 0.017. For
p2 < 0.4 their largest deviation is 0.084; notably, this is less
than half of the largest deviation for the smaller experiment.
Furthermore, the interquartile interval is narrower, indicating
that the square inner product is concentrating.
We stress here that, while we have not proved this relationship, experimental evidence suggests that the asymptotic prediction in Theorem 1 underestimates the mean square inner
product obtained in finite size experiments, and is therefore
a conservative or pessimistic estimate of the performance of
PCA. When determining how much to trust the subspace
estimated by PCA, this conservatism might be preferable to
overestimating the reliability of PCA.
Finally, note that the square inner product from the simulation is sometimes above one. This is because the true ũ is a
random vector and may not have unit norm. However, as the
dimension grows, the norm of ũ will concentrate around one
(see [15] for a treatment of the concentration of the norm).
IV. P ROOF OF MAIN RESULT
This section proves the main result (Theorem 1). The proof
has eight main parts. In IV-A, we first apply previous
results in the literature to obtain an initial expression for the
asymptotic prediction. This expression is difficult to evaluate

and analyze because it involves an integral transform of the
(nontrivial) limiting singular value distribution for a random
(noise) matrix as well as the corresponding limiting largest
singular value. In the remaining seven parts (IV-B-IV-H),
we find a simple equivalent expression by exploiting the
structure of the prediction.
A. Obtain an initial expression
Rewriting the model in (2) in matrix form yields
Y := θũz T + EH ∈ Rd×n
where E ∈ Rd×n is a matrix with columns ε1 , . . . , εn ∈ Rd
and H ∈ Rn×n is a diagonal matrix with diagonal entries
η1 , . . . , ηn ∈ R.
Recall that the subspace basis û estimated by PCA for
Y is the first left singular vector of Y. PCA is invariant to
scaling, so û is also the first left singular vector of
1
Ỹ := √ Y.
n
The matrix Ỹ matches the low rank (here rank one) perturbation of a random matrix model considered in [12] because
Ỹ = P + X
where

P := θũ

1
√ z
n

T


X :=


1
√ E H,
n

and P is generated according to the “i.i.d model” and
satisfies Assumption 2.4 of [12], and X satisfies Assumptions
2.1-2.3 of [12] (X here matches the random matrix in [14],
which [12] refers to as an example of a random matrix that
satisfies the assumptions).
Thus under the condition ϕ0 (b+ ) = −∞ (we will show
in subsection IV-C that it is indeed satisfied), Theorems 2.10
and 2.11 from [12] yield
(
−2ϕ(ρ)
2
θ2 > θ̄2
a.s.
T
θ 2 D 0 (ρ)
ũ û
−→
(4)
n, d → ∞
0
otherwise
n/d = c
where

−1
• ρ := D
1/θ2
2
+
• θ̄ := 1/D (b )

−1
−1
• D (z) := ϕ (z) c
ϕ (z) + 1−cz
z>b
Rb z
• ϕ (z) :=
dµX (t)
z>b
a z 2 −t2
• a and b are, respectively, the infimum and the supremum
of the support of µX (so b > a ≥ 0), and
• µX is the limiting singular value distribution of X
(compactly supported by Assumption 2.1 of [12]).
We use the notation f (b+ ) := limz→b+ f (z) as a convenient
shorthand for the limit of a function f (z).
Evaluating this asymptotic prediction would then consist
of evaluating the above intermediates from bottom to top.
These steps are challenging because because they involve an
integral transform of the limiting singular value distribution
for the random (noise) matrix as well as the corresponding

limiting largest singular value. The following sections eventually lead to the simpler expression in (3) that is easier to
both evaluate and analyze.
B. Carry out a change of variables

For all z and ξ > 0, z 2 c + iξ ∈ C+ and so combining (6)-(8)
yields that for all z > b


 −1
ψ (z) = − lim m z 2 c + iξ
ξ→0+

We begin by introducing the function
" Z
#−1
1
cz
1 b
dµX (t)
, z > b (5)
ψ (z) :=
=
ϕ (z)
c a z 2 − t2
because it turns out to have several nice properties that
simplify all of the following analysis.
Rewriting (4) using ψ (z) instead of ϕ (z) yields

−2c
θ2 > θ̄2
2
a.s.
θ 2 ψ(ρ)D 0 (ρ)/ρ
ũT û
−→
n, d → ∞
0
otherwise

= lim+ z 2 c + iξ − c
ξ→0

= z2c − c
= z2c − c

L
X
`=1

L
X
`=1
L
X
`=1

1+

σ`2

p` σ`2
1 + σ`2 m (z 2 c + iξ)

p` σ`2
limξ→0+ m (z 2 c + iξ)

p` σ`2
.
1 − σ`2 /ψ (z)

Rearranging yields

n/d = c

where now
D (z) =

cz 2
2

(ψ (z))

+

c−1
,
ψ (z)

z > b.

∀z > b, 0 =

Property 1. We show that ψ (z) satisfies a certain rational
equation for all z > b. For this, observe that the square
singular values of the noise matrix X are exactly the eigenvalues of cXXT divided by c (since X has more columns
than rows, namely, c = n/d > 1). Thus we first consider the
limiting eigenvalue distribution µcXXT of cXXT , and then
relate its Stieltjes transform m (ζ) to ψ (z).
Theorem 1 in [14] establishes that the random matrix



T
1
1
T
2
cXX = √ E H √ E
d
d
has a limiting eigenvalue distribution µcXXT whose Stieltjes
transform
Z
1
T (t) ,
m (ζ) :=
dµ
ζ ∈ C+
(6)
t − ζ cXX
satisfies the condition
!−1
L
X
p` σ`2
+
∀ζ ∈ C
m (ζ) = − ζ − c
(7)
1 + σ`2 m (ζ)
`=1

+

where C is the set of all complex numbers with positive
imaginary part.
Since the square singular values of X are exactly the
eigenvalues of cXXT divided by c, we have for all z > b
" Z
#−1
1 b
1
ψ (z) =
dµX (t)
c a z 2 − t2
" Z 2
#−1
1 bc
1
T (t)
=
dµ
c a2 c z 2 − t/c cXX
"Z 2
#−1
b c
1
=−
dµcXXT (t)
.
(8)
2
a2 c t − z c

1
c X
p` σ`2
. (9)
2 − ψ (z) − ψ (z)
ψ (z) − σ`2
(ψ (z))
`=1

where the last term is

C. Find useful properties of ψ (z)
Establishing some properties of ψ (z) aids simplification
significantly. Furthermore, these properties help us show that
ϕ0 (b+ ) is indeed −∞, as stated above in subsection IV-A.

L

cz 2

−



L
L
c X
p` σ`2
c X
ψ (z)
p
=
1
−
`
ψ (z)
ψ (z) − σ`2
ψ (z)
ψ (z) − σ`2
`=1

`=1

L

=

X
c
p`
−c
ψ (z)
ψ (z) − σ`2
`=1

because p1 + · · · + pL = 1. Substituting back into (9) finally
yields 0 = Q (ψ (z) , z) for all z > b, where
L

Q (s, z) :=

X p`
c−1
cz 2
+
−c
.
2
s
s
s − σ`2

(10)

`=1

Thus ψ is an algebraic function with associated rational
function Q (a polynomial can be formed by clearing the
denominator).
Property 2. We show that ψ (b+ ) is finite and ψ 0 (b+ ) = ∞.
For this, note first that ψ (b+ ) is a multiple root of Q (·, b)
and hence is finite. This follows from the observation in [16]
that non-pole boundary points of compactly supported distributions like µcXXT occur where the polynomial defining
the Stieltjes transform has multiple roots.
Differentiating 0 = Q (ψ (z) , z) with respect to z and
rearranging yields
∂Q
∂z
ψ 0 (z) = − ∂Q
∂s

(ψ (z) , z)
(ψ (z) , z)

.

Since ψ (b+ ) is a multiple root of Q (·, b),
 
∂Q
ψ b+ , b = 0
∂s
while on the other hand
 
2cb
∂Q
ψ b+ , b =
2 > 0.
∂z
(ψ (b+ ))
Thus ψ 0 (b+ ) = ∞, where the sign is necessarily positive
because ψ (z) is an increasing function.

Summarizing, we have shown that
1) ψ satisfies the equation 0 = Q (ψ (z) , z) for all z > b
2) ψ (b+ ) is finite and ψ 0 (b+ ) = ∞
As an immediate consequence of these properties, we also
have that indeed



ψ 0 (b+ )
c
1
−
z
= −∞.
ϕ0 b+ =
ψ (b+ )
ψ (b+ )

where
A (x) := 1 − c
Thus
ψ 0 (z) =

D (z) = Q (ψ (z) , z) + c

`=1

D (z) = c

`=1

0

D (z) = −cψ (z)

`=1

0

Thus the condition θ2 > θ̄2 is equivalent to
0>1−

L
X

θ2
p`
2
= B ψ b+
=
1
−
cθ
2
+
2
ψ (b ) − σ`
θ̄
`=1

where
B (x) := 1 − cθ2

L

X 2p` ψ (z)
2cz 2
= −2 (c − 1) + c
ψ (z)
ψ (z) − σ`2
`=1

because 0 = Q (ψ (z) , z) for z > b. Substituting into γ(z)
and forming a common denominator yields
L
L
2
X
X
2p` ψ (z)
p` (ψ (z))
γ (z) = 1 − c + c
−
c
2
ψ (z) − σ`2
(ψ (z) − σ 2 )
`

2

p`

`=1

(ψ (z)) − 2ψ (z) σ`2
2

(ψ (z) − σ`2 )

Dividing the summand with respect to ψ(z) and recalling
that p1 + · · · + pL = 1 yields
!
L
X
p` σ`4
γ(z) = 1 − c + c
p` −
2
(ψ (z) − σ`2 )
`=1
=1−c

L
X
`=1

p` σ`4
(ψ (z) −

2
σ`2 )

= A (ψ (z))

L
X
`=1

p`
.
x − σ`2

Using (13) yields
r :=

=

−2c

θ2 ψ (ρ) D0

(ρ) /ρ

=

A (ψ (ρ))
PL
p`
`=1 (ψ(ρ)−σ 2 )2
`

ψ (ρ) cθ2

A (ψ (ρ))
ψ (ρ) B 0 (ψ (ρ))

where we note that
B 0 (x) = cθ2

Note that

=1−c+c

(13)

`

L
X

p`
1
+
.
=
D
b
=
c
+
2
ψ (b ) − σ`2
θ̄
`=1

2
σ`2 )

`=1

(12)

`=1

(11)

and so we need to express ψ (z) in terms of ψ (z).
To do this, differentiate both sides of 0 = Q (ψ (z) , z)
with respect to z and solve for ψ 0 (z), obtaining
2cz
ψ 0 (z) =
γ (z)
where the denominator is
L
2
X
2cz 2
p` (ψ (z))
γ (z) := c − 1 +
−c
.
2 2
ψ (z)
`=1 (ψ (z) − σ` )

`=1
L
X

2cz
A (ψ (z))
L

p`
(ψ (z) −

2.

This subsection uses (11) and (13) to express the asymptotic
prediction in terms of ψ (b+ ) and ψ (ρ). Using (11) yields

E. Express D0 (z) /z in terms of only ψ (z)
This subsection uses the properties of ψ (z) to find a simple
expression for D0 (z) /z in terms of ψ (z). Differentiating (11) with respect to z yields
0

(x − σ`2 )

F. Express the prediction in terms of only ψ (b+ ) and ψ (ρ)

p`
.
ψ (z) − σ`2

p`
.
ψ (z) − σ`2

L
X

`=1

X
D0 (z)
2c2
p`
=−
2.
z
A (ψ (z))
(ψ (z) − σ 2 )

Recalling that 0 = Q (ψ (z) , z) for z > b, we have
L
X

p` σ`4

and so

D. Express D (z) in terms of only ψ (z)
This subsection uses the properties of ψ (z) to find a simple
expression for D (z) in terms of ψ (z). Observe that
L
X

L
X

L
X

p`

`=1

(x − σ`2 )

2.

Summarizing, the asymptotic prediction is now expressed as
(
A(ψ(ρ))
2
B (ψ (b+ )) < 0
a.s.
T
ψ(ρ)B 0 (ψ(ρ))
ũ û
−→
n, d → ∞
0
otherwise.
n/d = c
G. Express the prediction algebraically
This subsection finds an algebraic description of the asymptotic prediction. We first use the properties of ψ (z) to show
that ψ (b+ ) and ψ (ρ) are, respectively, the largest real roots
of A and B.
Using (12) yields

2cz
A ψ b+ = 0 + = 0
ψ (b )
because ψ 0 (b+ ) = ∞. Thus ψ (b+ ) is a real root of A.
For θ > θ̄, we have ρ := D−1 1/θ2 and so
0 = 1 − θ2 D (ρ) = 1 − cθ2

L
X
`=1

p`
= B (ψ (ρ)) .
ψ (ρ) − σ`2

PL

where α and β are, respectively, the largest real roots of A
and B.

p`
`=1 x−σ`2

σ32

σ22

σ12

β

1
cθ 2

x

Fig. 2: Illustration of the real roots of B(x) for L = 3 levels. They
occur where the sum (blue curve) intersects 1/(cθ2 ) (red line). The
largest is β.

Thus ψ (ρ) is a real root of B.
The functions A and B both have several real roots. To
show that ψ (b+ ) and ψ (ρ) are the largest ones, consider

ψ (ρ) = ψ D−1 1/θ2
as a function of θ as θ increases from θ̄ to infinity. Note that
" Z
#−1
1 b
1
ψ (z) =
dµX (t)
, z>b
c a z 2 − t2
continuously and monotonically increases from ψ (b+ ) to
infinity as z increases from b to infinity. Thus
D (z) = c

L
X
`=1

p`
,
ψ (z) − σ`2

z>b

continuously and monotonically decreases from 1/θ̄2 to
zero as z increases from b to infinity, and so D−1 1/θ2
continuously and monotonically increases from b to infinity
as θ increases from θ̄ to infinity.
As a result, ψ (ρ) continuously and monotonically increases from ψ (b+ ) towards infinity as θ increases from
θ̄ towards infinity. This is possible only if both ψ (b+ ) and
ψ (ρ) are larger than all the noise levels σ`2 . To see this, recall
that ψ (ρ) is a real root of B and note that the real roots of
B satisfy the following equation (illustrated in Figure 2):
L

X p`
1
=
.
2
cθ
x − σ`2
`=1

If either ψ (b+ ) or ψ (ρ) were less than any of the noise
levels σ`2 , then ψ(ρ) would change discontinuously as θ
varies. Thus ψ (b+ ) and ψ (ρ) are indeed both larger than
all the noise levels.
To the right of all the noise levels (i.e., for x larger than
all the noise levels σ`2 ), both A and B continuously and
monotonically increase from negative infinity to one, and so
each has exactly one real root larger than all the noise levels
(namely, the largest real root). Thus ψ (b+ ) is the largest real
root of A and when θ > θ̄, ψ (ρ) is the largest real root of B.
Using this yields the algebraic form of the prediction:
(
A(β)
2
B (α) < 0
a.s.
βB 0 (β)
ũT û
−→
n, d → ∞
0
otherwise
n/d = c

H. Further simplify the asymptotic prediction
We further simplify the asymptotic prediction by showing
that B (α) < 0 is equivalent to A (β) / (βB 0 (β)) > 0.
To do this, observe that both α and β are larger than all
the noise levels σ`2 , and note that A (x) and B (x) are both
monotonically increasing in this regime. Thus it follows that
B (α) < 0 ⇐⇒ α < β ⇐⇒ 0 < A (β)
because B (β) = 0 and A (α) = 0.
Furthermore B 0 (β) > 0 (since B is increasing in this
regime) and β > 0. Thus
A (β) > 0 ⇐⇒

A (β)
>0
βB 0 (β)

Using this equivalence finally leads to the main result in (3).
V. Q UALITATIVE ANALYSIS
This section applies the main result (Theorem 1) in several
settings to gain some insights into the performance of PCA.
A. Dependence on balance of noise variances
Here we would like to understand how the balance of the
noise variances affects the performance. We consider a case
with two noise levels where we sweep the noise variances
σ12 and σ22 while holding fixed the average noise variance
σ̄ 2 = p1 σ12 + p2 σ22 .
In particular, consider
c = 10

p = (0.7, 0.3)

θ=1

σ̄ = 1.3

where we sweep over λ ∈ [0, 1] and set

λ
σ̄ 2
p1 λ + p2 (1 − λ)
(1 − λ)
σ22 =
σ̄ 2 .
p1 λ + p2 (1 − λ)
σ12 =

As intended, this fixes the average noise variance at
n

1X 2
η = p1 σ12 + p2 σ22 = σ̄ 2 .
n i=1 i
Sweeping over λ adjusts the breakdown of the average noise
variance σ̄ 2 across the two noise levels. It is not initially
obvious whether better performance will occur halfway when
λ = 1/2 or at the extremes when λ = 0 or λ = 1.
When λ = 1/2 both noise levels are the same and so all
the samples have noise variance σ̄ 2 and this reduces to the
previously considered homoscedastic case analyzed as an
example in [12]. When λ = 0 or λ = 1, some of the points
have no noise (and so PCA may do better), but the rest have
noise larger than σ̄ 2 (and so PCA may do worse).
Figure 3a shows that the asymptotic prediction has a peak
at λ = 1/2; recovery is best when the two noise levels
are the same. In other words, having samples with more

σ22

p1
1
= σ̄ 2 − σ12 ,
p2
p2

which has slope −p1 /p2 = −7/3. Along this line, the
prediction does indeed decrease away from the diagonal
σ12 = σ22 .
Figure 3b illustrates that the prediction seems to depend
primarily on the larger of the two noise variances. This is
initially surprising but might be understood by considering
the root β of B. Recall that it is the largest value x satisfying

σ22

0.4
0.2
0

θ=1

p1 = 1

σ12

= σ̄

2

and we sweep over c > 1 and σ̄ 2 > 0, as shown in Figure 5a.
Note that this is the homoscedastic case analyzed as an
example in [12] and Figure 5a illustrates the predicted phase
transition at c = σ̄ 4 . In fact, the asymptotic prediction in (3)
specializes to the prediction in [12] for the case where there
is only one noise level and hence the noise is homoscedastic.

0.8

3

0.6

2

0.4

1

0.2
1

(a) Sweeping over noise levels
while keeping the average noise
variance fixed at σ̄ 2 = 1.69.
When λ = 1, σ12 is large and
σ22 = 0. When λ = 0, σ12 = 0
and σ22 is even larger.

2
σ12

3

0

4

(b) Sweeping over both noise
levels independently. The solid
black curves are contours. On
the dotted cyan line, the average
noise variance is σ̄ 2 = 1.69.

Fig. 3: Asymptotic prediction under various noise levels for c = 10,
p = (0.7, 0.3), θ = 1. Namely, 70% of samples have noise variance
σ12 and 30% have noise variance σ22 .

PL

p`
`=1 x−σ`2

`=1

Now consider holding everything constant except for the
sample-to-dimension ratio and average noise variance. We
first suppose that there is only one noise level (alternatively,
two noise levels that are equal). In particular we consider

4

1

1/2
λ

L

B. Dependence on sample-to-dimension ratio and average
noise variance

1

0

X p`
1
=
2
cθ
x − σ`2
as illustrated in Figure 4. This figure suggests that the largest
root is heavily influenced by the largest noise variance (it
is the nearest pole). As a result, changing the other noise
variances has much less impact on β and on the prediction.
The precise relative impact does depend on the proportions,
as seen in Figure 3b, where the shape of the level curves
are not symmetric around σ12 = σ22 (if the performance
depended exclusively on the max of σ12 and σ22 , it would
necessarily be symmetric). Nevertheless, for any proportion,
large noise variances can drown out the influence of small
noise variances.
This insight gives a rough explanation of why it may be
generally preferable to have equal noise variances (λ = 1/2)
for some fixed average noise variance. Imbalance (λ 6= 1/2)
means that one of the noise variances will be larger and
cause the performance to decline even though the other noise
variance is smaller.

5

|ũT û|2

0.6
|ũT û|2

noise hurts more than having samples with corresponding
less noise helps, regardless of which set is larger. This seems
to be a general phenomenon; the same has occurred for other
choices of parameters we tried.
To further investigate, we use the same parameters c, p, θ
as before but sweep over both σ12 and σ22 (independently) and
produce a heatmap of the asymptotic prediction, shown in
Figure 3b. On this figure, adjusting λ corresponds to moving
along the line (shown as light blue dashes):

σ22

1
cθ 2

x

σ12 β

Fig. 4: Location of the largest root β of B(x) for λ = 0.7 (i.e.,
σ12 = 2.04 and σ22 = 0.874), c = 10, p = (0.7, 0.3) and θ = 1.

We now consider an analogous setting where the noise is
imbalanced (i.e., heteroscedastic) with
1 2
p1 = 0.9
σ12 =
σ̄
2p1
1 2
p2 = 0.1
σ22 =
σ̄
2p2
so that

1 2 1 2
σ̄ + σ̄ = σ̄ 2 .
2
2
Figure 5b illustrates a similar behavior with a phase transition
further left. Namely, more samples are needed than in the
homoscedastic setting for the same average noise variance.
This agrees with the previous observation that performance
for a given average noise variance is best when all the points
have the same noise variance (i.e., are homoscedastic).
p1 σ12 + p2 σ22 =

C. Dependence on sample proportions
Finally we revisit the sweep carried out in the numerical
experiments of Section III. Recall that everything but the
proportions were fixed. In particular,
c = 10

θ=1

σ = (1.8, 0.2)

and p2 varied from 0 to 1 with p1 = 1 − p2 . Figure 1 shows
the prediction as a red curve (identical in both sub-figures).
As expected, the performance is best when p2 = 1 and all
the samples have the lower noise variance; it is preferable to

1

10

0.4

0.6

0.2

5

1

0

1

1

2

3
σ̄

4

5

0.8

10

5
0

1

15
c

0.6

|ũT û|2

15
c

20

0.8

2

(a) Homoscedastic (i.e., identically distributed) noise.

0.4

|ũT û|2

20

0.2
0

1

2

3
σ̄

4

5

0

2

(b) Heteroscedastic (i.e., imbalanced) noise.

Fig. 5: Asymptotic prediction as a function of average noise
variance σ̄ 2 and sample-to-dimension ratio c. Contours are overlaid
in black. Note that the phase transition in (b) is further left than in
(a); more samples are needed to tolerate the same amount of noise.

have a larger proportion of low noise samples. Interestingly,
the benefit of having more low noise samples is not uniform
through the range. The slope for small p2 (close to zero) is
less steep than that for high p2 (close to one). Hence, having
a larger proportion of low noise samples is not as helpful
when there are only a few low noise samples otherwise.
A careful investigation of this phenomenon would be an
interesting area of further study.
VI. D ISCUSSION AND EXTENSIONS
This paper considered PCA when noise is heteroscedastic
and provided a step towards the analysis of the recovery of
the subspace. In particular, we provided a simple asymptotic
prediction for the recovery of a one-dimensional subspace
by PCA from noisy heteroscedastic samples. We provided an
example, illustrating how the simple form enables easy and
efficient calculation of the asymptotic prediction, as well as
an experimental verification of the prediction in simulation.
Next, we used the simple form to reason qualitatively about
the asymptotic prediction and gain new insights about the
performance of PCA. Namely, we found that the performance
seems to often be most heavily influenced by the largest noise
variance present in the data. Hence, heteroscedasticity tends
to have a negative impact on the performance of PCA.
There are many avenues for potential extensions and further work. A natural direction is to extend this work to multidimensional subspaces. Another avenue of future work will
be to consider a weighted version of PCA, where the samples
are first weighted in the objective function (1) to reduce the
impact of very noisy points. Unfortunately, applying these
weights violates the construction of the subspace coefficients

as being identically distributed and so this is a challenging
extension. Other avenues include further investigation of the
phenomena discussed in the qualitative studies above as well
as further study of the algebraic structure of the expressions
in the prediction.
R EFERENCES
[1] J. Wright, A. Ganesh, S. Rao, Y. Peng, and Y. Ma, “Robust principal
component analysis: Exact recovery of corrupted low-rank matrices via
convex optimization,” in Advances in Neural Information Processing
Systems 22, 2009, pp. 2080–2088.
[2] S. Chatterjee, “Matrix estimation by universal singular value thresholding,” The Annals of Statistics, vol. 43, no. 1, pp. 177–214, 02 2015.
[3] B. A. Ardekani, J. Kershaw, K. Kashikura, and I. Kanno, “Activation
detection in functional MRI using subspace modeling and maximum likelihood estimation,” IEEE Transactions on Medical Imaging,
vol. 18, no. 2, pp. 101–114, 1999.
[4] A. Lakhina, M. Crovella, and C. Diot, “Diagnosing network-wide
traffic anomalies,” in Proceedings of the 2004 Conference on Applications, Technologies, Architectures, and Protocols for Computer
Communications, ser. SIGCOMM ’04, 2004, pp. 219–230.
[5] N. Sharma and K. Saroha, “A novel dimensionality reduction method
for cancer dataset using PCA and feature ranking,” in Advances in
Computing, Communications and Informatics (ICACCI), 2015 International Conference on, 2015, pp. 2261–2264.
[6] R. Basri and D. W. Jacobs, “Lambertian reflectance and linear
subspaces,” IEEE Transactions on Pattern Analysis and Machine
Intelligence, vol. 25, no. 2, pp. 218–233, 2003.
[7] S. Papadimitriou, J. Sun, and C. Faloutsos, “Streaming pattern discovery in multiple time-series,” in Proceedings of the 31st International
Conference on Very Large Data Bases, ser. VLDB ’05, 2005, pp. 697–
708.
[8] G. S. Wagner and T. J. Owens, “Signal detection using multi-channel
seismic data,” Bulletin of the Seismological Society of America, vol. 86,
pp. 221–231, 1996.
[9] I. Jolliffe, Principal Component Analysis. Springer-Verlag, 1986.
[10] F. D. la Torre and M. J. Black, “Robust principal component analysis
for computer vision,” in Computer Vision, 2001. ICCV 2001. Proceedings. Eighth IEEE International Conference on, vol. 1, 2001, pp.
362–369 vol.1.
[11] J. He, L. Balzano, and A. Szlam, “Incremental gradient on the
grassmannian for online foreground and background separation in
subsampled video,” in Computer Vision and Pattern Recognition
(CVPR), 2012 IEEE Conference on, 2012, pp. 1568–1575.
[12] F. Benaych-Georges and R. R. Nadakuditi, “The singular values
and vectors of low rank perturbations of large rectangular random
matrices,” Journal of Multivariate Analysis, vol. 111, pp. 120 – 135,
2012.
[13] G. Anderson, A. Guionnet, and O. Zeitouni, An introduction to random
matrices. Cambridge university press, 2009, vol. 118.
[14] G. Pan, “Strong convergence of the empirical distribution of eigenvalues of sample covariance matrices with a perturbation matrix,” Journal
of Multivariate Analysis, vol. 101, no. 6, pp. 1330 – 1338, 2010.
[15] T. Vincent, L. Tenorio, and M. Wakin, “Concentration of measure:
fundamentals and tools,” Lecture Notes. [Online]. Available:
http://www.stat.rice.edu/∼jrojo/PASI/lectures/TyronCMarticle.pdf
[16] R. R. Nadakuditi and A. Edelman, “The polynomial method for
random matrices,” Foundations of Computational Mathematics, vol. 8,
no. 6, pp. 649–702, 2008.

