Machine Learning for OR & FE

Dimension Reduction Techniques

Martin Haugh
Department of Industrial Engineering and Operations Research
Columbia University
Email: [Link]@[Link]

Additional References: Christopher Bishop’s PRML, David Barber’s BRML, Hastie, Tibshirani and
Friedman’s ESL, Hastie, Tibshirani and Wainwright’s SLS and Chiang’s Networked Life
Principal Components Analysis
The Maximum Variance Formulation
The Minimum Reconstruction Error Formulation
The Singular Value Decomposition (SVD )
High-Dimensional Data
Some PCA Applications in Machine Learning & OR/FE
Data Compression: Digit Recognition
Data Compression: Eigen Faces
Reduced-Rank LDA
Financial Modeling
Nearest-Neighbor Classification
Latent Semantic Analysis (LSA)
Kernel PCA
Problems with Kernel PCA
Matrix Completion and Collaborative Filtering
Baseline Estimators
Neighborhood Methods
Matrix Factorization (Revisited)
The Netflix Prize
Reputation Systems & PageRank
Other Applications of PageRank
2 (Section 0)
Appendix: Non-Negative Matrix Factorization
Appendix: Probabilistic Latent Semantic Analysis (PLSA)
Appendix: Conditional PLSA

3 (Section 0)
 Principal Components Analysis
Let x = (x1 , . . . , xd )> denote a d-dimensional random vector with
variance-covariance matrix, Σ.

The goal of PCA is to construct linear combinations

d
X
pi = wij xj , for i = 1, . . . , d
j=1

in such a way that:

(1) The pi ’s are orthogonal so that E[pi pj ] = 0 for i 6= j

(2) The pi ’s are ordered in such a way that:

(i) p1 explains the largest percentage of the total variance

(ii) each pi explains the largest percentage of the total variance that has not
already been explained by p1 , . . . , pi−1 .

4 (Section 1)
 The Eigen Decomposition of a Matrix
In practice it is common to apply PCA to the normalized random variables so
that E[xi ] = 0 and Var(xi ) = 1
– achieved by subtracting the means from the original random variables and
dividing by their standard deviations.
– done to ensure that no one component of x can influence the analysis by
virtue of that component’s measurement units.

Key tool of PCA is the eigen decomposition of a square matrix.

The eigen decomposition implies that any symmetric matrix, A ∈ Rd×d can be
written as
A = Γ ∆ Γ> (1)
where:
(i) ∆ is a diagonal matrix, diag(λ1 , . . . , λd ), of the eigen values of A
- without loss of generality ordered so that λ1 ≥ λ2 ≥ · · · ≥ λd .
(ii) Γ is an orthogonal matrix with ith column of Γ containing ith standardized
eigen-vector, γ i , of A.
- “standardized” means γ >i γi = 1
- orthogonality of Γ implies Γ Γ> = Γ> Γ = Id .
5 (Section 1)
 PCA and the Factor Loadings
Since Σ is symmetric can take A = Σ in (1). The positive semi-definiteness of
Σ implies λi ≥ 0 for all i = 1, . . . , d.

The principal components of x then given by p = (p1 , . . . , pd ) satisfying

p = Γ> x. (2)

Note that:
(a) E[p] = 0 since E[x] = 0
(b) Cov(p) = Γ> Σ Γ = Γ> (Γ ∆ Γ> ) Γ = ∆
– so components of p are uncorrelated as desired
– and Var(pi ) = λi .

The matrix Γ> is called the matrix of factor loadings. We can invert (2) to obtain

x = Γp (3)

– so easy to go back and forth between x and p.

6 (Section 1)
 Explaining The Total Variance
We can measure the ability of the first few principal components to explain the
total variance:
d
X d
X d
X
Var(pi ) = λi = trace(Σ) = Var(xi ). (4)
i=1 i=1 i=1

Pd Pd
If we take i=1 Var(pi ) = i=1 Var(xi ) to measure the total variability then by
(4) can interpret
Pk
i=1 λi
Pd
i=1 λi
as the percentage of total variability explained by first k principal components.

7 (Section 1)
 The Maximum Variance Formulation
Let p1 = γ >
1 x be the 1
st
principal component. Can easily show that γ 1 solves

max Var(a> x)
a
subject to a> a = 1.

More generally, let pi = γ >

i x be the i
th
principal component for i = 1, . . . , d.

Then γ i solves

max Var(a> x)
a
subject to a> a = 1 (5)
a> γ j = 0, j = 1, . . . , i − 1. (6)

So each successive principal component finds the linear combination of the

components of x that has maximal variance subject to the normalization
constraint (5) and the orthogonal constraints (6).

8 (Section 1)
 The Minimum Reconstruction Formulation
We can also obtain the principal components of x by posing the problem as one
of minimizing the sum of squared differences between each sample vector, xi and
it’s representation, x̃i , in some lower M -dimensional sub-space
– here we let x1 , . . . , xn denote (centered) sample observations of x.
We have the following problem:

n X
d
" M
#2
X X
min xj,i − bj,k zk,i
B,Z
i=1 j=1 k=1
2
≡ min kX − BZkF (7)
B,Z

where Bj,k := bj,k is the basis in the M -dimensional subspace and Zk,i := zk,i is
the weighting of the ith sample on the k th element in the basis.

9 (Section 1)
 The Minimum Reconstruction Formulation
Problem (7) is sufficient to identify a unique reconstruction point for each xi .
But it is not sufficient to identify a unique optimal B and Z. Why?
- hence we impose the constraint that B> B = IM .
But rotations still (why?) a problem
- recall U a rotation matrix if U> U = I
Choosing directions of maximal variance resolves this issue
- and then end up with solution where B = first M columns of Γ.

10 (Section 1)
 u1
x2
xn

x
en

x1

Figure 12.2 from Bishop: Principal component analysis seeks a space of lower
dimensionality, known as the principal subspace and denoted by the magenta line,
such that the orthogonal projection of the data points (red dots) onto this
subspace maximizes the variance of the projected points (green dots). An
alternative definition of PCA is based on minimizing the sum-of-squares of the
projection errors, indicated by the blue lines.
 The Singular Value Decomposition (SVD)
The singular value decomposition (SVD) of a d × n matrix X is given by

X = U D V> (8)
where:
1. U is d × d and U> U = Id
2. V is n × n and V> V = In
3. D is an d × n diagonal matrix of the positive singular values
We assume the decomposition has ordered the singular values so that the upper
left diagonal element of D has the largest singular value etc.

Now note that

XX> = U D V> V D> U>

= U D̃ U>

where D̃ := D D> is a d × d diagonal matrix with the n squared singular values

on the diagonal.

12 (Section 1)
 The Singular Value Decomposition (SVD)
We can also use the singular value decomposition of a matrix to determine the
principal components. Why?

In practice we don’t know the true covariance matrix and have to make do with
the sample covariance matrix based on (centered) observations, x1 , . . . , xn .

Let Σ be the sample variance-covariance matrix so that

1
Σ = X X>
n−1
where X = [x1 , . . . , xn ].

The SVD is more general than the eigen decomposition and so PCA is often
performed via an SVD decomposition:
1. The eigen vectors are given by the U matrix of the SVD
2. The eigen values are the squares of the singular values from the SVD
Note also that we can use the SVD to approximate X using only the first M
singular values
X = U D V> ≈ UM DM V> M.

13 (Section 1)
Figure 14.21 from HTF: The best rank-two linear approximation to the
half-sphere data. The right panel shows the projected points with coordinates
given by U2 D2 , the first two principal components of the data.
 High-Dimensional Data
In practice it is often the case that n << d
- e.g. have n = 500 images of d = 1000 × 1000 = 106 pixels
- can cause computational difficulties since computing eigen-decomposition
takes O(d3 ) operations.

But number of non-zero eigen values cannot exceed n. Why?

Can exploit this fact as follows:
1. Let x1 , . . . , xn be the n samples of the d-dimensional x
1 >
2. Recall we can write Σ = n−1 XX where X = [x1 , . . . , xn ].
3. Let E be the d × n matrix of eigen-vectors of Σ with corresponding
eigen-values in the n × n diagonal matrix ∆.

This implies
1
XX> E = E∆
n−1

15 (Section 1)
 High-Dimensional Data
Therefore have
1
X> XX> E = X> E∆. (9)
n−1

If we set Ẽ := X> E then (9) can be written as

1
X> XẼ = Ẽ∆
n−1
1
so that Ẽ is the matrix of eigen vectors for the n × n matrix n−1 X> X
- can be calculated in O(n3 ) operations
1 >
- so much easier to calculate than eigen-vectors of the d × d matrix n−1 XX .

So we simply solve for E by finding Ẽ and then setting (why?)

1
E= XẼ∆−1
n−1
– eigen values are unchanged.
But note O(n2 d) work required to compute X> X so could instead use
‘‘economical” SVD.
16 (Section 1)
 Some PCA Applications in Machine Learning & OR/FE
Will consider several application domains for PCA:
Data compression
- hand-writing
- eigen faces

Reduced-rank LDA
Nearest neighbor classification
Financial modeling
Latent semantic analysis
Missing data and collaborative filtering / recommender systems.

There are many, many more application domains.

17 (Section 2)
Data Compression: Digit Recognition Part I
Figure 14.22 from HTF displays 130 handwritten 3’s from a total of 658
- digitized as 16 × 16 grayscale images so can view them as points in R256 .
Figure 14.23 from HTF plots the first 2 principal components along with
images that are close to the vertices of the grid
- vertices are placed at 5%, 25%, 50%, 75% and 95% quantile points.
First component accounts for “lengthening of the lower tail of the three”.
Second component accounts for “character thickness”.
Equation (14.55) from HTF yields the 2-component model

Figure 14.24 from HTF displays the singular values (from the SVD)
compared with a randomized version of the data.
18 (Section 2)
Figure 14.22 from HTF: A sample of 130 handwritten 3’s shows a variety of
writing styles.
Figure 14.23 from HTF: Left panel: the first two principal components of the
hand- written threes. The circled points are the closest projected images to the
vertices of a grid, defined by the marginal quantiles of the principal components.
Right panel: The images corresponding to the circled points. These show the
nature of the first two principal components.
Figure 14.24 from HTF: The 256 singular values for the digitized threes,
compared to those for a randomized version of the data (each column of X was
scrambled).
 Data Compression: Eigen Faces

Figure 15.5 from Barber: 100 training images. Each image consists of 92 × 112 = 10304 greyscale pixels.
The train data is scaled so that, represented as an image, the components of each image sum to 1. The
average value of each pixel across all images is 9.70 × 10−5 . This is a subset of the 400 images in the full
Olivetti Research Face Database.
22 (Section 2)
Figure 15.6 from Barber: (a): SVD reconstruction of the images in fig(15.5) using a combination of the 49
eigen-images. (b): The eigen-images are found using SVD of the images in fig(15.5) and taking the mean
and 48 eigenvectors with largest corresponding eigenvalue. The images corresponding to the largest
eigenvalues are contained in the first row, and the next 7 in the row below, etc. The root mean square
reconstruction error is 1.121 × 10−5 , a small improvement over PLSA (see fig(15.16)).
 Reduced-Rank LDA
Recall the log-ratio in LDA was given by
P(G = k|X = x) π̂k 1 > −1
log = log − (x − µ̂k ) Σ̂ (x − µ̂k )
P(G = l|X = x) π̂l 2
1 > −1
+ (x − µ̂l ) Σ̂ (x − µ̂l ) (10)
2
We let Σ̂ = U DU > be the eigen decomposition of Σ̂ where U is m × m
orthornormal and D is a diagonal matrix of the positive eigen values of Σ̂.
Then used a change of variable with x∗ := D −1/2 U > x so that Cov (X ∗ ) = I.
Could then rewrite (10) as
P(G = k|X = x) π̂k 1 >
log = log − (x∗ − µ̂∗k ) (x∗ − µ̂∗k )
P(G = l|X = x) π̂l 2
1 >
+ (x∗ − µ̂∗l ) (x∗ − µ̂∗l ) (11)
2
where µ̂∗k and µ̂∗l are the centroid estimates for class k and l respectively under
the new coordinate system.
24 (Section 2)
 Reduced-Rank LDA
The K centroids lie in a K − 1-dimensional subspace, HK−1 say, of Rm so we
can project points to HK−1 and then classify the projected points according to
the nearest centroid, after adjusting for log π̂k ’s.

Then looked for an L-dimensional subspace HL ⊆ HK−1 for L < K − 1 where

HL is chosen to maximize between-class variance relative to within-class
variance.

Particular procedure is:

1. First perform a change of variables x∗i := D −1/2 U > xi for i = 1, . . . , N .
2. Let M ∗ be the m × K matrix of class centroids under these new coordinates.
3. Treat the K centroids as “observations”, compute their cov. matrix, B ∗ .
4. Let B ∗ = V DB V > be its eigen decomposition and let v l be the lth
column of V .
−1/2 >
5. Then Zl := v >l D U x = v> ∗
l x is the l
th
discriminant variable.

So central computational tool for reduced-rank LDA is simply PCA (applied to

the centroids).
25 (Section 2)
 1.5
1
0.5
0
−0.5
−1
−1.5
−2
−5 0 5

Figure 12.7 from Bishop: A comparison of principal component analysis with

Fisher’s linear discriminant for linear dimensionality reduction. Here the data in
two dimensions, belonging to two classes shown in red and blue, is to be
projected onto a single dimension. PCA chooses the direction of maximum
variance, shown by the magenta curve, which leads to strong class overlap,
whereas the Fisher linear discriminant takes account of the class labels and leads
to a projection onto the green curve giving much better class separation.
 PCA and Financial Modeling
There are many applications of PCA in finance:
1. Representing movements in term-structure of interest-rates
– useful for interpretation
– building factor models
– and hedging.
2. Representing movements in futures strips
– useful for interpretation
– building factor models
– and hedging.
3. Building factor models for equities.
4. Scenario generation in risk-management.
5. Estimation of risk measures such as VaR and CVaR.
6. And others ...

27 (Section 2)
 Nearest-Neighbor Classification
In nearest-neighbor classification it can be expensive to compute the distance
between data-points when dimensionality is high.
Can overcome this problem by using PCA to approximate distances.
Let Γ1 contain the first k principal components with k large enough so that

x = Γp ≈ Γ1 p1

provides a good approximation. If xa and xb are any two points in Rd then

>
xa − xb xa − xb (Γpa − Γpb )> (Γpa − Γpb )


=
≈ (Γ1 pa1 − Γ1 pb1 )> (Γ1 pa1 − Γ1 pb1 )
= (pa1 − pb1 )> Γ> a b
1 Γ1 (p1 − p1 )
= (pa1 − pb1 )> (pa1 − pb1 ). (12)

– can be much cheaper to compute (12).

And data is also often noisy so projecting onto a lower dimensional sub-space can
produce superior classification: see Figure 15.7 from Barber.

28 (Section 2)
Figure 15.7 from Barber: Finding the optimal PCA dimension to use for
classifying hand-written digits using nearest neighbours. 400 training examples
are used, and the validation error plotted on 200 further examples. Based on the
validation error, we see that a dimension of 19 is reasonable.
 Latent Semantic Analysis (LSA)
In the document analysis literature, PCA is called latent semantic analysis (LSA).
We assume:
1. There are d words in our dictionary, D
2. And n documents in the corpus.

The j th document could then be represented by xj = (xj1 , . . . , xjd )>

- the so-called bag-of-words representation
- where xji refers to the # occurrences of the ith word in the j th document.
But more common to normalize this number:

Definition. The term-frequency, tfji , is the number of times that the word i
appears in document j divided by the number of words in document j, i.e.
#i,j
tfji := P
i #i,j

where #i,j is the number of times that word i appears in document j.

30 (Section 2)
 The TF-IDF Representation
Definition. The inverse-document-frequency, idfi , is defined as
 
n
idfi := log .
# of documents that contain word i

It is possible to use different definitions of idfi

- as long as rarely occurring words are given more weight when they do occur.

Can now define the TF-IDF representation:

xji := tfji × idfi

A corpus of documents is then represented by

X = [x1 · · · xn ]

– a d × n matrix which is typically very large!

31 (Section 2)
 Latent Semantic Analysis (LSA)
LSA is simply a matter of applying PCA to (the variance-covariance matrix of) X
- with the interpretation that the principal directions, i.e. eigen vectors, now
define topics.

One weakness of LSA is that eigen vectors can have negative components
- but how can a word contribute negatively to a “topic”?

Probabilistic latent semantic analysis (PLSA) overcomes this problem

- based on non-negative matrix factorization.

Topic modeling has been a “hot topic” in the machine learning community.

32 (Section 2)
 Examples 15.4 from Barber
The dictionary, D, contains 10 words:

influenza, flu, headache, nose, temperature,

bed, cat, dog, rabbit, pet

The document corpus contains 2000 documents:

1. some articles discuss ailments, and some of these discuss influenza
2. some articles related to pets
3. some articles are background articles unrelated to ailments

Some of the ailment articles are informal and use the word “flu” whereas others
use the more formal “influenza”.

Each document is therefore represented by a 10-dimensional vector

– xji = 1 if ith word occurs in j th document
– xji = 0 otherwise.

See Figures 15.8 and 15.9 from Barber.

33 (Section 2)
Figure 15.8 from Barber: (Top) Document data for a dictionary containing 10
words and 2000 documents. Black indicates that a word was present in a
document. The data consists of two distinct topics and a random background
topic. The first topic contains two sub-topics which differ only in their usage of
the first two words, ‘influenza’ and ‘flu’. (Bottom) The projections of each
datapoint onto the three principal components.
 Information Retrieval and Example 15.5 from Barber
Consider a large collection of documents from which a dictionary, D, is created.

Given a document, xf , how do we find the “closest” document to it in the

collection?

First need a measure of dissimilarity between documents

- could use d(xf , xi ) := (xf − xi )> (xf − xi ).
And then select the document that solves

min d(xf , xi ).
i

In this case it’s a good idea (why?) to scale each vector so that it has unit length
- leads to the equivalent cosine similarity

s(xf , xi ) := cos(θ)

where θ is the angle between xf and xi .

35 (Section 2)
Figure 15.10 from Barber: (a): Two bag-of-word vectors. The Euclidean
distance between the two is large. (b): Normalised vectors. The Euclidean
distance is now related directly to the angle between the vectors. In this case two
documents which have the same relative frequency of words will both have the
same dissimilarity, even though the number of occurrences of the words is
different.

Problem with bag-of-words representation is that TD matrix will have mainly

zeros
- so differences may be due to noise.

LSA can help (why?) solve this problem

- consider Example 15.4 from Barber.
 Kernel PCA
Given x ∈ Rd , can define x̃ := φ(x) ∈ RD where D > d and possibly D = ∞.
Letting X̃ := [x̃1 · · · x̃n ], the corresponding variance-covariance matrix is
1 >
Σφ := X̃X̃ .
n−1
1 >
Let e be an eigen vector of n−1 X̃X̃ with Lagrange multiplier λ. Then have

1 >
X̃X̃ e = λe.
n−1
>
After pre-multiplying by X̃ obtain
1 >
X̃ X̃f = λf (13)
n−1
>
where f := X̃ e
- an eigen equation that we have seen before!

37 (Section 3)
 Kernel PCA
>
Now note that [X̃ X̃]i,j = φ(xi )> φ(xj ) := k(xi , xj ) := Ki,j .

The kernel trick therefore applies and we can write (13) as

1
Kf = λf. (14)
n−1
1
Having found f by solving (14) we can then recover e = λ(n−1) X̃f.

Corresponding (kernel) principal component then given by

p = x̃> e
1
= x̃> X̃f
λ(n − 1)
n
1 X
= k(x, xl )fl
λ(n − 1)
l=1

where fl is the lth component of f.

38 (Section 3)
 Centering the Data
Have implicitly assumed (where?) that the data, i.e. the φ(xi )’s, have already
been centered.
But even if the xi ’s have been centered in general the φ(xi )’s will not. Why?

Straightforward to check that we can resolve this problem, i.e. center the φ(xi )’s,
by replacing K in (14) with K̂ where
n n
1X 1X
K̂i,j := Ki,j − k(xl , xj ) − k(xi , xd )
n n
l=1 d=1
n
1 X
+ 2 k(xl , xk ). (15)
n
l,k=1

39 (Section 3)
 Problems with Kernel PCA
In regular PCA we can use the principal directions, i.e. eigen vectors
corresponding to the largest eigen values, to construct good approximations to x:

x = Γp ≈ Γ1 p1

Can also do this in kernel PCA but the approximation will lie in the higher
dimensional space RD
- hard to visualize or interpret this approximation
- would prefer to project it back to the original space, Rd
- but this is generally hard since the mapping x → φ(x) is non-linear.
One possible solution is to solve
>
min
0
(φ(x0 ) − z) (φ(x0 ) − z)
x

where z ∈ RD is the approximation to φ(x) given by the kernel PCA

- generally a computationally expensive, non-convex optimization problem.

40 (Section 3)
 x2 v1
φ1

x1
φ2

Figure 12.16 from Bishop: Schematic illustration of kernel PCA. A data set in
the original data space (left-hand plot) is projected by a nonlinear transformation
φ(x) into a feature space (right-hand plot). By performing PCA in the feature
space, we obtain the principal components, of which the first is shown in blue
and is denoted by the vector v1 . The green lines in feature space indicate the
linear projections onto the first principal component, which correspond to
nonlinear projections in the original data space. Note that in general it is not
possible to represent the nonlinear principal component by a vector in x space.
Figure 12.17 from Bishop: Example of kernel PCA, with a Gaussian kernel applied to a synthetic data set in
two dimensions, showing the first eight eigenfunctions along with their eigenvalues. The contours are lines
along which the projection onto the corresponding principal component is constant. Note how the first two
eigenvectors separate the three clusters, the next three eigenvectors split each of the cluster into halves, and
the following three eigenvectors again split the clusters into halves along directions orthogonal to the previous
splits.
 Matrix Completion and Collaborative Filtering

Table 7.2 from Hastie, Tibshirani and Wainwright’s Statistical Learning with Sparsity (2015): Excerpt of
the Netflix movie rating data. The movies are rated from 1 (worst) to 5 (best). A symbol shaded grey
circles represents a missing value: a movie that was not rated by the corresponding customer.

43 (Section 4)
 Matrix Completion and Collaborative Filtering
In many applications we are missing many / most of the elements of a matrix X.

e.g. Consider a d × n movies-ratings matrix with d movies and n users

- then xui represents the rating (on some scale) of user u for movie i
- X will then be very sparse in that very few of the xui ’s will be known.

Goal then is to somehow “fill in” or impute the missing values

- has obvious applications in collaborative filtering and recommender systems
- e.g. the Netflix prize.

Two main approaches to solving this problem:

1. Neighborhood methods
2. Matrix factorization methods

Will discuss each of them but first we establish some baseline estimators.

44 (Section 4)
 Baseline Estimators
A couple of obvious baseline estimators:
1. Set x̂ui = x̄, average of all ratings in the (training) data-set.
2. Introduce biases and set
x̂ui = x̄ + bu + bi (16)
where P
u xui
bi := − x̄ (17)
Mi
P
i xui
bu := − x̄ (18)
Mu
Mi = # users that rated movie i and Mu = # movies rated by user u.
3. Use (16) but choose bi ’s, bu ’s via
!
X 2
X X
min (xui − x̂ui ) + λ b2i + b2u .
bi ,bu
(u,i) i u

(16) seems sensible (why?) so will henceforth work with the residual matrix X̃
X̃u,i := x̃ui := xui − x̂ui .
45 (Section 4)
 Neighborhood Methods
Will continue to use the movie-ratings matrix application to develop ideas.
Neighborhood methods rely on calculating nearest neighbors:
1. Two rows are near neighbors if corresponding users have similar taste in
movies.
2. Two columns are near neighbors if the corresponding movies obtained
similar ratings from users.
To compute nearest neighbors we need a distance / similarity metric.
Can use cosine similarity so

x̃>i x̃j
dij := cos(θij ) = p
||x̃i ||2 ||x̃j ||2
P
x̃ui x̃uj
= qP u P (19)
2 2
u x̃ui u x̃uj

where summation is only over users u that rated both movies i and j.
Can define duv for users u and v similarly.
46 (Section 4)
 Neighborhood Methods
Given a fixed movie i can rank all other movies in descending order of |dij |.
Then top L movies in the ranking are the L nearest neighbors of movie i.
Yields new estimator
P
j∈L dij x̃uj
x̂N
ui = x̄ + bu + bi + P i (20)
j∈Li |dij |
| {z }
collaborative filtering term

where Li denotes the neighborhood of movie i.

Question: Why rank movies according to |dij | rather than dij when defining Li ?

Remarks:
1. Could just as easily use duv ’s rather than dij ’s and adjust (20) appropriately.
2. Still need to choose value of L and other obvious tweaks can be made.
3. Can use (20) to define new residual matrix X̃u,i
- and then start working on predicting missing elements of X̃u,i .
47 (Section 4)
 Matrix Factorization Methods
General matrix factorization problem formulated as

X̂ = arg min ||X − M||2F (21)

M∈Rd×n

subject to
Φ(M) ≤ r (22)
where || · ||F denotes Frobenius norm and Φ(·) a constraint function used to
encourage sparsity or regularization etc.
- many interesting problems obtained by varying Φ.

48 (Section 4)
 PCA Revisited
e.g. Take Φ(M) = rank(M). Then solution to (21) is

X̂ = UDr V> (23)

where UDV> = SVD(X) and Dr is constructed from D by keeping just the first
r singular values.
Really just an equivalent version of (7) once you know the following:
Fact: Any rank r matrix M can be factored as M = BZ where B ∈ Rd×r is
orthogonal and Z ∈ Rr×n .

49 (Section 4)
 Back to Matrix Completion
Let Ω ⊆ {1, . . . , d} × {1, . . . , n} denote observed entries of X. Then one version
of matrix completion problem is
X 2
min (xui − mui ) (24)
rank(M)≤r (u,i)∈Ω

r
!2
X X
≡ min xui − buk zki
B,Z
(u,i)∈Ω k=1

where we used fact from previous slide.

Unfortunately a non-convex problem and solution not known
- but many algorithms available for finding local minima.

Algorithm #1: Find optimal Z in terms of B and then try to minimize over B.

50 (Section 4)
 Matrix Completion: Alternating Least Squares
Algorithm #2 (ALTERNATING LEAST SQUARES):
Select an initial B̂ and then iterate the following two steps until convergence:
(i) Optimize over Z for given B̂. Let Ẑ be the optimal solution.
(ii) Optimize over B for given Ẑ. Let B̂ be the optimal solution.

Guaranteed to converge to a local min but may do so slowly in practice.

Can find principal directions or eigen vectors by performing SVD on final
X ≈ B̂Ẑ.
Should use cross-validation or validation / test set to choose r.

Figure 15.11 from Barber based on this method.

51 (Section 4)
Figure 15.11 from Barber: Top: original data matrix X. Black is missing, white present. The data is
constructed from a set of only 5 basis vectors. Middle : X with missing data (80% sparsity). Bottom :
reconstruction found using svdm.m, SVD for missing data. This problem is essentially easy since, despite
there being many missing elements, the data is indeed constructed from a model for which SVD is
appropriate. Such techniques have application in collaborative filtering and recommender systems where one
wishes to ‘fill in’ missing values in a matrix.
 Matrix Completion: Hard-Impute
Algorithm #3 (HARD-IMPUTE):
1. Use initial guess of missing values to obtain initial complete matrix Xcomp .
2. Compute a rank-r SVD of Xcomp , i.e. solve (23) to obtain X̂comp .
3. Use X̂comp to obtain new guess, Xcomp , of complete matrix.
4. Iterate steps 2 and 3 until convergence (which usually occurs but is not
guaranteed).

As before, should use cross-validation or validation / test set to choose r.

53 (Section 4)
 A Convex Relaxation to Matrix Completion
Alternative formulation is based on the following convex relaxation of (24):
X 2
min (xui − mui ) + λ||M||∗ (25)
M
(u,i)∈Ω

where || · ||∗ is the nuclear norm or trace norm with

||M||∗ := sum of singular values of M.

Algorithm #4 (to follow) easy to implement and solves (25).

But first need some definitions ...

54 (Section 4)
 A Convex Relaxation to Matrix Completion
Definition. Let W = UDV> be the SVD of a rank-r matrix W. Then the
soft-threshold operator Sλ (·) is given by

Sλ (W) := UDλ V>

where
Dλ := diag (s1 − λ)+ , . . . , (sr − λ)+ .
 

where s1 ≥ s2 · · · ≥ sr are the singular values on the diagonal matrix D.

Definition. The projection operators, PΩ : Rd×n → Rd×n and

PΩ⊥ : Rd×n → Rd×n are defined as

xui , if (u, i) ∈ Ω
[PΩ (X)]ui :=
0, otherwise.

and 
 ⊥  xui , if (u, i) 6∈ Ω
PΩ (X) ui :=
0, otherwise.

55 (Section 4)
 Matrix Completion: Soft-Impute
Algorithm #4 (SOFT-IMPUTE):
Initialize Xold = 0, set X to be any completion of the matrix and iterate until
convergence:
1. Compute
X̂λ := Sλ PΩ (X) + PΩ⊥ (Xold )

2. Update Xold := X̂λ .

λ chosen via cross-validation or validation / test set approach.

SOFT-IMPUTE requires full SVD of a potentially large dense matrix in each

iteration. Such a matrix may be too large to even store but note

PΩ (X) + PΩ⊥ (Xold ) = PΩ (X) − PΩ (Xold ) + Xold (26)

| {z } |{z}
sparse low rank
Each component on r.h.s of (26) easy to store (why?) for sufficiently large λ.
Implemented in R in softImpute package.
56 (Section 4)
 Other Results and Formulations
Considerable theoretical work on determining how large |Ω| must be in order to
recover X with “high probability”.
- results can depend on solution method, rank(X), distribution of X in matrix
population, “coherence” / “incoherence” of X etc.

Many other interesting instances (and generalizations) of (21) and (22) including
sparse PCA, robust PCA.
e.g. Decompose X into sum of low-rank and sparse components by solving
1
min ||X − (L + S)||2F + λ1 ||L||∗ + λ2 ||S||1
L ∈ Rd×n 2
S ∈ Rd×n

Chapters 7 and 8 of Statistical Learning with Sparsity by Hastie, Tibshirani and

Wainwright a good (but quite technical) source.

57 (Section 4)
 The Netflix Prize
Define the root mean-squared error (RMSE) to be
v
(xui − x̂ui )2
uX
RMSE := t .
u
C
(u,i)

In 2006 Netflix launched a competition with a $1m prize – the Netflix Prize.
Goal of competition was to develop a recommender system that improved
Cinematch (Netflix’s proprietary system) by at least 10%.
If goal not achieved in 3 years, consolation prize of $50k awarded to best team.
Progress prizes of $50k also awarded each year.
Competition sparked huge interest — over 5k teams and 44k submissions!

Data from period 1999 to 2005 on 480k users and 17,770 movies.
Each movie rated by more than 5,000 users on average.
Each user rated more than 200 movies on average.
58 (Section 4)
 The Netflix Data-Set

Figure 4.1 from Chiang’s Networked Life: 20 Questions and Answers (2012): The Netflix Prize’s four
data sets. The training set and probe set were publicly released, whereas the quiz set and test set were
hidden from the public and known only to Netflix. The probe, quiz, and test sets had similar statistical
properties, but the probe set could be used by each competing team as often as they want, and the quiz set
at most once a day. The final decision was based on comparison of the RMSE on the test set.

Improving Cinematch by 10% meant achieving an RMSE of

0.8563 on the quiz set
0.8572 on the test set.
Winner determined by performance on test set.
59 (Section 4)
 Matrix Factorization Results on the Netflix Data-set

Figure 7.2 from Hastie, Tibshirani and Wainwright’s Statistical Learning with Sparsity (2015): Left:
Root-mean-squared error for the Netflix training and test data for the iterated-SVD (HARD-IMPUTE) and
the convex spectral-regularization algorithm (SOFT-IMPUTE). Each is plotted against the rank of the
solution, an imperfect calibrator for the regularized solution. Right: Test error only, plotted against training
error, for the two methods. The training error captures the amount of fitting that each method performs.
The dotted line represents the baseline “Cinematch” score.
60 (Section 4)
 And the Winner Is ...
Turns out that 10% was an excellent choice by Netflix:
Relatively easy to improve on Cinematch by approx 8%
But much(!) harder to improve by 10%.

Most successful submissions used combinations or ensembles of nearest

neighbors and matrix factorization
- enough to get approx 8% improvement.

But many “tricks” required to get additional 2%.

Winning team’s tricks included temporal effects and implicit feedback
- hard to develop without domain knowledge!

Final leader-board available at at

[Link]
- $1m prize determined by 20 minute time differential in submission of top
two entries!

61 (Section 4)
Source: Matrix Factorization Techniques for Recommender Systems by Koren, Bell and Volinsky, Computer
(2009).
 Reputation Systems
The ever increasing amount of online activity and data has been accompanied by
the decentralization of control.
For example:
1. IMDB and Rotten Tomatoes help us choose what movies to watch.
2. How do we decide what sellers to trust on Ebay?
3. Who should we follow on Twitter?
4. How does Google figure out the importance of web-pages?
5. Why is Yelp so popular?

63 (Section 5)
 Reputation Systems
An important problem is in figuring out who or what pages are the “best”
- clearly this is also a dimension reduction problem!

Here we will consider Google’s PageRank system:

- clearly of enormous importance
- has applications to networks beyond evaluating the importance of web-pages
- can use it to highlight the ongoing challenges with reputation systems.

Question: How do we define reputation? How about this definition?

“You’re great if other great people think you’re great.”

Any problems with this? Yes ... but how do reputations in the real world work?

64 (Section 5)
 Google’s Naive PageRank Algorithm
Consider the web and suppose there are a total of d pages.
Say i → j if page i links to page j and define c(i) := # of out-links from page i.

Let Xt ∈ {1, . . . , d} denote the page that a web-surfer visits at time t.

Initial model assumes Xt a Markov chain with transition matrix Q where
Qi,j := P (Xt+1 = j | Xt = i)
 1
= c(i) , if i → j
0, otherwise.

Would like to find a stationary distribution, µ, of the Markov chain so that

µ = µQ.
Could then use µ(i) to measure the importance of the ith page. Why?
But µ will not be unique unless X is irreducible
- and X will not be irreducible in general because (for example) many
web-pages will not have out-links.
65 (Section 5)
 Naive PageRank

Figure 14.46 from HTF: Page-Rank algorithm: example of a small network

Another (less important) issue with defining web-page reputations in this way is
that it is easily “gamed” by just a couple of web-pages that collude.
See for example “Making Eigenvector-Based Reputation Systems Robust to
Collusion” (2004) by Zhang, Goel, Govindan, Mason and Van Roy.
 Google’s PageRank Algorithm
We resolve the uniqueness problem (and make collusion harder!) by instead
assuming a transition matrix Q̄ where

Q̄i,j := (1 − )Qi,j + >0 (27)
d
for 0 <  < 1. Can write (27) equivalently as
 >
Q̄ := (1 − )Q + 11 (28)
d
where 1 is a d × 1 vector of 1’s.

Note that (28) implies that Q̄ is irreducible and therefore has a unique
stationary distribution, µ, satisfying
µ = µQ̄. (29)
Can interpret the resulting Markov chain as one where with probability  we
choose a new web-page randomly from the entire universe of web-pages and with
probability 1 −  we click randomly on a link on the current page
- the random walk interpretation.
67 (Section 5)
 Google’s PageRank Algorithm
µ gives the page-ranks of all the web-pages.
If we run the Markov chain for a sufficiently long time then proportion of time we
spend on page i is µ(i).
Page-rank or “importance” of a page is therefore captured by importance of the
pages that link to it
- so still a circular definition of reputation!

PageRank has been generalized, manipulated, defended, etc. and reputation

systems in general have been the cause of much litigation
- should give you an idea of just how important they are!

Computing PageRank is computationally intensive. Rather than solving (29)

directly there are two popular approaches:
1. Power iteration – the most commonly used approach.
2. Monte-Carlo which exploits our earlier random walk interpretation
- an approach that has been proposed for the fast updating of PageRank.

68 (Section 5)
 Using Power Iteration to Estimate PageRank
Power iteration is very simple:
1. Choose an initial vector, µ0 . (Easy choice is µ0 = 1/d.)
2. Then iterate µt+1 = µt Q̄ for t = 1, . . . until convergence.
Question: Is convergence guaranteed?
Answer: Yes!
Pd
Proof: Let ∆t := i=1 |µt (i) − µ(i)| where µ is the solution to (29).
Using the triangle inequality it is easy to see that
X |µt (j) − µ(j)|
|µt+1 (i) − µ(i)| ≤ (1 − ) (30)
c(j)
j: j links to i

Summing (30) over i yields ∆t+1 ≤ (1 − )∆t .

Therefore (why?) ∆t → 0 as t → ∞ and so µt (i) converges to µ(i).2

Question: Is the convergence fast?

69 (Section 5)
 PageRank
Google originally used PageRank to index the web and to provide an ordering of
web-pages according to their reputation / importance.

This ordering was vital to search and advertising

- but also provided an obvious motive for companies to try and increase the
PageRank of their web-pages.

Over the years the PageRank system has adapted and heuristics added
- in part to counter attempts at collusion and manipulation.

While the specific details are private it is very likely that PageRank still plays the
key role in Google’s search algorithms.

70 (Section 5)
 PageRank and Search
Question: How might you use PageRank as part of a search system?
Here’s how ...
1. First compute the PageRank of all web-pages and store it in descending
order of importance
- just like a book with the 1st page having the highest PageRank, the 55th
page having the 55th highest PageRank etc.

2. A reverse index for all commonly searched terms is also computed and stored
- this is exactly like an index at the back of a book telling you what pages
important terms appear on.

Steps 1 and 2 are updated periodically.

3. When a user searches for “dog” say, the reverse index tells you exactly what
web-pages the term “dog” appears on and returns the first n of them.
- This is easy and “correct”. Why?

Question: How do you handle a multiple term search, e.g. “dog” and “pony”?

71 (Section 5)
 Other Applications / Extensions of PageRank
The ordering returned by PageRank can easily be adapted to reflect a preference
for a specific subset, C say, of web-pages.
Can do this by changing the algorithm so that w.p.  we choose a new web-page
randomly from C and w.p. 1 −  we click randomly on a link on the current page
- is known as personalized PageRank.

PageRank can be applied in most network settings where links between nodes
indicate some sort of preference. Other examples include:

1. Recommendation systems: form a bipartite graph with movies on one side

and users on the other. Add a link between a movie and a user if the user
liked that movie. Now PageRank gives an ordering of movies and users.
Question: How would you recommend movies to a specific user?

2. Networks can also be created out of users and tweets on Twitter with links
added to identify creators of tweets and followers. Can use such a system to
identify influential people and make recommendations to users.

72 (Section 5)
 Appendix: Non-Negative Matrix Factorization
Non-negative matrix factorization is an alternative to PCA where the matrix
components are required to be non-negative
- useful for modeling non-negative data, e.g. images, text, trading volume etc.

Given a d × n matrix, X ≥ 0, the goal is to find W ≥ 0 and H ≥ 0 such that

X ≈ WH

where W is d × r, H is r × n and r ≤ max(d, n).

One approach is to maximize

d X
X n h i
L(W, H) := xij log (WH)ij − (WH)ij (31)
i=1 j=1

- the log-likelihood from a model in which xij ∼ Poisson ((WH)ij )

73 (Section 6)
 Appendix: Non-Negative Matrix Factorization
L(W, H) is not convex and therefore does not have a unique maximum
- instead must search for a good local maximum.

Can be shown the following updates converge to a local maximum:

Pn
j=1 hkj xij / (WH)ij
wik ← wik Pp (32)
j=1 hkj
Pd
i=1 wik xij / (WH)ij
hkj ← hkj Pd (33)
i=1 wik

Figure 14.33 from HTF shows NNMF, vector quantization and PCA applied to a
database of n = 2, 429 images
- each image represented by a 19 × 19 matrix of pixels so d = 381
- positive values are represented by black pixels
- negative values are represented by red pixels.
Note how NNMF represents faces with a set of basis images each of which
resembles a face part.
74 (Section 6)
Figure 14.33 from HTF: Non-negative matrix factorization (NMF), vector quantization (VQ, equivalent to
k-means clustering) and principal components analysis (PCA) applied to a database of facial images. Details
are given in the text. Unlike VQ and PCA, NMF learns to represent faces with a set of basis images
resembling parts of faces.
Figure 14.34 from HTF: Non-uniqueness of the non-negative matrix
factorization. There are 11 data points in two dimensions. Any choice of the
basis vectors h1 and h2 in the open space between the coordinate axes and data,
gives an exact reconstruction of the data.

Often claimed that the basis in NNMF will be much easier to interpret in problem
applications where it makes sense for basis elements to be non-negative

- but the non-uniqueness of the basis as demonstrated by Figure 14.34 makes

this argument harder to make
- also implies solution to the NNMF problem depends on the starting values.
 Appendix: Probabilistic Latent Semantic Analysis (PLSA)
PLSA is closely related to non-negative matrix factorization
- again the matrix elements must be non-negative
- but there is now a probabilistic interpretation of the matrix components
which yields additional constraints on the factorization.

Again given a d × n matrix, X ≥ 0. Goal is to find W ≥ 0 and H ≥ 0 such that

X ≈ WH

where W is d × r, H is r × n and r ≤ max(d, n).

Suppose we have a count matrix C where Cij = # of times U = i and V = j

where U and V are random variables. Then
Cij
Xij := Pd Pn
i=1 j=1 Cij

may be interpreted as P(U = i, V = j).

77 (Section 6)
 Appendix: Probabilistic Latent Semantic Analysis (PLSA)
May then seek a decomposition where
X
P(U = i, V = j) ≈ P̃(U = i | Z = k) P̃(V = j | Z = k)P̃(Z = k) (34)
| {z }| {z }
k
Wik Hkj

≡ P̃(U = i, V = j)

where Z is a hidden or latent random variable and we’re assuming U and V are
conditionally independent given Z.
Note each column of W now represents a possible value of Z.
There are EM-style algorithms for finding W and H but only convergence to a
local minimum is assured
- see Section 15.6 of Barber for details.

Note also that PLSA does not yield a sequential basis

- i.e. the k = 2 solution, W2 say, will not coincide with the first 2 columns of
the k = 3 solution, W3
- not true of PCA.
78 (Section 6)
Figure 15.14 from Barber: PLSA for the document data in fig(15.8). (a): Hinton diagram for two basis
vectors. (b): Hinton diagram for three basis vectors. (c): The projections for the three basis vectors case.
The solution is quite satisfactory since the first 1000 documents are clearly considered to be from the similar
‘ailments’ topics, the next 500 from some other non-specific ‘background’ topic, and the last 500 from a
separate ‘pet’ topic.

Results of applying PLSA to our example term-document matrix with k = 3:

(b) is the W matrix and (c) is the H matrix.
 Appendix: Conditional PLSA
Sometimes it is more natural to model Xij as P(U = i | V = j).
May then seek a decomposition of the form
X
P(U = i | V = j) ≈ P̃(U = i | Z = k) P̃(Z = k | V = j) (35)
| {z }| {z }
k
Wik Hkj

≡ P̃(U = i, | V = j)

Conditional PLSA assumes that U is independent of V conditional on Z.

Again there are EM-style algorithms for performing conditional PLSA.

80 (Section 6)
 Appendix: Discovering the Basis (E.G. 15.7 in Barber)
Images displayed in Fig. 15.15(a) created from 4 base images in Fig. 15.15(b).
Each base image is positive and scaled so that the sum of the pixels is 1; i.e. for
k = 1, . . . , 4 X
P(U = i | Z = k) = 1
i

where U indexes the pixels and Z indexes the 4 base images.

Training images in (a) then constructed as random convex combinations of the 4
base images.
Each training image has elements P(U = i | V = j)
U indexes the pixels
V indexes the training images.

Goal: given only the training images, reconstruct the 4 basis images.
Results are shown in Figure 15.15(c) and are very good
- but number of elements in basis, i.e. 4, was assumed to be known!
81 (Section 6)
Figure 15.15 from Barber: (a) Training data, consisting of a positive (convex)
combination of the base images. (b): The chosen base images from which the
training data is derived. (c): Basis learned using conditional PLSA on the
training data. This is virtually indistinguishable from the true basis. (d):
Eigenbasis (sometimes called ‘eigenfaces’).

Eigenbasis from the first 4 principal components is given in Figure 15.15(d)

- not nearly as interpretable as the conditional PLSA basis
- but of course we “cheated”.

Figure 15.16 provides a fairer comparison.

Figure 15.16 from Barber: (a): Conditional PLSA reconstruction of the images in Fig 15.5 using a positive
convex combination of the 49 positive base images in (b). The root mean square reconstruction error is
1.391 × 10−5 . The base images tend to be more ‘localised’ than the corresponding eigen-images in Fig
15.6b. Here one sees local structure such as foreheads, chins, etc.
 Appendix: Other Apps of NNMF & (Conditional) PLSA
Two further applications of NNMF and (conditional) PLSA include:
1. Modeling Citations
Let d ∈ {1, . . . , D} index a corpus of research articles and let c ∈ {1, . . . , C}
index citations. Then d and c have the same domain, the index of a research
article.
If document d = i cites article c = j then set Cij = 1. Otherwise set Cij = 0.
We can form a distribution over research articles with
Cij
P(d = i, c = j) = P
C
ij ij

Can now perform PLSA to discover citation-topics. See e.g. 15.9 in Barber.
2. Modeling the Web
Have a collection of web-sites indexed by i.
If web-site j points to web-site i then set Cij = 1. Otherwise set Cij = 0
– yields a directed graph of web-site to web-site links.
A typical web-site will only “discuss” a small number of topics so a PLSA
decomposition might be useful for determining latent topics and authoritative
web-sites.
See Section 15.6.3 of Barber for further details and references.
84 (Section 6)