Ch.

7 Low-rank approximation & multidimensional scaling

Contents (class version)
7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.4
7.2 Low-rank approximation via Frobenius norm . . . . . . . . . . . . . . . . . . . . 7.5
7.2.1 Eckart-Young-Mirsky theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.6
7.2.2 Subspace approximation perspective . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.9
7.2.3 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.9
7.2.4 Choosing rank via permutation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.12
7.2.5 Non-uniqueness of SVD and low-rank approximation . . . . . . . . . . . . . . . . . . . 7.15
7.2.6 1D example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.17
7.2.7 Generalization to other norms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.21
[Link] Proof for spectral norm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.22
M ×N
7.2.8 Bases for F . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.23
[Link] Canonical basis for FM ×N . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.23
[Link] Orthonormal bases for FM ×N . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.25

7.2.9 Low-rank approximation summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.26

[Link] Generalizations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.27

7.2.10 Rank and stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.27

7.2.11 Example: Photometric stereo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.29

7.1
© J. Fessler, September 2, 2023, 07:03 (class version) 7.2

7.3 Sensor localization application: multidimensional scaling . . . . . . . . . . . . . . 7.31

7.3.1 Derivation (analysis) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.32
7.3.2 MDS method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.37
7.3.3 Practical implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.39
7.3.4 Extensions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.42
7.4 Proximal operators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.43
7.4.1 Soft thresholding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.44
7.4.2 Hard thresholding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.47
7.5 Alternative low-rank approximation formulations . . . . . . . . . . . . . . . . . . 7.49
7.5.1 Unconstrained / regularized formulation . . . . . . . . . . . . . . . . . . . . . . . . . . 7.49
7.5.2 General unitarily invariant formulations . . . . . . . . . . . . . . . . . . . . . . . . . . 7.50
7.5.3 Singular value hard thresholding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.51
7.5.4 Singular value soft thresholding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.52
7.5.5 Other extensions of low-rank approximation . . . . . . . . . . . . . . . . . . . . . . . . 7.54
7.6 Choosing the rank or regularization parameter . . . . . . . . . . . . . . . . . . . . 7.56
7.6.1 Stein’s unbiased risk estimate (SURE) . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.58
7.6.2 OptShrink . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.61
[Link] Illustration of OptShrink’s robustness to rank estimate . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.64

7.7 Related methods: autoencoders and PCA . . . . . . . . . . . . . . . . . . . . . . . 7.68

© J. Fessler, September 2, 2023, 07:03 (class version) 7.3

7.7.1 Relation to autoencoder with linear layers . . . . . . . . . . . . . . . . . . . . . . . . . 7.68

7.7.2 Relation to principal component analysis (PCA) . . . . . . . . . . . . . . . . . . . . . . 7.71
[Link] Maximum variance projections . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.74
[Link] Units and normalization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.76
[Link] PCA extensions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.77

7.8 Subspace learning for classification . . . . . . . . . . . . . . . . . . . . . . . . . . 7.78

7.8.1 Subspace clustering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.81
7.9 Subspace tracking and streaming PCA . . . . . . . . . . . . . . . . . . . . . . . . 7.83
7.9.1 Incremental SVD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.83
7.9.2 Streaming PCA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.85
7.10 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.86
© J. Fessler, September 2, 2023, 07:03 (class version) 7.4

7.1 Introduction L§8.1

In many applications, dimensionality reduction is important. Uses of dimensionality reduction include visu-
alization, removing noise, and decreasing compute and memory requirements, such as image compression.
( HW ) This chapter focuses on low-rank approximation of a matrix. There are theoretical models for why
big matrices should be approximately low rank [1]. Low-rank approximations are also used to compress large
neural-network models to reduce computation and storage [2].
The chapter begins with the classic approach to approximating a matrix by a low-rank matrix, using a non-
convex formulation that has a remarkably simple SVD solution. §7.3 applies this approach to the source
localization application via the multidimensional scaling (MDS) method. §7.5 turns to convex formulations
of low-rank approximation based on proximal operators from §7.4 that involve singular value shrinkage. §7.6
discusses methods for choosing the rank of the approximation, and describes the optimal shrinkage method
called OptShrink. §7.7 discusses related dimensionality-reduction methods including (linear) autoencoders
and principal component analysis (PCA). §7.8 applies the methods to learning low-dimensionality subspaces
from training data for subspace-based classification problems. Finally, §7.9 extends the method to streaming
applications with time-varying data.
Source material for this chapter includes [3, §8.1].
© J. Fessler, September 2, 2023, 07:03 (class version) 7.5

7.2 Low-rank approximation via Frobenius norm

We are given a matrix A ∈ FM ×N (often large), having rank r ≤ min(M, N ).
To perform dimensionality reduction we want to approximate A by another matrix ÂK having rank K ≤ r.
Often we pick K ≪ r.
To find the “best” ÂK we must define how closely ÂK approximates A. The simplest metric is the Frobenius
norm of the difference. This criterion leads to the follow low-rank approximation problem:
×N
LM B ∈ FM ×N : rank(B) ≤ K .

ÂK ≜ arg min ||B − A|| F = PLM ×N (A), K ≜ (7.1)
M ×N | K
B∈LK {z }
,→ Frobenius norm
×N
This is a non-convex optimization problem because the set LM
K of rank-K matrices is not convex.
Why rank is a non-convex function
M ×N
Example.
 To  is a non-convex function, and LK
 see why rank(·) is a nonconvex
 set, consider:

1 0 0 0 α 0
A= , B= , α ∈ (0, 1) =⇒ rank(αA + (1 − α)B) = rank =2
0 0 0 1 0 1−α
≰ α rank(A) +(1 − α) rank(B) = α · 1 + (1 − α) · 1 = 1.

Remarkably, despite this non-convexity, there is a simple solution for ÂK based on any SVD of A [4].
© J. Fessler, September 2, 2023, 07:03 (class version) 7.6

7.2.1 Eckart-Young-Mirsky theorem

Fact 7.1. (Eckart-Young-Mirsky theorem [4, 5].) (Earlier by Schmidt [6] in infinite-dimensional spaces.)
The best rank (at most) K approximation uses the first (largest) K singular components of A:

K min(M,N )
X X
ÂK = σk uk vk′ = UK ΣK VK′ , where A = U ΣV = ′
σk uk vk′ . (7.2)
k=1 k=1

The approximation error depends on the singular values of the remaining (discarded) terms:
r
Xrank(A)
×N
||ÂK − A|| F = σk2 ≤ ||B − A|| F , ∀B ∈ LM
K . (7.3)
k=K+1

• Clearly the approximation error will be small if the remaining singular values are small compared to the
first K singular values.
• Once again we see that the SVD is key tool. Here, the SVD is used both to construct the approximation
and to quantify the approximation error.
• The original matrix A has M N values. The approximation ÂK is formed from K(M + N + 1) values.
This saves memory if K ≪ min(M, N ). Picture of rectangle and outer products...
© J. Fessler, September 2, 2023, 07:03 (class version) 7.7

Q7.1 What is the rank of K ′

P
k=1 σk uk vk when 1 ≤ K ≤ r = rank(A) ? (Choose best answer.)
A: Always 1 B: Always r C: Usually r D: Always K E: Usually K ??
© J. Fessler, September 2, 2023, 07:03 (class version) 7.8

Proof . (cf. [wiki]) The approximation error expression is simple: (Read)

Xr XK Xr
A − ÂK = σk uk vk′ − σk uk vk′ = σk uk vk′
k=1 k=1 k=K+1
Xr 2 Xr
=⇒ ||ÂK − A|| 2F = σk uk vk′ = σk2 .
k=K+1 F k=K+1

M ×N
For any B ∈ LK , we have σK+1 (B) = 0 because B has rank at most K.
Now apply Weyl’s inequality for singular values (6.71) with m = i − 1 and n = K to see that

σK+i (A) = σK+i (A − B + B) ≤ σi (A − B) + σK+1 (B) = σi (A − B) (7.4)

for i = 1, . . . , min(M, N ) − K. Thus B provides no better approximation than ÂK :

min(M,N ) min(M,N )−K
X X
||A − B|| 2F = σi2 (A − B) ≥ σi2 (A − B)
i=1 i=1
min(M,N )−K min(M,N )
X X
≥ 2
σK+i (A) = σk2 (A) = ||A − ÂK ||2F . 2 (7.5)
i=1 k=K+1

In words: the best rank (at most) K approximation to a matrix is given by its truncated SVD (to K terms).
© J. Fessler, September 2, 2023, 07:03 (class version) 7.9

Explore 7.1. Prove that ÂK = PUK A

(Read)
7.2.2 Subspace approximation perspective
 
Often we arrange a collection of N data points x1 , . . . , xN ∈ FM into a matrix X = x1 . . . xN ∈
FM ×N . From (7.2), the best rank-K approximation to this matrix is X ≈ X̂K = UK ΣK VK′ , based on the
SVD X = U ΣV ′ . Now look at this approximation from the persective of an individual data point, the nth
column of X:
XK
xn = X:,n = Xen ≈ x̂n ≜ X̂K en = UK ΣK VK′ en = ∗
(σk vnk ) uk ∈ R(UK ) . (7.6)
k=1

In words, data point xn has an approximation that lies in the K-dimensional subspace R(UK ). Geometrically,
we often think of the data points as being in a cloud and here we are approximating that cloud by a K-
dimensional hyperplane.

7.2.3 Implementation
For implementation, one must choose K ≤ r = rank(A) ≤ min(M, N ).
The approximation error is given by the 2-norm of the excluded singular values (σK+1 , . . . , σr ),
so it is useful to plot all the singular values versus k and look for a “knee in the curve.”
Low-rank approximation is related to factor analysis, and a scree plot shows the singular values [7].
© J. Fessler, September 2, 2023, 07:03 (class version) 7.10

Example. Fig. 7.1 shows scree plots for the MRI data in Fig. 7.3 and for a synthetic example for r = 3 with
noise. Hopefully your scree plot has a roughly linear section corresponding to the noise, shown as red squares
in Fig. 7.1.

Singular values for MRI coil compression example

1
<k

0.5
0.2
0
1 2 3 4 5 6 7 8
k

Figure 7.1: Scree plots from MRI data (left) and synthetic data (right).

For moderate sized matrices, low-rank approximation needs just a few lines of J ULIA code.
One approach uses the sum of outer products solution (7.2) literally:
/ 7.1
using LinearAlgebra: svd
(U, s, V) = svd(A)
B = +([U[:,k] * s[k] * V[:,k]' for k in 1:K]...)
© J. Fessler, September 2, 2023, 07:03 (class version) 7.11

Alternatively, one can use matrix multiplication corresponding to the compact SVD form in (7.2):
/ 7.2
using LinearAlgebra: svd, Diagonal
(U, s, V) = svd(A) # default full=false saves memory
B = U[:,1:K] * Diagonal(s[1:K]) * V[:,1:K]'

Explore 7.2. Why use small s instead of capital S above?

In both versions, the final B is ÂK .

Quality code would also include a bounds check that 1 ≤ K ≤ min(M, N ).

Yet another option is to compute only the first K SVD components using svds :
/ 7.3
using LinearAlgebra: Diagonal
using Arpack: svds
tmp = svds(A, nsv=K, ritzvec=true)[1] # special SVD data structure
(U, s, Vt) = (tmp.U, tmp.S, [Link]) # extract needed components
B = U * Diagonal(s) * Vt

Explore 7.3. Why use Vt above?

© J. Fessler, September 2, 2023, 07:03 (class version) 7.12

7.2.4 Choosing rank via permutation

Looking for a jump in the scree plot seems subjective. A more quantitative approach uses a permutation
method [8].
Consider a matrix that is the sum of a low-rank component and IID random noise. If we permute (randomly)
each of the columns (or rows? or all?) of that matrix, the noise part will remain the same (statistically) but
the low-rank structure will be destroyed. The singular values of this shuffled matrix will basically correspond
to noise. The singular values of the original data that are larger than those of the shuffled data might be
considered to be the “significant” ones and can help guide rank choice. See [8] for theoretical analysis. One
might need to permute multiple times and average to get a reliable baseline.
Another approach uses random signflips [9]. See also [10].
Here is J ULIA code for such permutations.
/ 7.4
using Random: shuffle
Ys = mapslices(shuffle, Y, dims=1) # shuffle each column of Y
Ys = reshape(shuffle(vec(Y)), size(Y)) # or, shuffle all

Example. Fig. 7.2 illustrates a case where the matrix Y = X + Z where X is a matrix with rank r = 2 and
Z is a Gaussian noise matrix.
© J. Fessler, September 2, 2023, 07:03 (class version) 7.13

Figure 7.2: This plot uses synthetic data that is rank 2 plus noise. The two largest singular values stand out
above the shuffled case.
© J. Fessler, September 2, 2023, 07:03 (class version) 7.14

MRI coil compression application

Example. Modern MRI scanners use multiple receive coils to collect more image information at the same
time. Processing many sets of 3D images can require undesirably large computation times, and typically
the data recorded by multiple coils is close to being linearly dependent. A technique called coil compression
reduces data size [11, 12], essentially by making a low-rank approximation to the data and then just storing
and processing the K singular vectors. Fig. 7.3 shows MRI brain images recorded simultaneously by 8 coils
around the head. The scree plot in Fig. 7.1 suggests that most of the information is in the first 4 components.
We reshape the data into a (M = 2562 × N = 8) matrix, use an SVD to find the rank-4 approximation, and
reshape the first K = 4 vectors u1 , . . . , u4 into 256 × 256 “virtual coil” images for display in Fig. 7.3. (In
practice this operation is done on the raw data before making images.)

Figure 7.3: Top: MRI brain images from 8 different receive coils.
Bottom: virtual coil images based on low-rank approximation with K = 4.
© J. Fessler, September 2, 2023, 07:03 (class version) 7.15

7.2.5 Non-uniqueness of SVD and low-rank approximation

As noted previously, “the” SVD of a matrix A is not unique. If two singular values are the same, then one 
can swap the corresponding columns of U and V (or perform more general rotations) and have two equally
Pmin(M,N )
valid SVD expressions. When we write A = k=1 σk uk vk′ , the order of the terms in the sum is unique
if σ1 > σ2 > . . . > σmin(M,N ) , but if any of the singular values are identical (including two or more zero
singular values) then the ordering is not unique, and one can swap corresponding columns of U and V .
Example. Here are two different SVDs for a diagonal matrix with 7 = σ1 > σ2 = σ3 = 5:
         
7 0 0 1 0 0 7 0 0 1 0 0 1 0 0 7 0 0 1 0 0
0 5 0 = 0 1 0 0 5 0 0 1 0 = 0 0 1 0 5 0 0 0 1 .
0 0 5 0 0 1 0 0 5 0 0 1 0 1 0 0 0 5 0 1 0
| {z } | {z } | {z } | {z } | {z } | {z } | {z }
A U Σ V′ Ũ Σ Ṽ ′
A subsequent example shows that there are two distinct rank K = 2 approximations to this A.

Even if all the singular values are distinct, “the” SVD is still not unique because we can multiply both uk and
vk by eıϕ and get a “different” U and V . However, the low-rank approximation is still the same in this case
because σk (eıϕ uk )(eıϕ vk )′ = σk uk vk′ .
We have not focused on this issue too much yet because often uniqueness is relatively unimportant.
© J. Fessler, September 2, 2023, 07:03 (class version) 7.16

Now consider the question of whether the rank-K approximation of A is unique. This is simpler to answer.
• If σK > σK+1 then one can extend the proof in (7.5) to show that the rank-K approximation is unique.
• If σK = σK+1 , then the rank-K approximation is not unique because we could swap the Kth and (K + 1)
terms.
Example. Here are two distinct rank-1 approximations to a simple diagonal matrix:
         
5 0 1 0 5 0 1 0 0 1 5 0 0 1
A≜ = =
0 5 0 1 0 5 0 1 1 0 0 5 1 0

         
1 1 0 5 0 0 0 1 0 0
A ≈ B1 ≜ σ1 u1 v1′ =5 = , A ≈ B2 ≜ σ1 ũ1 ṽ1′ =5 = .
0 0 0 1 0 5

Q7.2 What is the approximation error ||A − B1 ||F in the preceding example?
√
A: 0 B: 5 C: 5 D: 10 E: 52 ??

In practice, any rank-K approximation (for suitable K) often is equally useful, so one infrequently worries
about non-uniqueness.
Furthermore, in practical finite precision computing, it is dubious to test whether two floating-point values
σK and σK+1 are “exactly” equal. Avoid code like if (a == b) with floating-point numbers!
For analysis of the numerical stability of LR matrix approximation, including consideration of the “singular
© J. Fessler, September 2, 2023, 07:03 (class version) 7.17

value gap” σK − σK+1 , see [13].

7.2.6 1D example

 a 2×9 matrix where each column is an (xn , yn ) coordinate for n = 1, . . . , N = 9,

Example. Here we consider
x . . . x9
i.e., A = 1 .
y1 . . . y 9

Explore 7.4. Based on the blue dots in Fig. 7.4, what is rank(A)?

Suppose we want to perform dimensionality reduction of the 2D data by reducing it to a 1D line.

• One option for processing this data is to perform a rank-1 approximation, leading to dashed black line
in Fig. 7.4. The black points are the nearest points in the subspace for each of the data points.
• An alternative is to perform a linear least-squares fit assuming yn ≈ αxn . This leads to the green line
with slope α̂ in Fig. 7.4.

Demo 7.1 provides details.

© J. Fessler, September 2, 2023, 07:03 (class version) 7.18

Figure 7.4: LS vs rank-1 approximation.

© J. Fessler, September 2, 2023, 07:03 (class version) 7.19

It might seem surprising that these two methods lead to different fits.
The rank-1 approximation is
  2 N   2
x1 . . . xN X xn
arg min −B = arg min − B:,n .
B∈LM ×N y1 . . . y N F B∈LM ×N
yn 2
1 1 n=1

The linear least-squares estimator is

    2
y1 x1 N
 ..   ..  X
arg min  .  −  .  α = arg min |yn − α xn |2 .
α α
yN xN n=1
2

These two methods measure approximation error differently, as illustrated in Fig. 7.5.
© J. Fessler, September 2, 2023, 07:03 (class version) 7.20

Figure 7.5: Error measures for rank-1 (left) and LS (right) approximations.
© J. Fessler, September 2, 2023, 07:03 (class version) 7.21

7.2.7 Generalization to other norms

Fact 7.2. (Eckart-Young-Mirsky theorem; see [5] and [14] for a proof.)
For any unitarily invariant matrix norm ||·|| UI , the low-rank approximation
Pr problem has the same solution
using the first (largest) K singular components of A = U ΣV = k=1 σk uk vk′ :
′

XK
ÂK ≜ arg min ||B − A|| UI = σk uk vk′ . (7.7)
M ×N | k=1
B∈LK {z }
,→ unitarily invariant norm

Although the proof in (7.5) was for the Frobenius norm, remarkably the same result holds (with a different
proof) for other unitarily invariant norms such as the spectral norm.
To be robust to data outliers, §6.2.5 showed that using ∥·∥p for 1 ≤ p ≪ 2 is preferable to using ∥·∥2 for ♦♦
regression problems. Similarly, to perform robust low-rank approximation, we can replace the Frobenious
norm, which is equivalent to ∥vec(B − A)∥2 , with a cost function like ∥vec(B − A)∥p [15, 16]. Usually
such alternatives are not unitarily invariant, so (7.7) is inapplicable. They require different solution methods
(typically iterative algorithms). SVD operations are still used as one part of many such algorithms [16].

qWhat
Q7.3 is ||ÂK − A|| 2 when A has rank r ? Assume K < r ≤ min(M, N ).
Pr 2
B: rk=K+1 σk
P
A: k=K+1 σk C: σK+1 D: 0 E: None of these ??
© J. Fessler, September 2, 2023, 07:03 (class version) 7.22

Explore 7.5. What is ||ÂK − A|| ∗ ? (Choose from same answer list.)
(Read)
[Link] Proof for spectral norm
Claim: rank(B) ≤ K =⇒ ||A − ÂK ||22 ≤ ||A − B|| 22 where ÂK ≜ K ′
P
k=1 σk uk vk .
 
Proof . Let W = v1 . . . vK+1 so dim(R(W )) = K + 1 because columns of W are orthonormal. Now
rank(B) ≤ K =⇒ dim(N (B)) ≥ N − K =⇒ dim(N (B)) + dim(R(W )) ≥ (N − K) + (K + 1) = N + 1.
But both N (B) and R(W ) are subspaces in FN , so they must have a nonzero intersection, by (4.12).
Thus there exists some x ̸= 0 such that x ∈ N (B) and x ∈ R(W ) . WLOG take ∥x∥2 = 1.
x ∈ R(W ) =⇒ x = W W ′ x =⇒ 1 = ∥x∥2 = ∥W W ′ x∥2 = ∥W ′ x∥2 by (2.39). To complete the proof:
Xr 2
2
||A − ÂK ||22 = σk uk vk′ = σK+1
2 2
= σK+1 ∥W ′ x∥2 because ∥W ′ x∥2 = 1
k=K+1 2
K+1
X K+1
X X
2 2 2
2
= σK+1 ∥vk′ x∥2 ≤ σk2 ∥vk′ x∥2 ≤ σk2 ∥vk′ x∥2 = ∥Ax∥22
k=1 k=1 k
2
= ∥(A − B) x∥2 because x ∈ N (B)
≤ ||A − B|| 22 ∥x∥22 = ||A − B|| 22 . 2

Alternative proof based on Weyl’s inequality (6.71). Taking i = 1 in (7.4), for rank(B) ≤ K we have
||A − ÂK ||2 = σK+1 (A) ≤ σ1 (A − B) = ||A − B|| 2 . 2
© J. Fessler, September 2, 2023, 07:03 (class version) 7.23
(Read)
7.2.8 Bases for FM ×N
Recall that a set of linearly independent vectors b1 , b2 , . . . is a basis for a vector space V iff
span({b1 , b2 , . . .}) = V, P
i.e., we can write every vector in V as a linear combination of the basis vectors: v ∈ V =⇒ v = k bk αk
for some αk ∈ F. Now we consider two types of bases for the vector space of matrices FM ×N .

[Link] Canonical basis for FM ×N

First, the obvious (canonical) basis for FM ×N is:

B = {em ẽ′n : m = 1, . . . , M, n = 1, . . . N } ,

where em denotes the mth unit vector of length M and ẽn denotes the nth unit vector of length N . These are
clearly linearly independent and span(B) = FM ×N because
XM XN
A ∈ FM ×N =⇒ A = amn em ẽ′n .
m=1 n=1

To attempt “dimensionality reduction” using this basis, we could seek an approximation of the form
XM XN
B= bmn em ẽ′n , where we limit the number of nonzero coefficients bmn .
m=1 n=1
© J. Fessler, September 2, 2023, 07:03 (class version) 7.24

Let ∥vec(B)∥0 denote the number of nonzero elements of matrix B. Then a natural optimization problem is:

B̂ = arg min ||A − B|| F . (7.8)

M ×N
B∈F , ∥vec(B)∥0 ≤ K

This problem has a trivial solution: choose the K elements of A having the largest magnitudes |amn | and set
all other elements of B̂ to zero. The approximation error is the square root of the sum of squared magnitudes
of all those other elements.

Example. Consider the following “0-norm” formulation: (Read)

   
3 5 7 0 5 7
B̂ = arg min −B = .
M ×N 6 0 2 6 0 0
B∈F , ∥B∥0 ≤ 3 F

√
 
3 0 0
The approximation error is ||A − B̂|| F = = 22 + 32 .
0 0 2 F

Explore 7.6. How large must K be here for this approximation to have zero error?

Explore 7.7. How large must K be here for the rank-K approximation to have zero error?
© J. Fessler, September 2, 2023, 07:03 (class version) 7.25

The low-rank approximation seems more parsimonious than approximation using B.

To state this rigorously, we must first show that our low-rank approximation also uses a basis for FM ×N .

[Link] Orthonormal bases for FM ×N

Let U denote any unitary M × M matrix and V denote any unitary N × N matrix.
Now endow FM ×N with the (Frobenius) inner product: ⟨A, B⟩ = trace{B ′ A}
p p
and corresponding (Frobenius) norm: ||A|| F = ⟨A, A⟩ = trace{A′ A}.
Claim: the following set of M N rank-1 M × N outer-product matrices is an orthonormal basis for FM ×N :
{um vn′ : m = 1, . . . , M, n = 1, . . . N } . (7.9)
Proof of orthonormality (for the Frobenius inner product):
′
⟨um vn′ , uk vl′ ⟩ = trace (uk vl′ ) (um vn′ ) = trace{vl u′k (um vn′ )}


′ ′ ′ ′ 1, k = m, l = n
= trace{uk um vn vl } = (uk um )(vn vl ) =
0, otherwise.
Being orthonormal, the set {um vn′ } spans a M N -dimensional
PM space
Pand is a basis for FM ×N . One can write
M ×N N
any X ∈ F in terms of that basis: i.e., X = U CV = m=1 n=1 cmn um vn′ , where C ≜ U ′ XV .
′

When working with a matrix A, using the basis for FM ×N that is formed from its own left and right singular
vectors, i.e., U and V , respectively, turns out to provide the most parsimonious representation. Note that the
original low-rank problem formulation (7.1) did not involve any SVD, but an SVD arose in the solution.
© J. Fessler, September 2, 2023, 07:03 (class version) 7.26

Q7.4 Does the canonical basis {em ẽ′n : m = 1, . . . , M, n = 1, . . . N } form an orthonormal basis for
FM ×N when using the Frobenius inner product?
A: Yes B: No C: Insufficient information ??

7.2.9 Low-rank approximation summary

The low-rank approximation problem (7.1) and its SVD-based solution (7.2) provides the best low-rank
representation (or approximation) of a given matrix (in
Pthe Specifically, if B is any
Frobenius norm sense). P
rank-K matrix having the same size as a matrix A = k=1 σk uk vk then for ÂK = K
r ′ ′
k=1 σk uk vk we have
v
u r ≤ ||A − B|| F , ∀B ∈ LM ×N
u X
K
t σk2 = ||A − ÂK ||F ×N (7.10)
< ||A − B|| F , ∀B ∈ LM K \ {ÂK }, σK > σK+1 .
k=K+1

Explore 7.8. The definition in (7.1) describes low-rank matrix approximation as a projection onto the (non-
convex) set of rank-(at most)-K matrices. In contrast, [17] defines the following operation A 7→ UK′ AVK .
Discuss whether that operation is a projection operator.
© J. Fessler, September 2, 2023, 07:03 (class version) 7.27
(Read)
[Link] Generalizations
What if the data is corrupted by noise? What if
Y =X +Z
where X is a low-rank matrix (or close to low rank). Due to the noise matrix Z, typically Y is not low-rank.
We really want to recover X here, not just represent or approximate Y . How do we do that? An answer
called OptShrink is given in [18] and discussed in §7.6.2.
Low-rank approximation is related closely to principal component analysis (PCA), discussed in §7.7.2.
There are many PCA generalizations, including robust methods (see §10.5.1), nonlinear models, multi-linear
(tensor) models, non-negative matrix factorization (NMF) (see §10.6) [19, 20], methods that use sparsity,
etc.
(Read)
7.2.10 Rank and stability
Although rank is a simple mathematical concept, it is not stable numerically. Precisely, it is not a continuous
function of the elements of a matrix.
Example. Consider the following two matrices for a ̸= 0:
   
a a a a
A= , B= , 0 < |ε| ≪ 1.
a a a(1 + ε) a
© J. Fessler, September 2, 2023, 07:03 (class version) 7.28

These two matrices are nearly indistinguishable numerically, yet rank(A) = 1 and rank(B) = 2.

An alternative matrix property used in some situations is the stable rank [21], defined (for A ̸= 0) as:

||A|| 2F ||A|| 2F 1 Xmin(M,N ) 2 Xmin(M,N ) σ 2

k
2 , where 1 ≤ 2 = 2
σ k = 1 + ≤ r, (7.11)
||A|| 2 ||A|| 2 σ1 k=1 k=2 σ12
√ 4a2 +2ϵa2
because ( HW ) ||A|| 2 ≤ ||A|| F ≤ r ||A|| 2 . For the above example, ||B|| 2F /|| B|| 22 ≈ 4a2
= 1 + 2ε .

||A|| 2F 22 + 12
 
1 0
Example. For A = , the stable rank is = = 1.25 ∈ [1, 2].
0 2 ||A|| 2
2
22

Q7.5 What is the stable rank of A = xy ′ , when x and y are nonzero vectors?
A: 0 B: 1 C: ∥x∥2 ∥y∥2 D: ∥x∥22 ∥y∥22 ??
© J. Fessler, September 2, 2023, 07:03 (class version) 7.29

7.2.11 Example: Photometric stereo

In computer vision, photometric stereo is a set of methods for estimating the surface normals of 3D objects
from a set of (often grayscale) images of that object taken with different lighting conditions. The original
method assumed knowledge of the lighting directions [22], whereas later methods relaxed this assumption
[23, 24]. Extensions include handling shadows and specular reflections [17, 25–27]. For a Lambertian
surface illuminated by a distance point light source, the reflected light intensity at a point on the object
is proportional to the inner product between the surface normal vector and the normalized light direction
vector. Given N images (for N light sources), each with M pixels, we create a M ×N data array Y where the
nth column of Y is the vec of the nth image. The Lambertian model is Y = XL where X is the (unknown)
matrix of M × 3 surface normals, and L is the 3 × N matrix of light directions. Ideally (in the absence of
noise and shadows, etc.), the matrix Y has (at most) rank 3, so one can use a rank-3 factorization to estimate
the surface normals and the lighting directions.

Demo 7.2 illustrates the input and output data for photometric stereo; see Fig. 7.6.
That demo also uses the Procrustes method from §6.7 to align coordinate systems.
© J. Fessler, September 2, 2023, 07:03 (class version) 7.30

Figure 7.6: Illustration of photometric stereo. Left: 12 input images of a static 3D object taken with different
lighting directions. Data from [26]. Right: x, y, and z components of true and estimated surface normals.
© J. Fessler, September 2, 2023, 07:03 (class version) 7.31

7.3 Sensor localization application: multidimensional scaling

We now turn to another application of matrix factorization: sensor localization via multidimensional scal-
ing (MDS). Suppose J ≥ d sensors are located at unknown locations c1 , . . . , cJ ∈ Rd , where typically d = 2
or d = 3. All we are given is the J × J distance matrix D having elements

dij = dji = ∥ci − cj ∥2 , i, j = 1, . . . , J. (7.12) c2

| {z }
,→ usual Euclidean distance c4
c3
Clearly by definition the diagonal elements of D are zero: djj = 0. d14
The goal of MDS is: given D, determine the locations {cj }. c1

Figure 7.7: Example sen-

sor set.
There is a fundamental ambiguity because translating or rotating or reflecting coordinates does not change
D. Specifically, if c̃j = Qcj + d where Q is a d × d orthogonal matrix (such as a rotation matrix) and d ∈ Rd
is a displacement vector, then

∥c̃i − c̃j ∥2 = ∥(Qci + d) − (Qcj + d)∥2 = ∥Q(ci − cj )∥2 = ∥ci − cj ∥2 , (7.13)

because the Euclidean norm is unitarily invariant. So D would be the same for {c̃j } and {cj }.
Thus we must be content with determining locations {cj } to within some translation and rotation factor.
© J. Fessler, September 2, 2023, 07:03 (class version) 7.32

Note the choice to use a matrix representation of the data!

And of course the locations are vectors.

7.3.1 Derivation (analysis)

To derive a solution (for noiseless data), define another matrix S by the square of the elements of D:
sij ≜ d2ij = ∥ci − cj ∥22 = ∥ci ∥22 + ∥cj ∥22 − 2 ⟨ci , cj ⟩ . (7.14)
Naturally we write S in matrix form:
∥c1 ∥22 . . . ∥c1 ∥22 ∥c1 ∥22 . . . ∥cJ ∥22
       
s11 . . . s1J ⟨c1 , c1 ⟩ . . . ⟨c1 , cJ ⟩
S =  ... . . . ...  =  ... ..  +  .. ..  − 2  .. ..
   
... .   . ... .   . ... . 
2 2 2 2
sJ1 . . . sJJ ∥cJ ∥2 . . . ∥cJ ∥2 ∥c1 ∥2 . . . ∥cJ ∥2 ⟨cJ , c1 ⟩ . . . ⟨cJ , cJ ⟩
∥c1 ∥22
   
c′1
=  ...  1′J + 1J ∥c1 ∥22 . . . ∥cJ ∥22 − 2  ...  c1 . . . cJ
      

∥cJ ∥22 c′J

∥c1 ∥22
 

= r1′J + 1J r ′ − 2C ′ C, r ≜  ...  ,
 
C ≜ c1 . . . cJ . (7.15)
 
2
∥cJ ∥2
| {z }
d × J unknown locations
Unfortunately we do not know the values in the vector r.
© J. Fessler, September 2, 2023, 07:03 (class version) 7.33

So we know the J × J matrix S and we want to solve for the d × J unknown location matrix C using

S = r1′J + 1J r ′ − 2C ′ C. (7.16)

Because S on the left is symmetric and zero along its diagonal, it contains (J 2 − J)/2 distinct known values.
The right-hand side depends on the dJ unknown values in C.
Because source localization has a translation ambiguity, we may as well use a coordinate system where the
centroid is 0, i.e., we assume C1J = Jj=1 cj = 0 without losing any further generality.
P

Following (6.103) define the following “de-meaning” operator that projects onto the orthogonal complement
of R(1J ):
1
P ⊥ ≜ I − 1J 1′J . (7.17)
J
Clearly P ⊥ 1J = 0, so 1′J P ⊥ = 0′ and CP ⊥ = C(I − J1 1J 1′J ) = C so from (7.16):

P ⊥ SP ⊥ = P ⊥ (r1′J + 1J r ′ − 2C ′ C) P ⊥ = −2C ′ C.

Rearranging leads to the following simple expression for the J × J Gram matrix:

1
G ≜ C ′ C = − P ⊥ SP ⊥ . (7.18)
2
In words, we remove the row and column means of S and divide by −2.
© J. Fessler, September 2, 2023, 07:03 (class version) 7.34

Now we have the Gram matrix G = C ′ C, but we really want the coordinates matrix C itself.
Although C ′ C is a J ×J matrix, d ≤ J so typically rank(C ′ C) = d (in the absence of noise). If the locations
happen to be linearly dependent, i.e., if the sensors are along one line through the origin in 2D or in the same
plane through the origin in 3D, then rank(C ′ C) < d.
To elaborate, consider d = 2 where
 
  x1 . . . xJ
C ≜ c1 . . . cJ = ,
y1 . . . y J

where cj = (xj , yj ). Because C has two rows, rank(C ′ C) = rank(C) ≤ 2.

Explore 7.9. Determine when rank(C) = 1.

The Gram matrix C ′ C consists of inner products, and it is invariant to rotations of the source locations:

C̃ = QC =⇒ C̃ ′ C̃ = (QC)′ (QC) = C ′ Q′ QC = C ′ C.
© J. Fessler, September 2, 2023, 07:03 (class version) 7.35

Given the (noiseless) J × J Gram matrix G ≜ C ′ C = − 21 P ⊥ SP ⊥ , we determine its compact SVD:

J
X d
X
G = V ΣV ′ = σk vk vk′ =
|{z} σk vk vk′ = Vd Σd Vd′ , (7.19)
k=1 (noiseless) k=1

 
σ1
 
| |
Σd ≜ 
 .. ,

Vd ≜ v1 . . . vd  .
.
σd | |
| {z } | {z }
d×d J ×d

• Because G is (symmetric) positive semidefinite (in the absence of noise), we know U = V .

• Because G has rank at most d (in the absence of noise), the above low-rank form is exact, not an approx-
imation like we used earlier in this chapter.
We define our estimate of the source locations to be the J columns of this d × J matrix:
√  
σ1 − v1′ −
1/2
Ĉ = Σd Vd′ =  ... ..
, (7.20)
  
.
√

′
σd − vd −
| {z }| {z }
d×d d×J
© J. Fessler, September 2, 2023, 07:03 (class version) 7.36

because this estimate Ĉ is consistent with C to within an unknown rotation and translation:
1/2 1/2
Ĉ ′ Ĉ = (Σd Vd′ )′ (Σd Vd′ ) = Vd Σd Vd′ =
|{z} V
| ΣV
′ ′
{z } = G = C C. (7.21)
(noiseless) !

Now it seems like we are done because we have an expression for G that is computable from S, so we can
use its SVD to find the source location estimates Ĉ. However, we assumed that C1J = 0 so we must verify
that Ĉ1J = 0, to ensure that our approach is self consistent.
Because P ⊥ 1J = 0, it follows from (7.17) that G1J = 0, so V ΣV ′ 1J = 0. In the noiseless case, from
1/2
(7.21) we then have Vd Σd Vd′ 1J = 0, which in turn implies that Ĉ1J = Σd Vd′ 1J = 0.

Q7.6 What is the size of the final 0 in the immediately preceding line?
A: J × J B: J C: J × d D: d × J E: d ??
© J. Fessler, September 2, 2023, 07:03 (class version) 7.37

7.3.2 MDS method

Summary of MDS method for finding locations from distances (in the possible presence of noise):
• Given distances arranged in a matrix D.
• Compute S from D by squaring elements.
• Compute Ĝ ≜ − 12 P ⊥ SP ⊥ ≈ C ′ C by de-meaning.
• Optional: force Ĝ to be symmetric by replacing it with (Ĝ + Ĝ′ )/2
• Compute the compact SVD Ĝ ≈ Ud Σd Vd′ . This is a rank-d matrix approximation in the presence of
noise (the topic of this chapter). (If Ĝ is symmetric then Ud = Vd .)
1/2
• Use the first d terms for source location estimate: Ĉ = Σd Vd′ .
This remarkably simple algorithm for finding Ĉ is called (classical) multidimensional scaling (MDS) [28,
Ch. 12], and using the low-rank compact SVD is a key step. Because G = C ′ C, matrix G is positive
semidefinite (in the absence of noise). If we force Ĝ to be likewise symmetric, then its eigenvalues are all
real, so alternatively one could use an eigendecomposition of Ĝ with suitably ordered eigenvalues.
Example. To see how symmetrizing a matrix does not ensure nonnegative eigenvalues in general (in the
presence of noise), consider the (asymmetric, but still de-meaned) matrix:
   
0 2 −2 0 3 −3
Ĝ =  4 0 −4 =⇒ (Ĝ + Ĝ′ )/2 =  3 0 −3 .
−4 −2 6 −3 −3 6
© J. Fessler, September 2, 2023, 07:03 (class version) 7.38

One can verify that the eigenvalues of (Ĝ + Ĝ′ )/2 are {−3, 0, 9}.

This situation is why the SVD is preferable to using an eigendecomposition for MDS.
Questions

Q7.7 Consider a given noiseless distance matrix D with J > d. The position estimates Ĉ returned
by the SVD-based MDS algorithm are unique (to within numerical precision), i.e., are the same for
any SVD of the Gram matrix.
A: True B: False ??

Q7.8 In the possible presence of noise, the product Ĉ ′ Ĉ is unique (to within numerical precision)
for any SVD of the calculated Gram matrix Ĝ.
A: never B: usually C: always ??
© J. Fessler, September 2, 2023, 07:03 (class version) 7.39

7.3.3 Practical implementation

J ULIA’s [Link] has a classical_mds command.
But the method is so simple that it is just as easy to code our own.
 
0 1 1
Example. Consider J = 3 sources that are all equally distant: S = D .^2 = 1 0 1 .
1 1 0
Just a few lines of J ULIA suffices (see demo code linked below):
/ 7.5
using Plots; using LinearAlgebra: svd, Diagonal
using Statistics: mean
D = [0 1 1; 1 0 1; 1 1 0] + 0e-9 * randn(3,3) # given distances
S = D.^2
tmp = S .- mean(S, dims=1)
G = -0.5 * (tmp .- mean(tmp, dims=2))
(_, s, V) = svd(G)
C = Diagonal(sqrt.(s[1:2])) * V[:,1:2]'
scatter(C[1,:], C[2,:], aspect_ratio=1, label="",
xaxis=("x", (-0.6,0.6), -.5:.5:.5),
yaxis=("y", (-0.6,0.6), -.5:.5:.5), title="s=$(s)")
© J. Fessler, September 2, 2023, 07:03 (class version) 7.40

Explore 7.10. What “shape” should appear?

The SVD has non-trivial non-uniqueness in this case:

   p p   
2 −1 −1 −p 2/3 p0 1/3 1/2 0 0
1
G = −1 2 −1 = V ΣV ′ , V =  p1/6
p
p1/2 p1/3 , Σ = 0 1/2 0
  
6
−1 −1 2 1/6 − 1/2 1/3 0 0 0
p p p 
1/2 1/6 p1/3
′
p
= Ṽ ΣṼ , Ṽ = 
p0 −2
p 1/6 p1/3 .

− 1/2 1/6 1/3

These decompositions lead to two (of many possible) different location estimates:

c2
c1
 p p p 
1/2
p − 2/3 p1/6 1/6 c3
Ĉ = Σ2 V2′= 1/2 p
0 1/2 − 1/2
p p  c3 c1
ˆ 1/2 ′
p 1/2 0 −p 1/2 . c2
C̃ = Σ2 Ṽ2 = 1/2 p p
1/6 −2 1/6 1/6
Figure 7.8: MDS results.
© J. Fessler, September 2, 2023, 07:03 (class version) 7.41

Demo 7.3 considers an unknown collection of J = 80 source points. See Fig. 7.9 for images of D and G.
Run the demo code to see Ĉ and consider noisy cases to examine robustness.

Figure 7.9: Images of D and G for source localization demo.

© J. Fessler, September 2, 2023, 07:03 (class version) 7.42

7.3.4 Extensions
• When the distance measurements are noisy then D may not be perfectly symmetric and P ⊥ SP ⊥ may
have (hopefully small) negative eigenvalues.
Use truncated SVD: truncate rank to dimension d of embedding, or shrink the eigenvalues [29].
• What if some sensors cannot reach some other sensors, i.e., if D has some missing elements?
• Is there a way to formulate the problem that considers the possibility of noise (errors in the distance values)
at the outset, instead of simply empirically investigating the effects of noise on our solution.
• What if there are outliers in the data, e.g., some sensors that give incorrect values (either due to errors or
because the sensor was hacked or otherwise compromised). Is MDS robust to such errors, or is a different
formulation needed, perhaps based on a robust distance measure like ∥·∥1 ?
• If we want to compare the location estimates Ĉ to some ground truth, we can use the Procrustes method
to align them, compensating for the inherent rotation/translation ambiguity.
Multidimensional scaling is used to track flu and corona virus mutations [30]. More generally, MDS is used
as a dimensionality reduction method, cf. §8.9.6.
Explore 7.11. What happens if D is geodesic distances or distances along roads between locations?

Explore 7.12. The derivation here was for the usual case where C ∈ Rd×J .
Suppose C ∈ Cd×J . Does the derivation extend to this case? If not, which step does not generalize?
© J. Fessler, September 2, 2023, 07:03 (class version) 7.43

7.4 Proximal operators

Our next topic is some alternative low-rank approximation methods, also formulated as optimization prob-
lems. To describe those alternatives, we first need to introduce an important operation.
Define. The proximal operator [31, 32] (also called the proximal mapping [33, Ch. 6]) associated with a
(typically convex) function f : FN 7→ R, is defined, for any v ∈ FN , as the following minimizer:
1
proxf (v) ≜ arg min ∥v − x∥22 + f (x). (7.22)
x∈FN 2
• The norm squared is a strictly convex function, so when f is convex the “arg min” is unique ∀v ∈ FN .
• Often the function f is additively separable, in which case the “arg min” separates into N individual 1D
minimization problems. See example below.
• The proximal operator transforms one function f (·) into another function proxf (·).

Q7.9 In general, the functions f and proxf have the same (domain?, co-domain?)
A: true,true B: true,false C: false,true D: false,false ??

??
© J. Fessler, September 2, 2023, 07:03 (class version) 7.44

7.4.1 Soft thresholding

Example. The case F = R and f (x) = β ∥x∥pp = β

PN
n=1 |xn |p , for β ≥ 0, is especially important.
The most important case is when p = 1, for which
N  
1 2
X 1 2
proxf (v) = arg min ∥v − x∥2 + β ∥x∥1 = arg min (vn − xn ) + β |xn | .
x∈RN 2 x∈RN n=1 2

This sum is additively separable, meaning each xn appears in just one term of the sum. So we can solve the
minimization problem separately for each xn term by solving
1
x̂n ≜ arg min g(xn ), gn (xn ) ≜ (vn − xn )2 + β |xn | , n = 1, . . . , N. (7.23)
xn ∈R 2

Often for minimization problems we take the derivative and set it to zero. That approach does not quite work
here because |·| is not differentiable at 0. So we need a “braces and cases” approach.
Suppose the minimizer satisfies x̂n > 0. In this regime the function g is differentiable with derivative

ġn (xn ) = xn − vn + β1.

Equating to zero yields x̂n = vn − β. However, we must keep in mind that we assumed here that x̂n > 0, so
this solution is valid only when vn > β.
© J. Fessler, September 2, 2023, 07:03 (class version) 7.45

Now suppose the minimizer satisfies x̂n < 0. In this regime the function g is differentiable with derivative

ġn (xn ) = xn − vn + β(−1).

Equating to zero yields x̂n = vn + β. Here we assumed that x̂n < 0, so this solution is valid only when
vn < −β. Because we have considered all cases where the minimizer is nonzero, for any other input value
vn , the minimizer must be zero. Thus we have shown that

soft(v; β)
1
arg min (vn − xn )2 + β |xn | = soft(vn ; β)
xn ∈R 2

 vn − β, vn > β −β β v
soft(vn ; β) ≜ vn + β, vn < −β (7.24)
0, otherwise.

Figure 7.10: Soft thresholding.
This operation is called soft thresholding and it induces sparsity because of the range of values set to zero.
Putting it all together, we have shown that the proximal operator for the 1-norm scaled by β is element-wise
soft thresholding:
proxβ∥·∥1 (v) = soft .(v; β). (7.25)
(The . here in soft. denotes Julia-style element-wise broadcast.)
© J. Fessler, September 2, 2023, 07:03 (class version) 7.46

Fig. 7.11 illustrates the cases where v < β and v > β, by plotting the following functions:
1
(x − v)2 + β |x|.
2
• When 0 < β ≤ v the minimizer is at x = v − β.
• When −β ≤ v ≤ β, the minimizer is at x = 0.

Figure 7.11: Soft thresholding component functions.

© J. Fessler, September 2, 2023, 07:03 (class version) 7.47

Fact 7.3. The generalization of the proximal operator (7.24) to complex numbers is ( HW ):
1
arg min |s − z|2 + β |z| = [|z| − β]+ eı∠z (7.26)
s∈C 2
[x]+ ≜ max(x, 0). (7.27)

7.4.2 Hard thresholding

Example. Another main case of interest here is the (nonconvex) 0-norm f (x) = β ∥x∥0 , for β ≥ 0, for
which the proximal operator is (element-wise) hard thresholding. Here is the derivation:
N  
1 2
X 1 2
proxβ∥·∥0 (v) = arg min ∥v − x∥2 + β ∥x∥0 = arg min |vn − xn | + βI{xn ̸=0} = hhard .(v; β),
x∈FN 2 x∈FN n=1 2

hhard (v; β)
1
hhard (v; β) ≜ arg min |v − x|2 + βI{x̸=0}
x∈F 2
v, 12 |v|2 > β
  p 
= = vH |v| − 2β , (7.28) √ v
0, otherwise 2β

Figure 7.12: Hard thresholding.

© J. Fessler, September 2, 2023, 07:03 (class version) 7.48

where H(t) = I{t>0} denotes the Heaviside step function. We cannot derive the minimizer by differentiation,
because the indicator function is not differentiable at zero. So we must apply “braces and cases.”
If the minimizer is nonzero, then the second term is β and the minimizer of the first term is v, and the overall
cost is 0 + β = β. If the minimizer is 0, then the first term is (v − 0)2 /2 and the second term is 0, so the cost
is v 2 /2. So 0 is the minimizer when v√
2
/2 < β, otherwise v is the minimizer.
The minimizer is√not unique if |v| = 2β; in such cases, both v and 0 are global minimizers.
The threshold is 2β due to the 1/2 in front of the norm in the proximal operator definition.

Now we return to low-rank approximation where we will use these results.

© J. Fessler, September 2, 2023, 07:03 (class version) 7.49

7.5 Alternative low-rank approximation formulations

If Y ∈ FM ×N is a given data matrix with SVD Y = U ΣV ′ and Y = X+Z where we believe rank(X) ≤ K
and Z denotes a M × N noise matrix, then the constrained formulation discussed so far is natural [4]:
K
X r
X
X̂ = arg min ||Y − X|| 2F = UK ΣK VK′ = σk uk vk′ = σk I{k≤K} uk vk′ . (7.29)
×N
X∈LM
K k=1 k=1

This problem is nonconvex, but still has a nice solution (7.2), reiterated in (7.29) in three different forms. The
final form reminds us that we have “set to 0” all the singular values for k > K.

7.5.1 Unconstrained / regularized formulation

When the rank is unknown, an alternative approach is to consider an unconstrained cost function that dis-
courages (penalizes) high rank (high complexity) rather than constraining it [34]:
r
1 X
X̂ = arg min ||Y − X|| 2F + β rank(X) = hhard (σk ; β) uk vk′ , (7.30)
X∈F M ×N 2 k=1

where (see proof sketch below) hhard (·; β) is the hard-thresholding function defined in (7.28).
© J. Fessler, September 2, 2023, 07:03 (class version) 7.50

7.5.2 General unitarily invariant formulations

There are several ways to form low-rank approximations that balance between data misfit and model com-
plexity, all having the general form:
1
X̂ = arg min ||Y − X|| 2UI + βR(X), (7.31)
X∈FM ×N 2
for some unitarily invariant matrix norm ||·|| and for some regularizer R(X), where R : FM ×N 7→ R, such
as R(X) = rank(X). The regularization parameter (aka hyperparameter) β > 0 controls the trade-off
between fit to the data Y , as measured by ||Y − X|| , and model complexity, as quantified by R(X).
We assume hereafter that the regularizer is also a unitarily invariant function, i.e., R(U XV ) = R(X) for
any suitably sized unitary matrices U and V . When both the data misfit norm and the regularizer are unitarily
invariant norms (or monotonically increasing functions thereof [14, Remark 3]), using the symmetric gauge
principles ( HW ) of [5, 35, 36] [14, Thm. 4], one can use the compact SVD Y = Ur Σr Vr′ to show that any ♦♦
minimizer of (7.31) has the form [14, Cor. 2]:
r
X
X̂ = ŵk uk vk′ = Ur Σ̂r Vr′ where Σ̂r = Diag{ŵk }, ŵk = hk (σ1 , . . . , σr ; β), (7.32)
k=1

for some shrinkage or thresholding function hk (·; β) that depends on the data misfit norm and the regular-
izer.
© J. Fessler, September 2, 2023, 07:03 (class version) 7.51

Q7.10 What is the rank of the X̂ in (7.32) when rank(Y ) = r ? (Choose best answer.)
A: Usually 1 B: Always 1 C: Usually r D: Always r E: None of these

??
Because (7.32) is the form of the solution, we can rewrite the general optimization problem (7.31) as
1
X̂ = Ur Σ̂r Vr′ , Σ̂r = arg min ||Σr − S|| 2UI + βR(S) = Diag{ŵ1 , . . . , ŵr } . (7.33)
S=Diag{s1 ,...,sr } 2

7.5.3 Singular value hard thresholding

Now return to the special case (7.30) with the Frobenius norm and the rank penalty. The equivalent mini-
mization problem in terms of the singular values requires minimizing over S = Diag{s1 , . . . , sr }:
r  
1 2
X 1 2
||Σr − S|| F + β rank(S) = (σk − sk ) + βI{sk ̸=0}
2 k=1
2

1 2
 p 
=⇒ ŵk = arg min (σk − s) + βI{s̸=0} = hhard (σk ; β) = σk H σk − 2β . (7.34)
s 2
We started with the cost function (7.30), simplified it to the equivalent problem with diagonal matrices (7.33),
and now we have r separate 1D problems that we solved earlier when deriving the proximal operator for the
0-norm that resulted in hard thresholding; see (7.28).
© J. Fessler, September 2, 2023, 07:03 (class version) 7.52

This analysis leads to ŵk = hhard (σk ; β) in the method (7.30) that is called singular value hard threshold-
ing.

7.5.4 Singular value soft thresholding

Any formulation involving rank(·) is non-convex, and it is somewhat remarkable that a nice solution to (7.30)
exists despite that non-convexity. Often it can be preferable to have a convex formulation. The convex
relaxation [37, 38] of rank is the nuclear norm that leads to soft thresholding of singular values:

1
X̂ = arg min ||Y − X|| 2F + β ||X|| ∗
X∈FM ×N 2
r 
X t, t > 0,
= [σk − β]+ uk vk′ , where [t]+ = (7.35)
0, otherwise.
k=1

So ŵk = hsoft (σk ; β), hsoft (σ; β) ≜ [σ − β]+ = max(σ − β, 0).

This approach is called singular value soft thresholding (SVST). See Fig. 7.13.
Proof sketch: 21 ||Σr − S|| 2F + β ||S|| ∗ = rk=1 12 (σk − sk )2 + βsk =⇒ ŵk = arg mins 21 (σk − s)2 + β |s| .
P
This is just the proximal operator of β |·| that we solved following (7.23), leading to ŵk = [σk − β]+ .
© J. Fessler, September 2, 2023, 07:03 (class version) 7.53

hsoft (σ; β)

β σ Figure 7.13: SVST function.

Q7.11 If Y has nonzero singular values 3, 5, 6, 7, 9 and β = 6, what is the rank of X̂ here?
A: 1 B: 2 C: 3 D: 4 E: 5 ??

Q7.12 This X̂ equals the best rank-K approximation to Y , where K is answer to previous problem.
A: True B: False ??

Q7.13 In general, what is ||X̂|| 2 in terms of the notation used throughout here?
√
A: β B: σ1 C: σ1 − β D: σr − β E: [σ1 − β]+ ??

Q7.14 In general, if ||Y ||2 ≥ β, then what is ||Y − X̂|| 2 in terms of the notation used throughout
here? √
A: β B: σ1 C: σ1 − β D: σr − β E: [σ1 − β]+ ??
© J. Fessler, September 2, 2023, 07:03 (class version) 7.54

One motivation for these unconstrained formulations is for solving matrix completion or matrix sensing
problems where we must replace ||Y − X|| 2F with something like ∥y − A vec(X)∥22 where A is a known
matrix. See §[Link]. These problems do not have closed-form solutions and require iterative methods, and
convergence of iterative methods is better understood for convex problems.

7.5.5 Other extensions of low-rank approximation

• Sometimes we might want to retain some columns of the data exactly, leading to a different type of
constrained low-rank approximation [39].
• Another variation replaces the Frobenius norm by the spectral norm: ( HW )
1
arg min ||Y − X|| 22 + β rank(X)
X 2
1
arg min ||Y − X|| 22 + β ||X|| ∗ .
X 2
• Another variation replaces the Frobenius norm by the nuclear norm:

arg min ||Y − X|| ∗ + β ||X|| ∗ .

X

The solution to this variation is left as an exercise. It seems not to be useful.

• For further generalizations using weighted norms, see [40–42].
© J. Fessler, September 2, 2023, 07:03 (class version) 7.55

Q7.15 If we use two Frobenius norms, the solution has the following general form:
r
1 1 X
X̂ = arg min ||Y − X|| 2F + β ||X|| 2F = ŵk uk vk′ , ŵk = hFF (σk , β). (7.36)
X 2 2 k=1

What is hFF (σ, β)?

√ σ
A: [σ − β]+ B: σH(σ − 2β) C: max(σ, β) D: E: (σ − β)2
1+β

??

Q7.16 If Y is M × N with rank r, then what is the rank of the solution X̂ here?
A: 0 B: 1 C: r D: min(N, M ) E: None of these.

??
© J. Fessler, September 2, 2023, 07:03 (class version) 7.56

7.6 Choosing the rank or regularization parameter

In all the above low-rank approximation formulations, we must either
◦ choose the rank K directly, or
◦ choose the regularization parameter β that influences the rank of X̂.
As illustrated in Fig. 7.14, this selection process is non-trivial because, when Y = X + Z and the latent
matrix X has low-rank:
• ||X̂ − Y ||F : the difference between the noisy data Y and the low-rank approximation X̂ always
◦ decreases monotonically as K increases (discretely)
◦ increases monotonically as β increases (continuously).
• ||X̂ − X|| F : the difference between the latent low-rank matrix X and the low-rank approximation X̂:
◦ usually decreases initially as K increases, then increases
◦ usually decreases initially as β increases, then increases.
Often we normalize the difference between X̂ and Y or that between X̂ and X and define the (unitless) nor-
malized root mean-squared difference (NRMSD) and normalized root mean-squared error (NRMSE):

||X̂ − Y ||F ||X̂ − X|| F

NRMSD = · 100%, NRMSE = · 100%. (7.37)
||Y ||F ||X|| F

(There are many variations of these definitions in the literature.)

© J. Fessler, September 2, 2023, 07:03 (class version) 7.57

Demo 7.4 considers a case where the 100 × 50 matrix X has rank 5 and Y = X + Z, where Z is M × N
matrix with IID Gaussian noise. See Fig. 7.14.

Figure 7.14: Trade-offs between NRMSE and NRMSD as one varies rank K or regularization parameter β.

Q7.17 What is the y-axis value of the second blue curve in Fig. 7.14 at β = 0?
A: ||Y − X|| F / ||Y ||F · 100% B: ||Y − X|| F / ||X|| F · 100% C: Neither ??

Ideally we would like to find K∗ or β∗ that minimizes ||X̂ − X|| F , the minimizers of the blue curves in
© J. Fessler, September 2, 2023, 07:03 (class version) 7.58

Fig. 7.14. Finding either of the choices exactly is impossible because X is unknown in practice!
All we have is Y = X + Z, and usually we are hoping that X is low rank, without knowing for sure.

7.6.1 Stein’s unbiased risk estimate (SURE)

Researchers have developed surrogates for ||X̂ − X|| F that help choose β.
One notable method is Stein’s unbiased risk estimate (SURE).
If x̂β = x̂(y; β) is a weakly differentiable estimate of a parameter x ∈ Rd from data with Gaussian noise:
y ∼ N(x, σ02 I), then an unbiased estimate of the risk, i.e., the mean-squared error (MSE), is:

divergence{x̂(·; β)} ≜ trace{∇y x̂(y; β)}

↰
z }| {
X ∂
SURE{β} ≜ ∥x̂(y; β) − y∥22 − d σ02 + 2 σ02 x̂i (y; β) (7.38)
i
∂y i
| {z }
Independent of x!
E[SURE{β}] = MSE{β} = E ∥x̂(y; β) − x∥2 .
 
(7.39)

For a proof, see [wiki] or [43]. Note that σ02 is a noise variance here, not a singular value!

Explore 7.13. Verify (7.39) for the case where x̂(y; β) = Uβ Uβ′ y, where Uβ is a unitary basis for a β-
dimensional subspace of Rd for β ∈ N.
© J. Fessler, September 2, 2023, 07:03 (class version) 7.59

The SURE approach uses the following approximation to select β:

β∗ ≜ arg min ∥x̂β − x∥2 ≈ arg min E ∥x̂β − x∥2 = arg min SURE{β} .
 
(7.40)
β β β

The hard thresholding function is not weakly differentiable, so we cannot apply the SURE method to rank
constrained or rank regularized cases. The soft thresholding function is weakly P differentiable. ′ Remarkably
the divergence expression derived in [34] for any method of the form X̂ = k uk hk (σk ; β)vk turns out to
depend only on the singular values, so it is practical to evaluate. When the singular values of Y are distinct:

SURE{β} = ||X̂ − Y ||2F − M N σ02

 Xmin(M,N ) h(σi ; β) Xmin(M,N ) Xmin(M,N ) σi hi (σi ; β) 
2
+ 2σ0 |M − N | + ḣi (σi ; β) + 2 . (7.41)
i=1 σi i=1 i̸=j σi2 − σj2

Example. See the demo notebook and Fig. 7.15. A practical challenge is that one must know the noise
variance σ02 to apply this method. As seen in Fig. 7.15, the minimizer β∗ of the SURE function is remarkably
close to the minimum MSE choice of β.
© J. Fessler, September 2, 2023, 07:03 (class version) 7.60

Figure 7.15: Left: Similarity of SURE to NRSME. Right: singular values of data, latent matrix, and SURE
estimate.
© J. Fessler, September 2, 2023, 07:03 (class version) 7.61

7.6.2 OptShrink
Both SVST and SVHT require the user to select a regularization parameter β. The OptShrink method,
based on random matrix theory, provides optimal singular value shrinkage for matrix denoising [18], seen
in Fig. 7.16.

Figure 7.16: Singular value shrinkage.

© J. Fessler, September 2, 2023, 07:03 (class version) 7.62

The model here is: Y = X + Z, where Y is the observed matrix, X is the latent matrix, and Z is noise. By
random matrix theory, for quite general noise models (bi-unitarily invariant), the best estimate of X is
Xr
X̂ = ŵk uk vk′ , (7.42)
k=1

where ŵk is the best “weight” for the kth SVD component [18].
The OptShrink algorithm has the following steps.
• Input: Y ∈ FM ×N signal+noise matrix
• Input: r̂ = estimate of rank of latent matrix X
Pmin(M,N )
• Compute SVD: Y = k=1 σk uk vk′
• Form tail singular value matrix:

Σr̂ = Diag σr̂+1 , . . . , σmin(M,N ) ∈ F(M −r̂)×(N −r̂)


(7.43)

• For k = 1, . . . , r̂, compute D(σk , Σr̂ ) and D′ (σk , Σr̂ ) using (7.45) and (7.46), then the optimal weight is

ŵk = −2D(σk , Σr̂ )/D′ (σk , Σr̂ ).

• Form final denoised estimate of latent matrix via

Xr̂
X̂ = ŵk uk vk′ . (7.44)
k=1

For estimates of the MSE of this estimate, see [18, eqn. (17)].
© J. Fessler, September 2, 2023, 07:03 (class version) 7.63

Let S denote a K × L rectangular diagonal matrix with diagonal entries s1 , . . . , smin(K,L) . The two key
expressions in this algorithm are
1 1
trace z(z 2 I − SS ′ )−1 trace z(z 2 I − S ′ S)−1
 
D(z, S) = (7.45)
K L
and
1 1
D′ (z, S) = trace z(z 2 I − SS ′ )−1 · trace −2z 2 (z 2 I − S ′ S)−2 + (z 2 I − S ′ S)−1
 
K L
1 1
+ trace z(z I − S ′ S)−1 · trace −2z 2 (z 2 I − SS ′ )−2 + (z 2 I − SS ′ )−1 .
 2 
(7.46)
L K

These equations are practical only if both M and N are sufficiently small. Here we adapt the approach to
′
allow one of M or N (but not both)
 to be large. To evaluate D(z, S), note that D(z, S) = D(z, S ) so WLOG
S
assume K ≥ L (tall). Then S = L and for z ̸= sk :
0
 2   2 
2 ′ 2 SL 0 z I − SL2 0
z I − SS = z I − =
0 0 0 z2I
(  −1 ) L
2 2 K −L X z
 2 ′ −1 z I − S L 0
=⇒ trace z(z I − SS ) = trace z 2 = + .
0 z I z k=1
z − s2k
2
© J. Fessler, September 2, 2023, 07:03 (class version) 7.64

Similarly, as long as z ̸= sk :
L L
! L
!
 2 ′ −1
X z 1 K −L X z X z
trace z(z I − S S) = , =⇒ D(z, S) = + .
k=1
z − s2k
2 KL z k=1
z − s2k
2
k=1
z − s2k
2

A similar simplification is feasible for D′ (z, S). ( HW )

Combining, one can implement OptShrink to calculate the optimal SVD component weights {ŵk } when at
most one of M or N is large.

[Link] Illustration of OptShrink’s robustness to rank estimate

Example. Fig. 7.18 shows X, a 100 × 30 logo image (aka 2D array aka matrix) with rank(X) = 4, and the
noisy data Y = X + Z, and the denoised matrix X̂ for two rank estimates. See demo:
[Link]
[Link]
Fig. 7.17 shows that the error of X̂ for conventional low-rank approximation increases faster (due to over-
fitting noise) than that of OptShrink. OptShrink finds the “threshold” for singular value shrinkage adaptively
from the data, namely from the tail singular values, using random matrix theory.

To extend optimal singular value shrinkage to heteroscedastic noise, see [44].

© J. Fessler, September 2, 2023, 07:03 (class version) 7.65

Figure 7.17: OptShrink demo plots.

© J. Fessler, September 2, 2023, 07:03 (class version) 7.66

Figure 7.18: OptShrink demo results. The NRMSE increases only slightly when r̂ = 9 ≫ r = 4.
© J. Fessler, September 2, 2023, 07:03 (class version) 7.67

Example. Fig. 7.19 illustrates why singular value shrinkage is needed when Y = X + Z where X ∈ R100,50
is rank 5 with singular values 5,4,3,2,1 and the noise is zero mean with standard deviation σ0 = 0.1. Fig. 7.19
shows histograms of σk (Y ) for 2000 realizations of Z. Noise causes a positive bias: E[σk (Y )] > σk (X), as
discussed in more detail in Ch. 12.

Figure 7.19: Singular value distributions for a signal+noise matrix Y = X + Z, where rank(X) = 5. The
bottom right plot shows all the singular values of the noise-only matrix Z.
© J. Fessler, September 2, 2023, 07:03 (class version) 7.68

(Read)
7.7 Related methods: autoencoders and PCA
This section discusses unsupervised dimensionality reduction methods that are closely related to low-rank
matrix approximation [45].

7.7.1 Relation to autoencoder with linear layers

An autoencoder is a type of artificial neural network (ANN) that is designed to perform dimensionality
reduction [46]. An autoencoder consists of an encoder and a decoder:

encoder decoder
x ∈ FM → → z ∈ FK → → x̂ = ψ(ϕ(x)) ∈ FM , (7.47)
ϕ ψ

where typically K ≪ M , i.e., the dimension of the code or latent variable z is much less than that of a data
point x.
The simplest possible case is where the encoder and decoder are each linear, and hence each can be repre-
sented via a (learnable) matrix:

ϕ(x) = Ex, E ∈ FK×M (wide)

ψ(z) = Dz, D ∈ FM ×K (tall). (7.48)

In this case, the output is x̂ = DEx.

© J. Fessler, September 2, 2023, 07:03 (class version) 7.69

The natural approach to training this ANN given training data x1 , . . . , xN ∈ FM , for N > K, is
N
X
(D̂, Ê) = arg min ∥xn − DExn ∥22
D∈FM ×K , E∈FK×M n=1

||X − DEX|| 2F ,
 
= arg min X ≜ x1 . . . xN . (7.49)
D∈FM ×K , E∈FK×M

The product DEX has at most rank K due to the dimensions of D and E.
Choices for D̂ and Ê such that D̂ ÊX = UK ΣK VK′ provide an optimal rank-K approximation to X, using
an SVD X = U ΣV ′ .
Based on matching the dimensions, a natural choice for D̂ looks to be D̂ = UK , so we start with that choice
and see if we can choose Ê such that ÊX = ΣK VK′ , i.e., for which ÊU ΣV ′ = ΣK VK′ .
 
I
Multiplying on the right by VK yields the goal ÊU Σ K = ΣK .
0
So evidently we want ÊU = IK 0 , which is achieved when Ê = UK′ .
 

Now we verify that a solution to (7.49) is

Ê = UK′ , D̂ = UK , (7.50)

where UK is the first K left singular vectors of U . For that choice of E and D, and assuming K ≤
© J. Fessler, September 2, 2023, 07:03 (class version) 7.70

min(M, N ), the product becomes

D̂ ÊX = UK UK′ |U ΣV ′ ′ ′
 
{z } = UK UK UK U:,(K+1):M ΣV
X  
ΣK 0K×(N −K) ′
VK V:,(K+1):N = UK ΣK VK′ = X̂K ,
  
= UK IK 0K×M
0(M −K)×K Σ(M −K)×(N −K)

namely the optimal rank at most K approximation to X.

So performing low-rank approximation of training data is closely related to learning a (linear) autoencoder
with a single hidden layer of reduced dimension. If the data dimension is too large to use an SVD, then one
can use backpropagation instead [47–49].
Nonlinear autoencoders are generalizations using nonlinear operations such as ϕ(x) = σ(W x + b) where
σ : FM → FM is some nonlinear function, W is a M × K matrix and b ∈ FM is bias vector, and one must
learn W and b from training data [50]. See ( HW ) for the case where the encoder and decoder are affine
instead of linear.
© J. Fessler, September 2, 2023, 07:03 (class version) 7.71
(Read)
7.7.2 Relation to principal component analysis (PCA)
So far we have focused on low-rank approximation of a matrix. Often the low-dimensional subspaces
span(UK ) and span(VK ) are more important to us than the overall approximation ÂK = UK ΣK VK′ , because
those subspaces provide the opportunity for dimensionality reduction. A classic unsupervised method for
dimensionality reduction is principal component analysis (PCA). Here we derive PCA using two perspec-
tives: a subspace learning perspective, and a variation maximization perspective. In both cases,
 we begin with

M
a set of N training data vectors x1 , . . . , xN ∈ F , organized into a M ×N data matrix X ≜ x1 . . . xN .
The subspace learning perspective makes an affine modeling assumption:

xn ≈ µ + Qzn , zn ∈ FK

where µ ∈ FM allows the data to be centered away from 0, and where Q ∈ VK (FM ) is a M ×K semi-unitary
matrix whose columns provide an orthonormal basis for the subspace R(Q). The zn values are coefficients
for vectors in that subspace. (In factor analysis, one also makes statistical assumptions about zn .) The
primary goal is to learn the basis Q, but we also must fit µ and {zn }. A natural approach is the following
least-squares formulation:
N
X
Q̂ = arg min min min ∥µ + Qzn − xn ∥22 . (7.51)
Q∈VK (FM ) z1 ,...,zN ∈FK µ∈FM
n=1
© J. Fessler, September 2, 2023, 07:03 (class version) 7.72

We start by minimizing w.r.t. zn , which is simply a LLS problem, so by §5.5.2:

ẑn = arg min ∥Qzn − (xn − µ)∥22 = Q+ (xn − µ) = Q′ (xn − µ),

zn ∈FK

because Q is semi-unitary. Substituting this solution into (7.51) yields the simplified minimization problem
XN 2
Q̂ = arg min min ∥µ + QQ′ (xn − µ) − xn ∥2
Q∈VK (FM ) µ∈F
M n=1
XN 2
= arg min min PQ⊥ (µ − xn ) 2
. (7.52)
Q∈VK (FM ) µ∈FM n=1

Next we minimize w.r.t. µ by zeroing the gradient, recalling (5.12):

XN XN
0= (PQ⊥ )′ (PQ⊥ )(µ − xn ) =⇒ N PQ⊥ µ̂ = PQ⊥ xn .
n=1 n=1

These normal equations have multiple solutions (when K < M ) because translating the origin of an affine
subspace along
Pthe subspace has no effect. The standard solution choice is simply the sample average of the
data: µ̂ = N1 Nn=1 xn . Substituting this solution into (7.52) yields a simpler subspace learning problem:
XN 2
Q̂ = arg min PQ⊥ yn 2
(7.53)
Q∈VK (FM ) n=1

where yn ≜ xn − µ denotes the de-meaned data.

© J. Fessler, September 2, 2023, 07:03 (class version) 7.73

To simplify further:
XN 2
XN 2 2
PQ⊥ yn 2
= ∥(I − QQ′ )yn ∥2 = ||Y − QQ′ Y ||F ,
n=1 n=1

where Y ≜ y1 . . . yN ∈ FM ×N . Thus (7.53) becomes

 

2
Q̂ = arg min ||Y − QQ′ Y ||F . (7.54)
Q∈VK (FM )

This problem is a special case of (7.49) where we require D = E ′ and D, E ′ ∈ VK (FM ). There we found
that the solution in (7.50) is simply UK , the first K left singular vectors of the data (here Y ), so we conclude
that Q̂ = UK .
An alternative approach is to expand the norm in (7.53):
2
PQ⊥ yn 2
= (PQ⊥ yn )′ (PQ⊥ yn ) = yn′ PQ⊥ yn = yn′ (I − PQ )yn = ∥yn ∥22 − trace{Q′ yn yn′ Q} .
The first term is independent of Q, so (7.53) simplifies to
XN
Q̂ = arg min (− trace{Q′ yn yn′ Q}) = arg max trace{Q′ Y Y ′ Q}, (7.55)
Q∈VK (FM ) n=1 Q∈VK (FM )

where Y Y ′ = N ′ M ×M
P
n=1 yn yn ∈ F is the sample covariance matrix of the de-meaned data Y . This trace
maximization problem must have the same solution, leading to the useful trace maximization property
2
arg max trace{Q′ Y Y ′ Q} = arg min ||Y − QQ′ Y ||F . (7.56)
Q∈VK (FM ) Q∈VK (FM )
© J. Fessler, September 2, 2023, 07:03 (class version) 7.74

In summary, we conclude that, from the SVD Y = U ΣV ′ of the de-meaned data, the subspace estimate is
2
Q̂ = arg min ||Y − QQ′ Y ||F = arg max trace{Q′ Y Y ′ Q} = UK . (7.57)
Q∈VK (FM ) Q∈VK (FM )

We have now solved (7.51) for {zn }, µ and Q. The fitting error (lowest cost) is similar to (7.3):
N
X N
X
min min min ∥µ + Qzn − xn ∥22 = ∥µ̂ + Q̂ẑn − xn ∥22
Q∈VK (FM ) z1 ,...,zN ∈FK µ∈FM
n=1 n=1
XM
= ||PQ̂⊥ Y ||2F = σk2 (Y ). (7.58)
k=K+1
(Read)
[Link] Maximum variance projections
The variance maximizing perspective on PCA focuses on learning a linear transform T ∈ FN ×K (having
orthonormal columns) that maximizes variance when applied to the data via T ′ X. The matrix T ′ will be
wide (K ≪ N ), so this product provides dimensionality reduction. PCA first subtracts from each the mean
of all training vectors to form a de-meaned data matrix:
1
X1 ≜ x1 − µ̂ . . . xN − µ̂ = X − X 1N 1′N ,
 
(7.59)
N
where X and µ̂ were defined above. Note that X1 1N = 0M . Now PCA finds the first column of T by seeking
the (unit norm) linear combination of rows of X (elements of x) that maximizes the (empirical) variance as
© J. Fessler, September 2, 2023, 07:03 (class version) 7.75

follows:
N
1 X ′ 2
T:,1 ≜ arg max |u X1 [:, n]| = arg max ∥(u′ X1 )′ ∥2 = arg max ∥X1′ u∥2 = u1 ,
u : ∥u∥=1 N n=1 u : ∥u∥=1 u : ∥u∥=1

where an SVD of X1 is X1 = U ΣV ′ . Another way of writing this uses a Rayleigh quotient [51, §6.4]:
u′ X1 X1′ u
arg max , (7.60)
u∈FM \{0} ∥u∥22
where we recognize X1 X1′ /N ∈ FM ×M as an empirical covariance matrix.
Having found that u1 is the (unit norm) linear combination that maximizes the (empirical) variance, we can
now remove the component along that direction and then seek another linear combination that maximizes the
variance of what is left. To remove the component along u1 we construct a new matrix as follows:
⊥
X2 ≜ PR(u1)
X1 = (I − PR(u1 ) )X1 = (I − u1 u′1 )X1 . (7.61)
Note that X2 1 = 0 so we still have a zero-mean array. Now PCA finds the second column of T as follows:
T:,2 ≜ arg max ∥(u′ X2 )′ ∥2 = arg max ∥X1′ (I − u1 u′1 )u∥2 = u2 .
u : ∥u∥=1 u : ∥u∥=1

To see why u2 is the solution here, use an SVD of X1 to write

X2 = (I − u1 u′1 ) X1 = (I − u1 u′1 ) U ΣV ′ = U Diag{0, σ2 , σ3 , . . .} V ′ ,
© J. Fessler, September 2, 2023, 07:03 (class version) 7.76

which is almost an SVD of X2 except for a permutation. Clearly the nonzero singular values of X2 are
(σ2 , σ3 , . . . , σr ) and the spectral norm of X2 is σ2 and the corresponding left singular vector, which is its
principal left singular vector, is u2 .
One can continue this reasoning to show that T = UK , the first K columns of U , provides a linear transform
T having orthonormal columns that maximizes the resulting variances of the scores T ′ X1 . We say that UK′
represents the linear combinations that “explain the most variance” in the data.
Using the same SVD of X1 , we have UK′ X1 = UK′ U ΣV ′ = ΣK VK′ , where VK denotes the first K columns
of V . Thus instead of storing all M N elements of X1 it suffices to store the M K and N K elements of UK
and VK (and the K diagonal elements of ΣK ) reflecting dimensionality reduction.
(Read)
[Link] Units and normalization
The requirement that T have orthonormal columns means that the elements of T must be unitless. When T
is unitless, then the matrix product T ′ X is valid only if the rows of X have the same units (or are unitless);
cf. §[Link]. When the features in xn have different units, it is essential to normalize the rows of X as a
preprocessing step, typically by dividing each row by its standard deviation, so that the product T ′ X is valid.

Demo 7.5 illustrates applying PCA to images of hand-written digits. Fig. 7.20 shows 60 images from each
of three hand-written digits (0, 1, 4) that are reduced to K = 2 dimensions. To apply PCA here, we vec each
28 × 27 image into a vector in R756 , and then arrange those vectors into a 756 × 180 matrix. After computing
the SVD, we reshape each of the first two left singular vectors into a 28 × 27 matrix for display in Fig. 7.20.
© J. Fessler, September 2, 2023, 07:03 (class version) 7.77

This approach of forming a matrix where each column is a vectorized version of one image in a set of images
is very common for SVD methods.

Figure 7.20: Results from PCA demo.

(Read)
[Link] PCA extensions
Many nonlinear dimensionality reduction methods exist that generalize PCA, including kernel PCA [52],
nonlinear PCA [53], Laplacian eigenmaps [54, 55] (see §8.9.6), t-distributed stochastic neighbor em-
bedding (t-SNE) [56], locally linear embedding [57] and PaCMAP [58].
Applications of PCA methods are innumerable, including face recognition using eigenfaces [59] and analysis
of macromolecules using cryo-EM [60].
© J. Fessler, September 2, 2023, 07:03 (class version) 7.78

7.8 Subspace learning for classification

This section discusses an application of low-rank matrix approximation to subspace learning. In this setting,
the low-dimensional subspace R(UK ) and its basis UK are more important to us than the overall approxima-
tion ÂK = UK ΣK VK′ .

In supervised classification, we are given N labeled training data

samples xj,1 , . . . , xj,N for each of J classes of objects, where
each feature vector xj,n ∈ FM . We would like to learn some-
thing about the nature of each class so that later when we get a
new sample vector x0 we can assess which class it is most like.
One basic approach to classification is the K-nearest neighbors
method. When K = 1, this method simply chooses the class of
the nearest neighbor to a test point x0 among the training sam-
ples. Fig. 7.21 illustrates this process that is expressed mathemat-
ically as:
ĵ = arg min min ∥x0 − xj,n ∥ . (7.62)
j∈{1,...,J} n=1,...,N

Figure 7.21: Decision boundaries for

nearest-neighbor classification.
© J. Fessler, September 2, 2023, 07:03 (class version) 7.79

This basic approach requires JN comparisons, so complexity grows with training data size. Also, its perfor-
mance may not improve as N increases when the marginal distributions overlap.
Instead, often we try to learn a model from the training data, and then use fit to the model as the basis for
classification. Here we focus on learning a subspace model for each class. Then to classify a test point x0 we
simply compute the distance of x0 to each of the J subspaces and see which is closest, as illustrated in (4.58)
and Fig. 4.12.
To simplify notation, we drop the class subscript j and focus on learning a subspace for a set of samples for
one class type, x1 , . . . , xN where each xn ∈ FM .
Saying that all the xn vectors lie in a subspace of dimension K means that there is some M × K matrix Q
having orthonormal columns such that xn ∈ R(Q). The columns of Q form an orthonormal basis for the
subspace. In practice the data will only lie approximately in a subspace, so xn ≈ Qzn forsome coefficient

vector zn ∈ FK . We want to learn a subspace basis matrix Q from the training data X = x1 . . . xN , a
M × N matrix. Writing xn ≈ Qzn in matrix form:
 
X ≈ Q |{z}
|{z} Z , where Z ≜ z1 . . . zN . (7.63)
|{z}
M ×N M ×K K×N

To find Q and Z we could pursue the following optimization problem on the Stiefel manifold of (4.45):

Q̂ = arg min min ||X − QZ|| F . (7.64)

K×N
Q∈VK (FM ) Z∈F
© J. Fessler, September 2, 2023, 07:03 (class version) 7.80

In this setting, typically K ≪ min(M, N ) so the product QZ is matrix with (at most) rank K.
Thus this problem is essentially a low-rank approximation problem except that here we really care only about
the subspace basis Q and not the coefficients Z (nor their product). The low-rank solution is
K
X
X̂K = σk uk vk′ = UK ΣK VK′ = UK ΣK VK′ . (7.65)
|{z} | {z }
k=1
Q̂ Z

So to learn a subspace basis Q̂ we simply use the first K left singular vectors of the training data X.
After learning an orthonormal basis Qj for a subspace approximation for each of the J classes, we classify
a test point x0 by finding the nearest subspace:

ĵ = arg min min ∥x0 − Qj z∥22 = arg min fj (x0 )

k
j∈{1,...,J} z∈F j∈{1,...,J}

2 2 2
fj (x0 ) ≜ x0 − Qj (Q′j x0 ) 2 = ∥x0 − Pj x0 ∥22 = Pj⊥ x0 2 = PR(U:,(K+1):M ) x0 , (7.66)
| {z } | {z } | {z } | {z }2
(a) (b) (c) (d)

where Pj x ≜ PR(Qj ) x denotes the orthogonal projection of x onto the subspace R(Qj ) and Pj⊥ ≜ I − Pj .
© J. Fessler, September 2, 2023, 07:03 (class version) 7.81

Q7.18 Which form above is the most efficient (per sample) in the usual case where K ≪
min(M, N )?
A: (a) B: (b) C: (c) D: (d) ??

??

Q7.19 For that efficient form, how many multiplies are needed per test sample?
A: JM (2K + 1) B: JM (M + 1) C: J2M K D: JM 2 E: None ??

A discussion task will explore this method with MNIST hand-written digit images.
Learning subspaces via (7.64) is optimal in the Frobenius norm sense of low-rank approximation, but is not
necessarily optimal for the purposes of classification. An extension called supervised PCA aims to find
subspaces that are good for both approximation and classification [61–63].
See [64] for a special issue on subspace learning and a generalization called submanifold learning.
(Read)
7.8.1 Subspace clustering
The subspace learning approach described above, especially (7.64), assumes that the training data is labeled
by class type, so it is supervised learning. There are also unsupervised subspace learning methods that
perform subspace clustering as part of the subspace learning process [65, 66] [14].
Given unlabeled training data x1 , . . . , xN ∈ FM , in a J-subspace clustering approach we want to learn a set
© J. Fessler, September 2, 2023, 07:03 (class version) 7.82

of J orthogonal bases Q1 , . . . , QJ ∈ VK (FM ) ⊂ FM ×K , such that each training point xn is close to one of
the J subspaces {R(Qj )}. In other words, the model here is that all of the data points lie on (or near) a union
of subspaces ∪Jj=1 R(Qj ) . One possible optimization formulation for learning the subspace bases is:

N
2
X
arg min min xn − PQjn xn 2
. (7.67)
Q1 ,...,QJ ∈VK (FM ) j1 ,...,jN ∈{1,...,J} n=1

An alternative way of writing this is to use weights where wnj indicates whether the nth training point is
assigned to the j class:
J
N X
( J
)
2
X X
arg min min wnj xn − PQj xn 2 , W ≜ W ∈ RN ×J : wnj ∈ {0, 1} , wnj = 1 .
Q1 ,...,QJ ∈VK (FM ) W ∈W n=1 j=1 j=1

These are challenging discrete optimization problems called integer programming problems. One can al-
ternate between updating the basis {Qj } and the cluster assignments ({jn } or W ).
An alternative formulation called sparse subspace clustering (SSC) relaxes the discrete problem [65].
A related concept for unlabeled data is called generalized PCA (GPCA) [67].
© J. Fessler, September 2, 2023, 07:03 (class version) 7.83

(Read)
7.9 Subspace tracking and streaming PCA
The methods described thus far in this chapter are suitable when an entire data matrix X is well modeled
as being low rank. This type of “global low rank” model can require a lot of memory when SVD-based
methods are used, and is poorly suited to applications with time-varying data. To reduce memory and/or to
accommodate time-varying data, one can employ adaptive / incremental / online / recursive / streaming /
tracking / updating methods [68, 69]. (The terminology varies in the literature.)
The specific algorithms used depend in part on whether one’s goal is to perform computation incrementally
or to track time-varying subspaces, in analogy with the RLS methods in §5.10. We summarize here one
example of each type. For a thorough survey, see [69].

7.9.1 Incremental SVD

Suppose that data vectors x1 , x2 , . . . arrive sequentially. The incremental
 SVD (ISVD) [70] provides a
method to update an SVD as each new data vector arrives. Let Xn ≜ x1 x2 . . . xn ∈ Fd×n denote
the data matrix at time n, for n ≪ d (high-dimensional data). Suppose we previously computed the SVD of
′
the data matrix at time n − 1: Xn−1 = Un−1 Σn−1 Vn−1 where Un−1 ∈ Fd×d , Σn−1 ∈ Fd×(n−1) , Vn−1 ∈
F(n−1)×(n−1) . After the nth data point arrives, we want to “update” the SVD by finding the SVD of the new
data matrix that has one additional column in it:

Xn = Xn−1 xn = Un Σn Vn′ .
 
(7.68)
© J. Fessler, September 2, 2023, 07:03 (class version) 7.84

The data covariance involves a rank-1 update:

′
Xn Xn′ = Xn−1 xn Xn−1 xn = Xn−1 Xn−1
′
+ xn x′n .
 
(7.69)
Furthermore, defining Dn−1 ≜ Σn−1 Σ′n−1 ∈ Fd×d we have:
′
Un−1 Xn Xn′ Un−1 = Dn−1 + zn zn′ , ′
zn ≜ Un−1 x n ∈ Fd , (7.70)
so, by ( HW ) and [71], the squared singular values of Xn are the (suitably ordered) roots of the rational
function
d
X |[zn ]i |2
1+ , (7.71)
i=1
di − λ
 
w.r.t. λ, where Dn−1 = Diag{di } . To update the left singular vectors, i.e., to find Un = u1 . . . ud , use
the equality Xn Xn′ ui = λi ui , where λi is the square of the ith singular value of Xn . Thus
′
+ xn x′n ui = λi ui


Un−1 Dn−1 Un−1
or equivalently (Dn−1 + zn zn′ ) bi = λi bi where bi ≜ Un−1
′
ui . Rearranging [68] yields (Dn−1 − λi Id ) bi =
′ −1
−(zn bi ) zn , so bi ∝ (Dn−1 − λi Id ) zn . Expressing in terms of ui and normalizing yields
Un−1 (Dn−1 − λi Id )−1 zn
ui = , i = 1, . . . , d. (7.72)
(Dn−1 − λi Id )−1 zn 2
For further refinements that reduce computation and accommodate missing data, see [69, 72].
© J. Fessler, September 2, 2023, 07:03 (class version) 7.85

Explore 7.14. After computing {ui } via (7.72), how would you compute the right singular vectors (if
needed)?

7.9.2 Streaming PCA

Streaming PCA methods are designed for cases where we believe each xn lies in (or near) a K-dimensional
subspace, where that subspace varies (somewhat slowly) with time instead of being static. The goal is to
estimate an orthonormal subspace basis Un ∈ VK (Fd ) for each time point n, from the data acquired up
until that point, i.e., from Xn , where the Un should evolve slowly. One must initialize the process with
some U0 ∈ VK (Fd ), e.g., by applying the Gram-Schmidt process to the first K observations. One popular
approach is Oja’s method [73] that updates the subspace basis recursively as follows:

Un = PVK (Fd ) (Un−1 + ηn xn x′n Un−1 ), (7.73)

where PVK (Fd ) (Y ) denotes the projection of the matrix Y onto the Stiefel manifold VK (Fd ). The step-size
parameter(s) {ηn } control how rapidly the subspace bases can change with time. See [74] for statistical
performance analysis and see [69] for numerical examples with code.

Explore 7.15. Relate the units of the step size ηn to the units of the data vectors {xn }.
© J. Fessler, September 2, 2023, 07:03 (class version) 7.86

7.10 Summary
This chapter described several methods for computing a low-rank approximation to a matrix, all of which
use a SVD. It also discussed methods for choosing the rank or regularization parameter. It applied the
methods to sensor localization (multidimensional scaling) and to subspace learning, which is useful for
classification problems.
There are two main perspectives in this chapter.
• We started by taking an arbitrary matrix A and finding a low-rank approximation ÂK ≈ A.
This perspective has applications in data compression and in dimensionality reduction.
In these applications, we want K to be small, but we also want the approximation error ||ÂK − A|| F to be
small, and there is a trade-off between the two desires as we vary K.
• We then considered matrix denoising applications where we are given a noisy matrix Y that we model as
Y = X + Z where we think X is a low-rank matrix. Then we process Y to make an estimate X̂ of the
latent matrix X.
In this setting, we want the estimation error ||X̂ − X|| F to be small, and we care less about the approxi-
mation error ||X̂ − Y ||F .
All the methods in this chapter require an SVD of an M × N matrix, which is impractical if both M and
N are very large. Extensions for “big data” are an active research area, e.g., using sketching [75, 76],
streaming/online methods [73, 74, 77] and nonconvex optimization [78], e.g., using the Burer-Monteiro
factorization [79–81]; see §10.4.9.
© J. Fessler, September 2, 2023, 07:03 (class version) 7.87

Bibliography

[1] M. Udell and A. Townsend. “Why are big data matrices approximately low rank?” In: SIAM J. Math. of Data Sci. 1.1 (2019), 144–60. DOI:
10.1137/18M1183480 (cit. on p. 7.4).
[2] E. J. Hu, Y. Shen, P. Wallis, Z. Allen-Zhu, Y. Li, S. Wang, L. Wang, and W. Chen. “LoRA: low-rank adaptation of large language models”. In: Proc. Intl.
Conf. on Learning Representations. 2022. URL: [Link] (cit. on p. 7.4).
[3] A. J. Laub. Matrix analysis for scientists and engineers. Soc. Indust. Appl. Math., 2005 (cit. on p. 7.4).
[4] C. Eckart and G. Young. “The approximation of one matrix by another of lower rank”. In: Psychometrika 1.3 (1936), 211–8. DOI:
10.1007/BF02288367 (cit. on pp. 7.5, 7.6, 7.49).
[5] L. Mirsky. “Symmetric gauge functions and unitarily invariant norms”. In: The quarterly journal of mathematics 11.1 (Jan. 1960), 50–9. DOI:
10.1093/qmath/11.1.50 (cit. on pp. 7.6, 7.21, 7.50).
[6] E. Schmidt. “Zur Theorie der linearen und nichtlinearen Integralgleichungen”. In: Math. Ann. 63.4 (1907), 433–76. DOI: 10.1007/BF01449770
(cit. on p. 7.6).
[7] R. B. Cattell. “The scree test for the number of factors”. In: Multivariate Behavioral Research 1.2 (1966), 245–76. DOI:
10.1207/s15327906mbr0102_10 (cit. on p. 7.9).
[8] E. Dobriban. “Permutation methods for factor analysis and PCA”. In: Ann. Stat. 48.5 (Oct. 2020), 2824–47. DOI: 10.1214/19-AOS1907 (cit. on
p. 7.12).
[9] D. Hong, Y. Sheng, and E. Dobriban. Selecting the number of components in PCA via random signflips. 2020. URL:
[Link] (cit. on p. 7.12).
[10] W. Leeb and E. Romanov. “Optimal spectral shrinkage and PCA with heteroscedastic noise”. In: IEEE Trans. Info. Theory 67.5 (2021), 3009–37. DOI:
10.1109/TIT.2021.3055075 (cit. on p. 7.12).
[11] M. Buehrer, K. P. Pruessmann, P. Boesiger, and S. Kozerke. “Array compression for MRI with large coil arrays”. In: Mag. Res. Med. 57.6 (June 2007),
1131–9. DOI: 10.1002/mrm.21237 (cit. on p. 7.14).
[12] T. Zhang, J. M. Pauly, S. S. Vasanawala, and M. Lustig. “Coil compression for accelerated imaging with Cartesian sampling”. In: Mag. Res. Med. 69.2
(Feb. 2013), 571–82. DOI: 10.1002/mrm.24267 (cit. on p. 7.14).
© J. Fessler, September 2, 2023, 07:03 (class version) 7.88

[13] P. Drineas and I. C. F. Ipsen. “Low-rank matrix approximations do not need a singular value gap”. In: SIAM J. Matrix. Anal. Appl. 40.1 (Jan. 2019),
299–319. DOI: 10.1137/18m1163658 (cit. on p. 7.17).
[14] Y-L. Yu and D. Schuurmans. “Rank/norm regularization with closed-form solutions: application to subspace clustering”. In: Proc. 27th Conf.
Uncertainty in AI. 2011, 778–85. URL: [Link] (cit. on pp. 7.21, 7.50, 7.81).
[15] P. Verboon and W. J. Heiser. “Resistant lower rank approximation of matrices by iterative majorization”. In: Comp. Stat. Data Anal. 18.4 (Nov. 1994),
457–67. DOI: 10.1016/0167-9473(94)90163-5 (cit. on p. 7.21).
[16] E. J. Candes, X. Li, Y. Ma, and J. Wright. “Robust principal component analysis?” In: J. Assoc. Comput. Mach. 58.3 (May 2011), 1–37. DOI:
10.1145/1970392.1970395 (cit. on p. 7.21).
[17] T-H. Oh, Y-W. Tai, J-C. Bazin, H. Kim, and I. S. Kweon. “Partial sum minimization of singular values in robust PCA: algorithm and applications”. In:
IEEE Trans. Patt. Anal. Mach. Int. 38.4 (2016), 744–58. DOI: 10.1109/TPAMI.2015.2465956 (cit. on pp. 7.26, 7.29).
[18] R. R. Nadakuditi. “OptShrink: an algorithm for improved low-rank signal matrix denoising by optimal, data-driven singular value shrinkage”. In: IEEE
Trans. Info. Theory 60.5 (May 2014), 3002–18. DOI: 10.1109/TIT.2014.2311661 (cit. on pp. 7.27, 7.61, 7.62).
[19] Y. Shitov. “The nonnegative rank of a matrix: hard problems, easy solutions”. In: SIAM Review 59.4 (2017), 794–800. DOI: 10.1137/16M1080999
(cit. on p. 7.27).
[20] L. T. K. Hien, D. N. Phan, N. Gillis, M. Ahookhosh, and P. Patrinos. “Block Bregman majorization minimization with extrapolation”. In: SIAM J. Math.
of Data Sci. 4.1 (2022), 1–25. DOI: 10.1137/21M1432661 (cit. on p. 7.27).
[21] M. B. Cohen, J. Nelson, and D. P. Woodruff. Optimal approximate matrix product in terms of stable rank. 2016. URL:
[Link] (cit. on p. 7.28).
[22] R. J. Woodham. “Photometric method for determining surface orientation from multiple images”. In: Optical Engineering 19.1 (1980), p. 191139. DOI:
10.1117/12.7972479 (cit. on p. 7.29).
[23] H. Hayakawa. “Photometric stereo under a light source with arbitrary motion”. In: J. Opt. Soc. Am. A 11.11 (Nov. 1994), 3079–89. DOI:
10.1364/JOSAA.11.003079 (cit. on p. 7.29).
[24] R. Basri and D. Jacobs. “Photometric stereo with general, unknown lighting”. In: Proc. IEEE Conf. on Comp. Vision and Pattern Recognition. Vol. 2.
2001, 374–81. DOI: 10.1109/CVPR.2001.990985 (cit. on p. 7.29).
[25] L. Wu, A. Ganesh, B. Shi, Y. Matsushita, Y. Wang, and Y. Ma. “Robust photometric stereo via low-rank matrix completion and recovery”. In: ACCV.
2011, 703–17. DOI: 10.1007/978-3-642-19318-7_55 (cit. on p. 7.29).
© J. Fessler, September 2, 2023, 07:03 (class version) 7.89

[26] S. Ikehata, D. Wipf, Y. Matsushita, and K. Aizawa. “Robust photometric stereo using sparse regression”. In: Proc. IEEE Conf. on Comp. Vision and
Pattern Recognition. 2012, 318–25. DOI: 10.1109/CVPR.2012.6247691 (cit. on pp. 7.29, 7.30).
[27] M. Zhang and M. S. Drew. “Robust surface normal estimation via greedy sparse regression”. In: EURASIP J. Image Vid. Proc. 2015.1 (2015), p. 42.
DOI : 10.1186/s13640-015-0098-x (cit. on p. 7.29).

[28] I. Borg and P. J. F. Groenen. Modern multidimensional scaling: Theory and applications. Springer, 2005. DOI: 10.1007/0-387-28981-X (cit. on
p. 7.37).
[29] E. Peterfreund and M. Gavish. “Multidimensional scaling of noisy high dimensional data”. In: Applied and Computational Harmonic Analysis 51
(2021), 333–73. DOI: 10.1016/[Link].2020.11.006 (cit. on p. 7.42).
[30] E. Waltz. “The algorithm that mapped omicron shows a path forward”. In: IEEE Spectrum (Feb. 2022). URL:
[Link] (cit. on p. 7.42).
[31] J. J. Moreau. “Proximité et dualité dans un espace hilbertien”. In: Bulletin de la Société Mathématique de France 93 (1965), 273–99. URL:
[Link] (cit. on p. 7.43).
[32] N. Parikh and S. Boyd. “Proximal algorithms”. In: Found. Trends in Optimization 1.3 (2013), 123–231. DOI: 10.1561/2400000003 (cit. on p. 7.43).
[33] A. Beck. First-order methods in optimization. Soc. Indust. Appl. Math., 2017. DOI: 10.1137/1.9781611974997 (cit. on p. 7.43).
[34] E. J. Candes, C. A. Sing-Long, and J. D. Trzasko. “Unbiased risk estimates for singular value thresholding and spectral estimators”. In: IEEE Trans. Sig.
Proc. 61.19 (Oct. 2013), 4643–57. DOI: 10.1109/TSP.2013.2270464 (cit. on pp. 7.49, 7.59).
[35] J. von Neumann. “Some matrix-inequalities and metrization of matrix-space”. In: Tomsk. Univ. Rev. 1 (1937), 286–300 (cit. on p. 7.50).
[36] K. Fan and A. J. Hoffman. “Some metric inequalities in the space of matrices”. In: Proc. Amer. Math. Soc. 6.1 (Feb. 1955), 111–6. DOI:
10.2307/2032662 (cit. on p. 7.50).
[37] M. Fazel, H. Hindi, and S. Boyd. “A rank minimization heuristic with application to minimum order system approximation”. In: Proc. Amer. Control
Conf. Vol. 6. 2001, 4734–9. DOI: 10.1109/ACC.2001.945730 (cit. on p. 7.52).
[38] E. J. Candes and T. Tao. “The power of convex relaxation: near-optimal matrix completion”. In: IEEE Trans. Info. Theory 56.5 (May 2010), 2053–80.
DOI : 10.1109/TIT.2010.2044061 (cit. on p. 7.52).

[39] G. Golub, A. Hoffman, and G. Stewart. “A generalization of the Eckart-Young-Mirsky matrix approximation theorem”. In: Linear Algebra and its
Applications 88 (Apr. 1987), 317–27. DOI: 10.1016/0024-3795(87)90114-5 (cit. on p. 7.54).
© J. Fessler, September 2, 2023, 07:03 (class version) 7.90

[40] S. Gu, L. Zhang, W. Zuo, and X. Feng. “Weighted nuclear norm minimization with application to image denoising”. In: Proc. IEEE Conf. on Comp.
Vision and Pattern Recognition. 2014, 2862–9. DOI: 10.1109/CVPR.2014.366 (cit. on p. 7.54).
[41] Y. Xie, S. Gu, Y. Liu, W. Zuo, W. Zhang, and L. Zhang. “Weighted Schatten p-norm minimization for image denoising and background subtraction”. In:
IEEE Trans. Im. Proc. 25.10 (Oct. 2016), 4842–57. DOI: 10.1109/tip.2016.2599290 (cit. on p. 7.54).
[42] X. Liu, X-Y. Jing, G. Tang, F. Wu, and Q. Ge. “Image denoising using weighted nuclear norm minimization with multiple strategies”. In: Signal
Processing 135 (June 2017), 239–52. DOI: 10.1016/[Link].2017.01.003 (cit. on p. 7.54).
[43] T. Blu and F. Luisier. “The SURE-LET approach to image denoising”. In: IEEE Trans. Im. Proc. 16.11 (Nov. 2007), 2778–86. DOI:
10.1109/TIP.2007.906002 (cit. on p. 7.58).
[44] W. E. Leeb. “Matrix denoising for weighted loss functions and heterogeneous signals”. In: SIAM J. Math. of Data Sci. 3.3 (2021), 987–1012. DOI:
10.1137/20M1319577 (cit. on p. 7.64).
[45] J. P. Cunningham and Z. Ghahramani. “Linear dimensionality reduction: Survey, insights, and generalizations”. In: J. Mach. Learning Res. 16 (2015),
2859–900. URL: [Link] (cit. on p. 7.68).
[46] G. E. Hinton and R. R. Salakhutdinov. “Reducing the dimensionality of data with neural networks”. In: Science 313.5786 (July 2006), 504–7. DOI:
10.1126/science.1127647 (cit. on p. 7.68).
[47] T. D. Sanger. “Optimal unsupervised learning in a single-layer linear feedforward neural network”. In: Neural Networks 2.6 (1989), 459–73. DOI:
10.1016/0893-6080(89)90044-0 (cit. on p. 7.70).
[48] P. J. B. Hancock, R. J. Baddeley, and L. S. Smith. “The principal components of natural images”. In: Network: Computation in Neural Systems 3.1
(1992), 61–70. DOI: 10.1088/0954-898X_3_1_008 (cit. on p. 7.70).
[49] I. Davies and S. Eglen. “The principal components of natural images [Re]”. In: ReScience C 6.3 (Oct. 2020), p. 6. DOI: 10.5281/zenodo.4139359
(cit. on p. 7.70).
[50] Y. Bengio, A. Courville, and P. Vincent. “Representation learning: A review and new perspectives”. In: IEEE Trans. Patt. Anal. Mach. Int. 35.8 (Aug.
2013), 1798–828. DOI: 10.1109/TPAMI.2013.50 (cit. on p. 7.70).
[51] L. Eldén. Matrix methods in data mining and pattern recognition. Errata: [Link]
[Link] Soc. Indust. Appl. Math., 2007. DOI: 10.1137/1.9780898718867 (cit. on p. 7.75).
[52] B. Scholkopf, A. Smola, and K-R. Muller. “Nonlinear component analysis as a kernel eigenvalue problem”. In: Neural Computation 10.5 (July 1998),
1299–319. DOI: 10.1162/089976698300017467 (cit. on p. 7.77).
© J. Fessler, September 2, 2023, 07:03 (class version) 7.91

[53] M. Scholz, F. Kaplan, C. L. Guy, J. Kopka, and J. Selbig. “Non-linear PCA: a missing data approach”. In: Bioinformatics 21.20 (2005), 3887–95. DOI:
10.1093/bioinformatics/bti634 (cit. on p. 7.77).
[54] M. Belkin and P. Niyogi. “Laplacian eigenmaps and spectral techniques for embedding and clustering”. In: NeurIPS. Vol. 14. 2001. URL:
[Link] (cit. on
p. 7.77).
[55] M. Belkin and P. Niyogi. “Laplacian eigenmaps for dimensionality reduction and data representation”. In: Neural Computation 15.6 (2003), 1373–96.
DOI : 10.1162/089976603321780317 (cit. on p. 7.77).

[56] L. van der Maaten and G. Hinton. “Visualizing data using t-SNE”. In: J. Mach. Learning Res. 9.86 (2008), 2579–605. URL:
[Link] (cit. on p. 7.77).
[57] S. T. Roweis and . L. K. Saul. “Nonlinear dimensionality reduction by locally linear embedding”. In: Science 290.5500 (2000), 2323–6. DOI:
10.1126/science.290.5500.2323 (cit. on p. 7.77).
[58] Y. Wang, H. Huang, C. Rudin, and Y. Shaposhnik. “Understanding how dimension reduction tools work: an empirical approach to deciphering t-SNE,
UMAP, triMap, and PaCMAP for data visualization”. In: J. Mach. Learning Res. 22.201 (2021), 1–73. URL:
[Link] (cit. on p. 7.77).
[59] N. Muller, Lourenco Magaia, and B. M. Herbst. “Singular value decomposition, eigenfaces, and 3D reconstructions”. In: SIAM Review 46.3 (2004),
518–45. DOI: 10.1137/S0036144501387517 (cit. on p. 7.77).
[60] P. A. Penczek, M. Kimmel, and C. M. T. Spahn. “Identifying conformational states of macromolecules by eigen-analysis of resampled cryo-EM
images”. In: Structure 19.11 (Nov. 2011), 1582–90. DOI: 10.1016/[Link].2011.10.003 (cit. on p. 7.77).
[61] E. Bair, T. Hastie, D. Paul, and R. Tibshirani. “Prediction by supervised principal components”. In: J. Am. Stat. Assoc. 101.473 (2006), 119–37. DOI:
10.1198/016214505000000628 (cit. on p. 7.81).
[62] E. Barshan, A. Ghodsi, Z. Azimifar, and M. Z. Jahromi. “Supervised principal component analysis: Visualization, classification and regression on
subspaces and submanifolds”. In: Pattern Recognition 44.7 (July 2011), 1357–71. DOI: 10.1016/[Link].2010.12.015 (cit. on p. 7.81).
[63] A. Ritchie, C. Scott, L. Balzano, D. Kessler, and C. S. Sripada. “Supervised principal component analysis via manifold optimization”. In: IEEE Data
Science Workshop (DSW). 2019, 6–10. DOI: 10.1109/DSW.2019.8755587 (cit. on p. 7.81).
[64] Y. Ma, P. Niyogi, G. Sapiro, and R. Vidal. “Dimensionality reduction via subspace and submanifold learning [From the guest editors]”. In: IEEE Sig.
Proc. Mag. 28.2 (Mar. 2011), 14–126. DOI: 10.1109/MSP.2010.940005 (cit. on p. 7.81).
© J. Fessler, September 2, 2023, 07:03 (class version) 7.92

[65] R. Vidal. “Subspace clustering”. In: IEEE Sig. Proc. Mag. 28.2 (Mar. 2011), 52–68. DOI: 10.1109/MSP.2010.939739 (cit. on pp. 7.81, 7.82).
[66] J. Lipor, D. Hong, Y. S. Tan, and L. Balzano. “Subspace clustering using ensembles of K-subspaces”. In: Information and Inference 10.1 (Nov. 2020),
73–107. DOI: 10.1093/imaiai/iaaa031 (cit. on p. 7.81).
[67] R. Vidal, Y. Ma, and S. Sastry. “Generalized principal component analysis (GPCA)”. In: IEEE Trans. Patt. Anal. Mach. Int. 27.12 (Dec. 2005), 1945–59.
DOI : 10.1109/TPAMI.2005.244 (cit. on p. 7.82).

[68] P. Comon and G. H. Golub. “Tracking a few extreme singular values and vectors in signal processing”. In: Proc. IEEE 78.8 (Aug. 1990), 1327–43. DOI:
10.1109/5.58320 (cit. on pp. 7.83, 7.84).
[69] L. Balzano, Y. Chi, and Y. M. Lu. “Streaming PCA and subspace tracking: the missing data case”. In: Proc. IEEE 106.8 (Aug. 2018), 1293–310. DOI:
10.1109/jproc.2018.2847041 (cit. on pp. 7.83, 7.84, 7.85).
[70] J. R. Bunch and C. P. Nielsen. “Updating the singular value decomposition”. In: Numer. Math. 31.2 (1978), 111–29. DOI: 10.1007/BF01397471
(cit. on p. 7.83).
[71] G. H. Golub. “Some modified matrix eigenvalue problems”. In: SIAM Review 15.2 (1973), 318–34. DOI: 10.1137/1015032 (cit. on p. 7.84).
[72] M. Brand. “Fast low-rank modifications of the thin singular value decomposition”. In: Linear Algebra and its Applications 415.1 (May 2006), 20–30.
DOI : 10.1016/[Link].2005.07.021 (cit. on p. 7.84).

[73] E. Oja. “Simplified neuron model as a principal component analyzer”. In: J. Math. Bio. 15.3 (Nov. 1982), 267–73. DOI: 10.1007/BF00275687
(cit. on pp. 7.85, 7.86).
[74] C. J. Li, M. Wang, H. Liu, and T. Zhang. “Near-optimal stochastic approximation for online principal component estimation”. In: Mathematical
Programming 167.1 (Jan. 2018), 75–97. DOI: 10.1007/s10107-017-1182-z (cit. on pp. 7.85, 7.86).
[75] N. Halko, P. Martinsson, and J. Tropp. “Finding structure with randomness: probabilistic algorithms for constructing approximate matrix
decompositions”. In: SIAM Review 53.2 (2011), 217–88. DOI: 10.1137/090771806 (cit. on p. 7.86).
[76] J. A. Duersch and M. Gu. “Randomized projection for rank-revealing matrix factorizations and low-rank approximations”. In: SIAM Review 62.3 (2020),
661–82. DOI: 10.1137/20M1335571 (cit. on p. 7.86).
[77] D. Zhang and L. Balzano. “Global convergence of a Grassmannian gradient descent algorithm for subspace estimation”. In: AISTATS. Vol. 51. 2016,
1460–8. URL: [Link] (cit. on p. 7.86).
[78] Y. Chi, Y. M. Lu, and Y. Chen. “Nonconvex optimization meets low-rank matrix factorization: an overview”. In: IEEE Trans. Sig. Proc. 67.20 (Oct.
2019), 5239–69. DOI: 10.1109/tsp.2019.2937282 (cit. on p. 7.86).
© J. Fessler, September 2, 2023, 07:03 (class version) 7.93

[79] S. Burer and R. D. C. Monteiro. “A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization”. In: Mathematical
Programming 95.2 (Feb. 2003), 329–57. DOI: 10.1007/s10107-002-0352-8 (cit. on p. 7.86).
[80] Y. Chen and Y. Chi. “Harnessing structures in big data via guaranteed low-rank matrix estimation: Recent theory and fast algorithms via convex and
nonconvex optimization”. In: IEEE Sig. Proc. Mag. 35.4 (July 2018), 14–31. DOI: 10.1109/msp.2018.2821706 (cit. on p. 7.86).
[81] Z. Li, Y. Luo, and K. Lyu. “Towards resolving the implicit bias of gradient descent for matrix factorization: greedy low-rank learning”. In: Proc. Intl.
Conf. on Learning Representations. 2021. URL: [Link] (cit. on p. 7.86).

Solutions to explorations

Explore 7.1. PUK A = UK UK′ A = UK UK′ U ΣV ′ = UK I 0 ΣV ′ = UK ΣK VK′ = ÂK .

 

Explore 7.2. As a reminder that the output is a vector, not a matrix.

Explore 7.3. It avoids an unnecessary transpose operation.

Explore 7.4. The blue dots do not lie on a line, so rank(A) = 2.

Pr
Explore 7.5. k=K+1 σk

Explore 7.6. K = 5 because A has 5 nonzero elements.

Explore 7.7. K = 2 because A has rank 2.

© J. Fessler, September 2, 2023, 07:03 (class version) 7.94

Explore 7.8. This operator maps M × N matrices into K × K matrices, so it cannot be a projection operation.

Explore 7.9. Only if the two rows of C are linearly dependent, i.e., if y = αx for some α ∈ R, i.e., if
yj = αxj for all j, which is the equation for points along a line through the origin with slope α.

Explore 7.10. The 3 vertices of an equilateral triangle.

Explore 7.11. Distortion in location estimates. One should consider other dimensionality reduction methods
in such situations.

Explore 7.12. P ⊥ SP ⊥ would be real{C ′ C}, because ⟨ci , cj ⟩ in (7.14) becomes real{⟨ci , cj ⟩} per (2.32).

Explore 7.13. In this linear case where y = x + ε and ε ∼ N(0, σ02 I):
2 2
SURE{β} = ∥U U ′ y − y∥2 − dσ02 + 2σ02 trace{U U ′ } = PU⊥ y 2 − (d − 2β)σ02
E[SURE{β}] = trace PU⊥ (xx′ + σ02 Id ) −(d − 2β)σ02 = x′ PU⊥ x + σ02 (d − β) − (d − 2β)σ02


= x′ PU⊥ x + σ02 β
h i h 2
i
2
MSE{β} = E ∥U U ′ y − x∥2 = E PU ε − PU⊥ x 2 = x′ PU⊥ x + E ∥PU ε∥22
 

= x′ PU⊥ x + trace{PU E[εε′ ]} = x′ PU⊥ x + σ02 β

© J. Fessler, September 2, 2023, 07:03 (class version) 7.95

Explore 7.14. Because Xn′ ui = σi vi we simply take vi = Xn′ ui / ∥Xn′ ui ∥2 for i = 1, . . . , r.

Explore 7.15. The formula (7.73) for Oja’s method has unit balance only when all elements of xn have the
same units, and in that case the units of ηn must be the inverse of the square of the units of xn . Thus,
normalizing the data seems desirable.