2 of 130

Eigenvectors and Eigenvalues

Given an n x n matrix A, a non-zero vector v is said to be an eigenvector of A with eigenvalue λ if Av = λv.

In other words, when transformed by matrix A, vector v either shrinks (|λ| < 1) or expands (|λ| > 1) (decreases or increases in magnitude) but its direction remains the same or flips over. Of course, when |λ| = 1, the magnitude of v undergoes no change.

3 of 130

Eigenvectors and Eigenvalues

The preservation of the direction of the vector v when transformed by A gives rise to the name “eigen” – which is German for “self”.

Eigenvectors and eigenvalues come in useful in several applications in physics, image and signal processing and machine learning.

Here we stick to real matrices, though almost all the theory is applicable to complex matrices as well.

5 of 130

Non-uniqueness

Eigenvectors and eigenvalues are not unique.
An n x n matrix A has n eigenvectors (not necessarily all distinct) and n eigenvalues (not necessarily all distinct).
Note that it can happen that two eigenvectors v and w have the same eigenvalue (say) 𝜆 even though v ≠ 𝛼w for any non-zero scalar 𝛼. In such a case any linear combination of v and w (say cv + dw for scalars c, d) is also an eigenvector with eigenvalue 𝜆 because A(cv + dw) = 𝜆(cv + dw).
An eigenvector when scaled by an arbitrary non-zero constant (including a negative one) remains an eigenvector with the same eigenvalue.
By convention, eigenvectors are normalized to unit magnitude, but they contain an inherent ambiguity: if v is an eigenvector of matrix A with eigenvalue 𝜆, then so is -v.

6 of 130

Characteristic equation

The eigenvalues of A can be obtained by finding the roots of the following equation called the characteristic equation:

This equation has n (not necessarily distinct or real) roots – the eigenvalues of A. The equation can be written as follows:

7 of 130

Characteristic equation

The characteristic equation of an n x n matrix A has degree n.
It always has n roots – this is as per a result called the fundamental theorem of algebra.
The roots may not be distinct.
The roots may not be real even if A is real.

8 of 130

Matrix of eigenvectors and eigenvalues

9 of 130

Physical example 1

In this transformation, the red arrow changes direction but the blue arrow does not. The blue arrow is an eigenvector of this shear mapping, and since its length is unchanged its eigenvalue is 1. The red direction in the second figure is also an eigenvector – it is the direction in which the cloth is stretched and it has an eigenvalue greater than 1.

http://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors

10 of 130

Physical example 1

This example helps us design transformation matrices that accomplish a certain task.
Let’s say you wanted to stretch by a factor of 2 in direction [1,1] and apply no stretch or shrinking in direction [1,0].
We treat these directions as eigenvectors with eigenvalues 2 and 1 respectively.
This gives us

11 of 130

Physical example 2

Consider the following matrix which you will often see in computer graphics:

It is a matrix that rotates a 3D point about the Z axis – the Z coordinate remains unchanged, but the X and Y coordinates change by a rotation in the XY plane.

In this case, the Z axis, i.e. vector (0 0 1)^tis an eigenvector with eigenvalue 1.

If the axis of rotation changes, the rotation matrix changes. But the axis of rotation will always be an eigenvector of the rotation matrix with eigenvalue 1.

12 of 130

Properties of Eigenvalues and Eigenvectors

Trace of a matrix = sum of its eigenvalues
Determinant of a matrix = product of its eigenvalues (hence a singular matrix has at least one 0 eigenvalue).
A matrix is invertible iff (= if and only if) all its eigenvalues are non-zero.
The eigenvalues of A^-1 are reciprocals of the eigenvalues of A (can you prove this?).
The eigenvectors of A and A^-1 are the same (can you reason why?).
Every eigenvalue of an orthonormal/unitary matrix has absolute value 1 (why?).
Complex eigenvalues always occur in conjugate pairs, i.e. if a+ib is an eigenvalue of real matrix A for eigenvector v, then a-ib will also be an eigenvalue of A for eigenvector v* (can you prove this?).

13 of 130

Properties of Eigenvalues and Eigenvectors

Given the eigenvalues and eigenvectors of A, what can you say about the eigenvalues and eigenvectors of A-σI for some scalar σ?
Given the eigenvalues and eigenvectors of A, what can you say about the eigenvalues and eigenvectors of A² or Aⁿ for some n > 0?
Note that Aⁿ =A*A*A*…*A (product with itself n-1 times). This is different from A.^n for which A.^n(i,j) = A(i,j)^n.

14 of 130

Properties of Eigenvalues and Eigenvectors

Trace of a matrix (size n x n) = sum of its eigenvalues
Determinant of a matrix = product of its eigenvalues

The product (𝜆-a₁₁)(𝜆-a₂₂)...(𝜆-a_nn) has a term 𝜆^n-1 with the coefficient -(a₁₁ + a₂₂+....+a_nn)_.

15 of 130

Eigenvectors and eigenvalues of special matrices

Diagonal matrix: eigenvalues are elements along the diagonal, eigenvectors are the column vectors normalized to unit magnitude.
Symmetric matrix: eigenvalues are real (proof), eigenvectors corresponding to distinct eigenvalues are orthogonal (proof).

The converse is also true, i.e. a matrix with orthonormal eigenvectors and real eigenvalues is symmetric.

16 of 130

Computation of Eigenvectors and Eigenvalues

The following MATLAB commands compute eigenvectors and eigenvalues:

[V,D] = eig(A) returns two outputs. D is a diagonal matrix containing the eigenvalues. V is a matrix whose columns are the corresponding right eigenvectors.

[V,D] = eigs(A) returns a diagonal matrix D of A's six largest magnitude eigenvalues and a matrix V whose columns are the corresponding eigenvectors.

[V,D] = eigs(A,k) returns a diagonal matrix D of A's k largest magnitude eigenvalues and a matrix V whose columns are the corresponding eigenvectors.

17 of 130

Computation of Eigenvectors and Eigenvalues

Eigenvalues can be computed by computing the roots of the characteristic equation.
This is easy for 2 x 2 matrices, but not for larger matrices, since polynomials in arbitrary degree (greater than 5) usually have no closed form solutions.
How do you compute eigenvectors and eigenvalues?

18 of 130

Computing eigenvectors

Assume n x n matrix A has a unique eigenvalue λ₁ of largest magnitude. Let the corresponding eigenvector be v₁.
Consider vector x₀ which is expressible as a linear combination of eigenvectors of A if we assume A has linearly independent eigenvectors, i.e.

Repeated multiplication of A with x₀ converges to a multiple of v₁. Why?

19 of 130

Each such term tends to 0, as k tends to ∞, because λ₁ is the largest magnitude eigenvalue of A.

This method of eigenvector computation is called as power iteration. It converges to the eigenvector with the largest eigenvalue.

Caution: The possibility that the starting vector x₀ has no component in the direction v₁ (i.e. if α₁ = 0) is very small, and can be eliminated by slightly and randomly perturbing the starting vector.

Caution: If there are two or more eigenvalues that are the largest and have the same magnitude, the starting vector will converge to some linear combination of the corresponding eigenvectors.

Caution: If the initial vector is real and A is real, it will not converge to a complex eigenvector!

20 of 130

Computation of Eigenvectors and Eigenvalues

How do you compute the largest eigenvalue given the eigenvector?

What if you wanted to compute the eigenvector corresponding to the smallest eigenvalue? Work with A^-1 instead of A.
What about other eigenvectors? There exist more advanced algorithms for those, which we skip in this course. But we handle one special case.

21 of 130

Computation of Eigenvectors and Eigenvalues

Let A be a symmetric matrix. So its eigenvalues are real and its eigenvectors are orthonormal.

To compute the second eigenvector, work with the following matrix instead of A

where v₁ and λ₁ are the estimated first eigenvector and first eigenvalue.

The eigenvector of A^(-v1)corresponding to its largest eigenvalue is the eigenvector of A corresponding to its second largest eigenvalue.

22 of 130

Computation of Eigenvectors and Eigenvalues

Several more advanced algorithms exist.
The algorithms have ready implementation in a famous library called LAPACK – written in Fortran 90, with C/C++/MATLAB interfaces.

23 of 130

Singular Value Decomposition (SVD)

24 of 130

Singular value Decomposition

For any m x n matrix A, the following decomposition always exists:

Diagonal matrix with non-negative entries on the diagonal – called singular values.

25 of 130

Singular value Decomposition

For any m x n real matrix A, the SVD consists of matrices U,S,V which are always real – this is unlike eigenvectors and eigenvalues of A which may be complex even if A is real.

The singular values are always non-negative, even thought the eigenvalues may be negative.

While writing the SVD, the following convention is assumed, and the left and right singular vectors are also arranged accordingly:

26 of 130

Singular value Decomposition

If only r < min(m,n) singular values are non-zero, the SVD can be represented in reduced form as follows:

27 of 130

Singular value Decomposition

This m by n matrix u_iv^T_iis the product of a column vector u_i and the transpose of column vector v_i. It has rank 1. Thus A is a weighted summation of r rank-1 matrices.

Note: u_i and v_i are the i-th column of matrix U and V respectively.

28 of 130

Singular value decomposition

29 of 130

Application: SVD of Natural Images

An image is a 2D array – each entry contains a grayscale value. The image can be treated as a matrix.
It has been observed that for many image matrices, the singular values undergo rapid decay (note: they are always non-negative).
An image can be approximated with the k largest singular values and their corresponding singular vectors:

30 of 130

Left to right, top to bottom:

Reconstructed image using the first i= 1,2,3,5,10,25,50,100,150 singular values and singular vectors.

Last image: original

31 of 130

Left to right, top to bottom, we display:

where i = 1,2,3,5,10,25,50,100,150.

Note each image is independently re-scaled to the 0-1 range for display purpose.

Note: the spatial frequencies increase as the singular values decrease

32 of 130

SVD: Use in Image Compression

Instead of storing mn intensity values, we store (n+m+1)r intensity values where r is the number of stored singular values (or singular vectors). The remaining m-r singular values (and hence their singular vectors) are effectively set to 0.

This is called as storing a low-rank (rank r) approximation for an image.

33 of 130

Properties of SVD: Best low-rank reconstruction

SVD gives us the best possible rank-r approximation to any matrix (it may or may not be a natural image matrix).
In other words, the solution to the following optimization problem:

is given using the SVD of A as follows:

Note: We are using the singular vectors corresponding to the r largest singular values.

This property of the SVD is called the Eckart Young Theorem.

34 of 130

Properties of SVD: Best low-rank reconstruction

Frobenius norm of the matrix (fancy way of saying you square all matrix values, add them up, and then take the square root!)

Why?

35 of 130

Geometric interpretation: Eckart-Young theorem

The best linear approximation to an ellipse is its longest axis.
The best 2D approximation to an ellipsoid in 3D is the ellipse spanned by the longest and second-longest axes.
And so on!

36 of 130

Properties of SVD: Singularity

A square matrix A is non-singular (i.e. invertible or full-rank) if and only if all its singular values are non-zero.
The ratio σ₁/σ_ntells you how close A is to being singular. This ratio is called condition number of A. The larger the condition number, the closer the matrix is to being singular.

37 of 130

Properties of SVD: Rank, Inverse, Determinant

The rank of a rectangular matrix A is equal to the number of non-zero singular values. Note that rank(A) = rank(S).
SVD can be used to compute inverse of a square matrix:

Absolute value of the determinant of square matrix A is equal to the product of its singular values.

38 of 130

Properties of SVD: Pseudo-inverse

SVD can be used to compute pseudo-inverse of a rectangular matrix:

39 of 130

Properties of SVD: Frobenius norm

The Frobenius norm of a matrix is equal to the square-root of the sum of the squares of its singular values:

40 of 130

Solution to equations of the form � Av = 0

Remember A is a m by n matrix.
The equation Av = 0 may have a solution apart from the zero-vector (also called as trivial solution).
To obtain the solution, we first compute the SVD of A giving A = USV^T. The solution is proportional to the column of V corresponding to the zero singular value.

41 of 130

Solution to equations of the form � Av = 0

To obtain the solution, we first compute the SVD of A giving A = USV^T. The solution is proportional to the column of V corresponding to the zero singular value.

42 of 130

Solution to equations of the form � Av = 0

The derivation on the previous slide shows that a right singular vector corresponding to a zero singular value lies in the nullspace of A.
Now, we will prove that if a vector v lies in nullspace(A), it is a linear combination of right singular vectors of A corresponding to zero singular values. Let A = U𝛴V^T be the SVD of A.

Thus v is a linear combination of right singular vectors corresponding to zero singular values.

We can express v = Vc where c is a vector of n coefficients since V is an orthonormal basis.

43 of 130

Solution to equations of the form � Av = 0

If more than one singular value is 0, then the corresponding right singular vector (and its scaled versions) are also solutions.
In fact, all such right singular vectors with zero singular values span the null-space of A.
Equations of this form arise in several applications – for example, in computer vision (camera calibration, image alignment).
In some cases, Av = 0 has no exact solution. Instead, we seek to find a solution w such that Aw is as small as possible in magnitude. In this case, w is given by the right singular vector of A corresponding to the smallest singular value.

44 of 130

Orthogonal Procrustes problem

Read the notes from:

http://www.cse.iitb.ac.in/~ajitvr/CS740_Fall2026/procrustes.pdf

45 of 130

Geometric interpretation of the SVD

Any m x n matrix A transforms a sphere Q of unit radius (called as unit sphere) in Rⁿ into a hyperellipse in R^m (assume m >= n).

46 of 130

Geometric interpretation of the SVD

Assume A has full rank for now.
The singular values of A are the lengths of the n principal semi-axes of the hyperellipse. The lengths are thus σ₁, σ₂, …, σ_n.
The n left singular vectors of A are the directions u₁, u₂, …, u_n (all unit-vectors) aligned with the n semi-axes of the hyperellipse.
The n right singular vectors of A are the directions v₁, v₂, …, v_n (all unit-vectors) in sphere S, which the matrix A transforms into the semi-axes of the hyperellipse, i.e.

47 of 130

Geometric interpretation of the SVD

Expanding on the previous equations, we get the reduced form of the SVD

n x n diagonal matrix - S

m x n matrix with orthonormal columns - U

n x n orthonormal matrix V

48 of 130

Geometric interpretation of the SVD

But why does A transform the sphere into a hyperellipse?
This is because A = USV^T.
V^T transforms the sphere into another (rotated/reflected) sphere.
S stretches the sphere into a hyperellipse whose semi-axes coincide with the coordinate axes as per V.
U rotates/reflects the hyperellipse without affecting its shape.
As any matrix A has an SVD decomposition, it will always transform the sphere into a hyperellipse.
If A does not not have full rank, then some of the semi-axes of the hyper-ellipse will have length 0!

49 of 130

Computation of the SVD

We will not explore algorithms to compute the SVD of a matrix, in this course.
SVD routines exist in the LAPACK library and are interfaced through the following MATLAB functions:

s = svd(X) returns a vector of singular values.

[U,S,V] = svd(X) produces a diagonal matrix S of the same dimension as X, with nonnegative diagonal elements in decreasing order, and unitary matrices U and V so that X = U*S*V'.

[U,S,V] = svd(X,0) produces the "economy size" decomposition. If X is m-by-n with m > n, then svd computes only the first n columns of U and S is n-by-n.

[U,S,V] = svd(X,'econ') also produces the "economy size" decomposition. If X is m-by-n with m >= n, it is equivalent to svd(X,0). For m < n, only the first m columns of V are computed and S is m-by-m.

s = svds(A,k) computes the k largest singular values and associated singular vectors of matrix A.

50 of 130

SVD Uniqueness

If the singular values of a matrix are all distinct, the SVD is unique – up to a multiplication of the corresponding columns of U and V by a sign factor (or a phase factor e^iθ for complex matrices).
Why?

51 of 130

SVD Uniqueness

A matrix is said to have degenerate singular values, if it has the same singular value for 2 or more pairs of left and right singular vectors.
In such a case any normalized linear combination of the left (right) singular vectors is a valid left (right) singular vector for that singular value.

52 of 130

Any other applications of SVD?

Face recognition – the popular eigenfaces algorithm.
Point matching: Consider two sets of points, such that one point set is obtained by an unknown rotation of the other. Determine the rotation!
Structure from motion: given a sequence of images of a object undergoing rotational motion, determine the 3D shape of the object as well as the rotation at every time instant!

53 of 130

Structure from Motion

CS 763

Ajit Rajwade

54 of 130

Problem definition

Structure from motion refers to the inference of the object’s 3D structure or shape (i.e. the X,Y,Z coordinates of several points on the object’s surface) given a sequence of the object’s images when the object is in relative motion w.r.t. a camera.

55 of 130

Human perception of motion

We humans have the ability to do this inference – see below:

https://www.youtube.com/watch?v=zdKX7Xo3Cb8&feature=player_detailpage#t=270

56 of 130

Contents of the lecture

We are going to study an interesting algorithm called as factorization.
It was developed by Tomasi and Kanade and was published in the early 90s.
The algorithm is simple and elegant.

Tomasi and Kanade, “Shape and motion from image streams under orthography: a factorization method”, International Journal of Computer Vision, 1992.�http://link.springer.com/article/10.1007%2FBF00129684#page-1

57 of 130

Algorithm input and assumptions

Given: A sequence of some F ≥ 3 images of a non-planar object acquired under an orthographic camera moving relative to the object.
The camera may be actually moving and the object could be still, or vice-versa, or both could be in motion.
For simplicity but without loss of generality, we will assume the former.

58 of 130

Algorithm input and assumptions

Let the object consist of n ≥ 3 non-coplanar points – P₁, P₂, …, P_n measured in some world coordinate system.
We will assume that (1) these n 3D points are visible in each of the F frames, and (2) the corresponding n image points are tracked and marked out in each of the F frames.

59 of 130

Algorithm: input and assumptions

We are assuming an orthographic camera.
We are assuming that the whole video sequence is obtained a priori with points tracked.

60 of 130

Algorithm: input and assumptions

Let p_ij = (x_ij, y_ij) = j-th image point (j = 1 to n) in the i-th frame (i = 1 to F).
Assemble matrix W (size 2F x n) as follows:

61 of 130

Algorithm: input and assumptions

Consider the following matrix (size 2F x n) as follows:

For each frame, compute the centroid of the 2D points. Deduct the centroid from the points in every frame to create the new matrix on the left. Why do we do this? We will see soon.

We will prove that this matrix actually has rank at the most 3 under ideal conditions (no noise in point coordinates). This is called the Rank Theorem.

62 of 130

Proof: Rank Theorem

The 3D object is stationary and the camera is moving (performing rotation and translation).
Each time the camera moves, its extrinsic parameters change, i.e. the rotation transformation between the camera axes and the world coordinate axes changes, and also the translation vector between the origin of the camera coordinate system and the origin of the world coordinate system changes.
World coordinate system = It is a coordinate system chosen by the user (say the designer of an environment for a robot to travel, the architect of a building, etc).
The world coordinate system will typically be different from the camera coordinate system.

63 of 130

Proof: Rank Theorem

In the i-th frame, let the translation vector be given as t_i. Let the axes of the camera as measured in the world coordinate system be given as i_i, j_i, k_i = i_i x j_i.

The image coordinates are thus given as follows:

point coordinates measured in camera’s coordinate system

The X and Y axes of the coordinate system of the camera in the i-th frame – represented as vectors in 3D. These axes are expressed in the world coordinate system (if they were expressed in the camera’s coordinate system, they would be (1,0,0) and (0,1,0).

P_j = (X_j,Y_j,Z_j) = coordinates of the j-th point in world coordinate system. {P_j} for j=1 to n constitute the object’s “structure”

64 of 130

Proof: Rank Theorem

Without loss of generality, we will assume that the origin of the world coordinate system is at the centroid of the 3D object.
In other words, we have

65 of 130

Proof: Rank Theorem

Now consider the following equations:

Combining them, we have

66 of 130

Proof: Rank Theorem

Reconsider the following equations:

R has size 2F x 3 and has rank 3 as F ≥ 3.

S has size 3 x n

and will have rank 3 if the points in S are non-coplanar.

So W̃ has rank 3.

67 of 130

What does the rank theorem tell you?

Given the matrix W̃, we compute its SVD as follows:

For i = 1 to F, the i^th and (F+i)^th rows of R give you the vectors i_i and j_i respectively. Since k_i = i_i x j_i, we have axes of the camera coordinate system in the i-th frame. Comparing the camera coordinate systems across consecutive frames tells you how much the camera rotated from one frame to another.
The columns of S give the 3D point coordinates.

Reduced form of the SVD as W̃ has rank only 3

68 of 130

Problem!

But the obtained R and S are not unique because for any invertible 3 x 3 matrix Q, we have:

How do we resolve this?

69 of 130

Problem solution

But the obtained R and S are not unique because for any invertible 3 x 3 matrix Q, we have:
Observe that the rows of a rotation matrix (here RQ) must have unit magnitude. Any two rows must be perpendicular to each other. So we solve for Q by observing that:

These 3 equations are true for all i from 1 to F (i.e. for each frame). These equations are called the metric properties or metric constraints on R. Recall that i_i and j_i are obtained from the R matrix that you get from the SVD of W̃.

70 of 130

Problem solution: not so soon!

We can solve for Q which will satisfy the following equations using Newton’s method (details later):

The final R and S matrices will be as follows:

But these solutions are also unique only up to some unknown orthonormal transformation R₀, i.e.

These 3 equations are true for all i from 1 to F (i.e. for each frame). Recall that i_i and j_i are obtained from the R matrix that you get from the SVD of W̃.

71 of 130

Problem solution: not bad after all!

Note that this R₀ cannot be uniquely obtained by exploiting the metric properties unlike the case of Q (why?).
All this means is that the if you assumed all the camera positions were rotated by some fixed R₀ in every frame, the object coordinates would rotate by a fixed (R₀)^-1 in every frame.
This can be resolved by assuming that in the first frame, the world coordinate system is aligned with the camera coordinate system.

72 of 130

What about camera translation from frame to frame?

This is orthographic projection: so we can never determine the Z component of the translation vector in any frame.
The X and Y components of the translation vector (in frame t) are obtained by the difference between the image centroids in frame t and those in frame t-1.

73 of 130

What about camera rotation from frame to frame?

Compare the i,j,k axes of the camera in frame t and frame t-1.

74 of 130

Measurement noise

The rank theorem says that W̃ has rank 3.

But that is true only when there is no noise in measuring the coordinates of the tracked points in every frame.

What if there is noise? One can attempt to “filter out” the noise in W̃ by considering its rank 3 approximation.

75 of 130

Measurement noise

Consider the SVD:

Due to noise, the rank exceeds 3, but we can create a rank-3 approximation by considering only the 3 largest singular values in D (and their corresponding columns in U and V).

This is the best rank-3 approximation to W̃ as per the well-known Eckart-Young Theorem on SVD.

76 of 130

How to estimate Q?

Look at the following equations (totally 3F in number):

This is a system of non-linear equations, the variables being the 9 entries of Q which we rearrange to yield vector q. We will label each equation as f_k(q) = 0 (k = 1 to 3F).
No closed-form solution unlike linear case ☹

77 of 130

How to estimate Q?

Start with an initial guess for q, for example q_t = vectorized form of identity matrix.
If q_t is the true solution, then f_k(q_t) = 0 for all k from 1 to 3F, and you stop (this won’t happen in the first step when t = 0!).
Instead we want to find vector δ such that f_k(q_t + δ) = 0 for all k.
We seek to find δ by approximating each f_k as a linear function in the neighborhood of q_t.

78 of 130

How to estimate Q?

5. The linear approximation is given as:

6. But we want f_k(q_t+ δ) = 0 for all k. Hence for a given k, we have

This is a 9 x 1 vector of first derivatives. Remember that δ is a 9 x 1 vector.

79 of 130

How to estimate Q?

7. Collecting together 3F such equations, we have:

8. One can solve for δ by pseudo-inverse.

9. But this solution will not exactly satisfy all the equations as we performed a linear approximation which was not fully accurate, and also because a least squares solution for δ is not guaranteed to yield f_k(q_t+ δ) = 0 for all k.

The yellow box contains a 9 x 3F matrix called the Jacobian. Again, δ is a 9 x 1 vector and f(q_t) is also a 9 x 1 vector.

80 of 130

How to estimate Q?

10. Hence we update our solution from q_t to q_t+1 = q_t+ δ.

We repeat the previous steps with t = 0, 1, 2,… and so on until we reach a time when f_k(q_t + δ) ≈ 0 for all k.
This overall method is called Newton-Raphson method of root-finding.

84 of 130

Principal Components Analysis

PCA has many applications – apart from face/object recognition – in image processing/computer vision, statistics, econometrics, finance, agriculture, and you name it!

Why PCA? What’s special about PCA? See the next slides!

85 of 130

PCA: what does it do?

It finds ‘k’ perpendicular directions (all passing through the mean vector) such that the original data are approximated as accurately as possible when projected onto these ‘k’ directions.

We will see soon why these ‘k’ directions are eigenvectors of the covariance matrix of the data!

86 of 130

PCA

Look at this scatter-plot of points in 2D. The points are highly spread out in the direction of the slanted red line.

87 of 130

PCA

This is how the data would look if they were rotated in such a way that the major axis of the ellipse (the red line) now coincided with the Y axis.
As the spread of the X coordinates is now relatively insignificant (observe the axes!), we can approximate the rotated data points by their projections onto the Y-axis (i.e. their Y coordinates alone!). This was not possible prior to rotation!

88 of 130

PCA

As we could ignore the X-coordinates of the points post rotation and represent them just by the Y-coordinates, we have performed some sort of lossy data compression – or dimensionality reduction.

The job of PCA is to perform such a rotation as shown on the previous two slides!

89 of 130

PCA

Aim of PCA: Find the unit vector u passing through the sample mean (i.e. ), such that the projection of any mean-deducted point

onto u, most accurately approximates it.

x_i

z_i

Projection error for y_i

Mean deducted point

Projection of y_i onto u

90 of 130

Total reconstruction error J(u)

The first term does not depend on u, so minimizing J(u) is equivalent to maximizing the following:

91 of 130

Scatter matrix (proportional to the sample covariance matrix)

As S is symmetric, V is orthonormal and eigenvalues in 𝛬 are non-negative

Express u as a linear combination of the eigenvectors – this is possible as V is an orthonormal basis. That is any vector in R^d can be expressed as a linear combination of the columns of V

92 of 130

The quantity u^TSu is maximized if a₁ = 1 and the other a_j values are all zero.
Note that 𝜆₁ ⩾ 𝜆₂ ⩾ … ⩾ 𝜆_d – then we will have u = v₁ which is the eigenvector corresponding to the largest eigenvalue 𝜆₁.

Thus the direction u such that projection onto it yields the least total projection error is given by the eigenvector of the scatter matrix S (equivalently the covariance matrix) with the largest eigenvalue.
This direction also captures the maximum variance of the data.

It is easy to show this if you express u = Va where a is a d x 1 vector

93 of 130

PCA

PCA thus projects the data onto that direction that minimizes the total squared difference between the data-points and their respective projections along that direction.

This equivalently yields the direction along which the spread (or variance) will be maximum.

Why? Note that the eigenvalue of a covariance matrix tells you the variance of the data when projected along that particular eigenvector:

This term is proportional to the variance of the data when projected along u.

94 of 130

PCA

But for most applications (including face recognition), just a single direction is absolutely insufficient!

We will need to project the data (from the high-dimensional, i.e. d-dimensional space) onto k (k << d) different mutually perpendicular directions.

What is the criterion for deriving these directions?

We seek those k directions for which the total reconstruction error of all the N images when projected on those directions is minimized.

95 of 130

PCA

We seek those k directions for which the total reconstruction error of all the N images when projected on those directions is minimized.

One can prove that these k directions will be the eigenvectors of the S matrix (equivalently covariance matrix of the data) corresponding to the k-largest eigenvalues. These k directions form the eigen-space.

If the eigenvalues of S are distinct, these k directions are defined uniquely (up to a sign factor)

96 of 130

PCA

One can prove that these k directions will be the eigenvectors of the S matrix (equivalently covariance matrix of the data) corresponding to the k-largest eigenvalues. These k directions form the eigen-space.
We will prove this for k = 2, other values for k will follow very similarly.

97 of 130

For any column vector z, z^tz = |z|² = trace(zz^t)

Here we are assuming that 𝜆₁ ⩾ 𝜆₂ ⩾ … ⩾ 𝜆_d.

Thus, the optimal U consists of the eigenvectors of S corresponding to the two largest eigenvalues. This derivation easily extends for any larger k as well.

98 of 130

PCA: Algorithm

Compute the mean of the given points:

Deduct the mean from each point:

Compute the covariance matrix of these mean-deducted points:

99 of 130

PCA: algorithm

4. Find the eigenvectors of C:

5. Extract the k eigenvectors corresponding to the k largest eigenvalues. This is called the extracted eigenspace:

There is an implicit assumption here that the first k indices indeed correspond to the k largest eigenvalues. If that is not true, you would need to pick the appropriate indices.

100 of 130

PCA: algorithm

6. Project each point onto the eigenspace, giving a vector of k eigen-coefficients for that point.

100

We are representing each face as a linear combination of the k eigenvectors corresponding to the k largest eigenvalues. The coefficients of the linear combination are the eigen-coefficients.

Note that α_ik is a vector of the eigencoefficients of the i-th sample point, and it has k elements. The j-th element of this vector is denoted as α_ik (j).

101 of 130

PCA and Face Recognition: Eigen-faces

Consider a database of cropped, frontal face images (which we will assume are aligned and under the same illumination). These are the gallery images.

We will reshape each such image (a 2D array of size H x W after cropping) to form a column vector of d = HW elements. Each image will be a vector x_i, as per the notation on the previous two slides.

And then carry out the six steps mentioned before.

The eigenvectors that we get in this case are called eigenfaces. Each eigenvector has d elements. If you reshape those eigenvectors to form images of size H x W, those images look like (filtered!) faces.

101

102 of 130

Example 1

102

A face database

http://people.ece.cornell.edu/land/courses/ece4760/FinalProjects/s2011/bjh78_caj65/bjh78_caj65/

103 of 130

103

Top 25 Eigen-faces for this database!

http://people.ece.cornell.edu/land/courses/ece4760/FinalProjects/s2011/bjh78_caj65/bjh78_caj65/

104 of 130

PCA and Face recognition: �Eigenfaces

For each gallery image, you compute the eigen-coefficients. You then store the eigen-coefficients and the identity of the person in a database.
You also store in the database.
During the testing phase, you are given a probe image (say) z_p in the form of a column vector of HW elements.
You deduct the mean image from z_p:

104

105 of 130

PCA and Face recognition: �Eigenfaces

You then project the mean-deducted face image onto the eigen-space:

Now, compare α_p with all the α_ik(eigen-coefficients of the gallery images) in the database.

Find the closest match in terms of the squared distance between the eigen-coefficients. That gives you the identity (see next slide).

105

Eigen-coefficients of the probe image z_p.

106 of 130

PCA and Face recognition: �Eigenfaces

106

Eigen-coefficients of the probe image z_p.

Eigen-coefficients of the l-th gallery image x_l.

Note: other distance measures (different from sum of squared differences) may also be employed. One example is sum of absolute differences, given as follows:

Another could be normalized dot product (and this distance measure should be maximized!):

107 of 130

PCA and Face recognition: eigenfaces

The eigen-face images contain more and more high frequency information as the corresponding eigen-values decrease.
Although PCA is a technique known for a long time, its application in face recognition was pioneered by Turk and Pentland in a classic paper in 1991.

M. Turk and A. Pentland (1991). Eigenfaces for recognition, Journal of Cognitive Neuroscience, 3(1): 71–86.

107

108 of 130

108

http://www.cs.ucsb.edu/~mturk/

http://web.media.mit.edu/~sandy/

109 of 130

PCA and Face recognition: eigenfaces

We can regard the k eigenfaces as key signatures.
We express each face image as a linear combination of these eigenfaces, i.e. the average face + (say) 3 times eigenface 1 + (say) 5 times eigenface 2 + (say) -1 times eigenface 3 and so on. (note: 3,5,-1 here are the eigen-coefficients, and some of them can be negative).

109

110 of 130

One word of caution: Eigen-faces

The algorithm described earlier is computationally infeasible for eigen-faces, as it requires storage of a d x d Covariance matrix (d – the number of image pixels - could be more than 10,000). And the computation of the eigen-vectors of such a matrix is a O(d³) operation!

We will study a modification to this that will bring down the computational cost drastically.

110

111 of 130

Eigen-faces: reducing computational complexity.

Consider the covariance matrix:

It will require too much memory if d is large, and computing its eigenvectors will be a horrendous task!

Consider the case when N is much less than d. This is very common in face recognition applications. The number of training images is usually much smaller than the size of the image.

111

112 of 130

Eigen-faces: reducing computational complexity.

In such a case, the rank of C is at the most N-1. So C will have at the most N-1 non-zero eigen-values.

We can write C in the following way:

112

113 of 130

Back to Eigen-faces: reducing computational complexity.

Consider the matrix X^TX (size N x N) instead of XX^T (size d x d). Its eigenvectors are of the form:

113

Xw is an eigenvector of C=XX^T! Computing all eigenvectors of C will now have a complexity of only O(N³) for computation of the eigenvectors of X^TX + O(N x dN) for computation of Xw from each w = total of O(N³ + dN²) which is much less than O(d³). Note that C has at most only min(N-1,d) eigenvectors corresponding to non-zero eigen-values (why?).

114 of 130

Eigenfaces: Algorithm (N << d case)

Compute the mean of the given points:

Deduct the mean from each point:

Compute the following matrix:

114

115 of 130

Eigen-faces: Algorithm (N << d case)

4. Find the eigenvectors of L:

5. Obtain the eigenvectors of C from those of L:

6. Unit-normalize the columns of V.

7. C will have at most only N eigenvectors corresponding to non-zero eigen-values*. Out of these you pick the top k (k < N) corresponding to the largest eigen-values.

* Actually this number is at most N-1 – this is due to the mean subtraction, else it would have been at most N.

115

116 of 130

Example 1

116

A face database

http://people.ece.cornell.edu/land/courses/ece4760/FinalProjects/s2011/bjh78_caj65/bjh78_caj65/

117 of 130

117

Top 25 Eigen-faces for this database!

http://people.ece.cornell.edu/land/courses/ece4760/FinalProjects/s2011/bjh78_caj65/bjh78_caj65/

118 of 130

Example 2

118

The Yale Face database

119 of 130

119

Top 25 eigenfaces from the previous database

Reconstruction of a face image using the top 1,8,16,32,…,104 eigenfaces (i.e. k varied from 1 to 104 in steps of 8)

120 of 130

What if both N and d are large?

This can happen, for example, if you wanted to build an eigenspace for face images of all people in Mumbai.
Divide people into coherent groups based on some visual attributes (eg: gender, age group etc) and build separate eigenspaces for each group.

120

121 of 130

How to pick k in an actual application?

Trial and error. Usually between 50 to 100 for images of size 200 x 200.
Divide your training set into two parts – A and B (B is usually called the validation set). Pick the value of k that gives the best recognition rate on B when you train on A.
Stick to that value of k.
Note: a larger k implies a better reconstruction but it may even cause a decrease in the recognition accuracy!

121

122 of 130

How to pick k in an actual application?

Note: a larger k implies a better reconstruction but it may even cause a decrease in the recognition accuracy!

Why – because throwing out some of the eigenvectors may lead to filtering of the data, removing some unnecessary artifacts for example.

122

123 of 130

Some observations about PCA for face images

The matrix V (with all columns) is orthonormal. Hence the squared error between any image and its approximation using just top k eigenvectors is given by:

123

This error is small on an average for a well-aligned group of face images – we will see why on the next slide.

Vector whose first k entries are identical to those in α_i, and the rest are 0.

124 of 130

124

The eigenvalues of the covariance matrix typically decay fast in value (if the faces were properly normalized). Note that the j-th eigenvalue is proportional to the variance of the j-th eigencoefficient, i.e.

What this means is that the data have low variance when projected along most of the eigenvectors, i.e. effectively the data are concentrated in a lower-dimensional subspace of the d-dimensional space.

X-axis = eigenvalue index (in sorted, descending order)

Y-axis = eigenvalue

125 of 130

Example

125

The Yale Face database

126 of 130

126

Top 25 eigenfaces from the previous database: each eigenface is an eigenvector of the covariance matrix reshaped to form an image

Reconstruction of a face image using the top 1,8,16,32,…,104 eigenfaces (i.e. k varied from 1 to 104 in steps of 8)

127 of 130

PCA Algorithm using SVD

Compute the mean of the given points:

Deduct the mean from each point:

Compute the covariance matrix of these mean-deducted points:

128 of 130

PCA Algorithm using SVD

4. Instead of finding the eigenvectors of C, we find the left singular vectors of X and its singular values

5. Extract the k eigenvectors in U corresponding to the k largest singular values to form the extracted eigenspace:

There is an implicit assumption here that the first k indices indeed correspond to the k largest eigenvalues. If that is not true, you would need to pick the appropriate indices.

U,S,V are obtained by computing the SVD of X.

129 of 130

LU Decomposition, QR decomposition and Cholesky decomposition

See lecture material here (SSO accessible).

QR factorization: lecture material (see board-work) taken from here (See part on Gram-Schmidt process).

Cholesky decomposition: see board work; a derivation similar to the one I gave in class is here.

129