SVD and PCA - Math for Machine Learning

SVD and PCA Math for Machine Learning 1 Adjoints How are a matrix A and its transpose A> related? Note that hAx, yi = (Ax)> y = x> A> y = hx, A> yi. So, when taking inner products, it seems that A> represents the action of A on the opposite argument. In fact, this relationship connects back to the underlying linear map T represented by A. We define the adjoint of T , T ∗ , to be the linear map represented by A> . Equivalently, B(T (v), w) = C(v, T ∗ (w)) for any appropriate nondegenerate bilinear forms B, C. (We will not define this term, as it is very far outside the scope of the class. Suffice to say that the inner products h·, ·irange(T ) , h·, ·irange(T ∗ ) can be filled in for B and C, although inner products are not the only examples of nondegenerate bilinear forms.)1 We now explore how the adjoint T ∗ connects back to A. Observe that any vector of the form Av is a linear combination of the columns of A. Likewise, any vector of the form A> w is a linear combination of the rows of A. Therefore, Im(T ) = range(A) = col(A), the column space of A; Im(T ∗ ) = range(A> ) = row(A), the row space of A. For any linear map T : Rn → Rm (or the associated matrix A ∈ Rm×n ), we therefore have the following four ”fundamental subspaces” associated with T : 1. Im(T ), a.k.a. col(A) 2. ker(T ), a.k.a. null(A) 3. Im(T ∗ ), a.k.a. row(A) 4. ker(T ∗ ), a.k.a. null(A> ). 2 Fundamental Theorem of Linear Algebra and SVD It turns out these subspaces are related in a way captured by the first part of the Fundamental Theorem of Linear Algebra: Theorem (FTLA, Part I) Let A ∈ Rm×n . Then (a) Rm = range(A) ⊕ ker(A> ). (b) rank(A) + dim ker(A) = n. 1 If you want to learn more about adjoints, make sure you can solve all of Q4 on Homework 2, as understanding dual spaces is essential to understanding adjoints. 1 Proof Proving (a) amounts to showing that ker(A> ) = range(A)⊥ . We have x ∈ ker(A> ) ⇐⇒ A> x = 0 ⇐⇒ a> i x = 0 for all i = 1, . . . , n (where ai is the i’th column of A) ⇐⇒ v> x = 0 for all v ∈ range(A) ⇐⇒ x ∈ Im(A)⊥ . Part (b) hinges on the fact that rank(A) = rank(A> ). Assuming that fact, apply (a) combined with that dim(U ⊕ V ) = dim U + dim V .  Part (b) is usually known as the Rank-Nullity Theorem. Singular Value Decomposition Given A ∈ Rm×n , the FTLA gives us two natural-looking orthogonal decompositions involving the four fundamental subspaces of A: (i) Rn = ker(A) ⊕ range(A> ) (ii) Rm = ker(A> ) ⊕ range(A) To dig deeper, we must examine the matrices A> A ∈ Rn×n and AA> ∈ Rm×m . From an intuitive perspec> tive, note that (A> A)i,j = a> i aj , where ai is the i’th column of A. Thus, A A gives some measure of the > similarity between the columns of A. Similarly, AA measures similarity between the rows of A. Lemma 1. If A has full column rank (i.e. A has n linearly independent columns), then A> A is invertible. Proof. Let x ∈ ker(A> A), i.e. A> Ax = 0. Then Ax ∈ ker(A> ). By FTLA, Ax ∈ range(A)⊥ . But clearly Ax ∈ range(A) as well, so Ax = 0 by the fact that range(A) ⊥ range(A)⊥ . Finally, since A is full-rank, Ax = 0 implies that x = 0. Thus, A> A is a square matrix with trivial kernel, so it is invertible.  Lemma 2. ker(A> A) = ker(A), so rank(A> A) = rank(A). Proof. Exercise. Lemma 3. A> A is positive semi-definite (PSD). Proof. Clearly (A> A)> = A> A, so it is symmetric. To show that it is PSD, see that x> A> Ax = hAx, Axi ≥ 0 with equality iff Ax = 0.  To sum up what we know about A> A: (i) It preserves the kernel and rank of A (ii) It is PSD. If we take A to be invertible (thus encoding a change of basis in Rn ), we can consider the transformation x 7→ Ax. A natural question is: what happens to our standard inner product under the transformation? If the standard inner product of x, y is hx, yi = x> y = x> In y, then the standard inner product of the transformed vectors Ax, Ay is hAx, Ayi = (Ax)> Ay = x> (A> A)y. So A> A can be thought of as a “scaling factor” by which we can recover the standard inner product in the transformed space under x 7→ Ax. 2 Existence of SVD Let’s examine the spectrum2 of A> A. Since A> A is PSD, its eigenvalues are all ≥ 0. Since A> A is symmetric, it has a spectral decomposition A> A = VΛV> More precisely, if A> A has r nonzero eigenvalues (with multiplicity), then write  0   0
Λ 0 V = A> A V0 V0 , 0 0 V0 where Λ0 is a diagonal matrix with the nonzero eigenvalues of A> A corresponding to eigenvectors in V0 and the eigenvectors with vanishing eigenvalue are in V0 . Next, define 1 U0 = AV0 Λ0− 2 ∈ Rm×r . Then we have 1 1 1 U0 Λ0 2 V0> = AV0 Λ0− 2 Λ0 2 V0> = AV0 V0> = A because V is unitary. Moreover, 1 1 U0> U0 = (AV0 Λ0− 2 )> AV0 Λ0− 2 1 1 = Λ0− 2 V0> A> AV0 Λ0− 2 1 1 = Λ0− 2 Λ0 Λ0− 2 = Ir , so the columns of U0 are orthonormal and can be extended to form an orthonormal basis for Rm . If we choose U0 containing these added columns, then
U = U0 U0 is unitary. Next, we form  Σ= 1 Λ2 0  so that Σ has m − r rows of zeros at the bottom and is thus in Rm×n . We arrive at   0
Λ 12 V UΣV> = U0 U0 V0 0  
Λ 12 V0> = U0 U0 0 1 = U0 Λ 2 V0> = A. The decomposition A = UΣV> is known as a singular value decomposition of A, and we have just proven its existence. To sum up, 2 A.k.a., the ”eigenstuff”. 3 Theorem (SVD) Given any matrix A ∈ Rm×n , there exist orthonormal bases {v1 , . . . , vn } of Rn and {u1 , . . . , um } of Rm such that if V = (v1 , . . . , vn ), U = (u1 , . . . , um ), and Σ ∈ Rm×n is a rectangular diagonal matrix containing the square roots of the nonzero eigenvalues of A> A, then A = UΣV> . The square roots of the nonzero eigenvalues of A> A are called the singular values of A.3 They are usually denoted (σi ). The corresponding columns of V, i.e. the eigenvectors (vi ) of A> A corresponding to the (σi ), are called right-singular vectors of A. The corresponding columns of U are called left-singular vectors of A. The singular values of A are unique, but the corresponding singular vectors are not. Note that the columns of U form an eigenbasis of Rm with respect to AA> : AA> = UΣV> (UΣV> )> = UΣΣ> U > = UΛU> for appropriately sized Λ. Similarly, we defined V via the spectral decomposition of A> A: A> A = VΛV> for appropriately sized Λ. From our discussion about A> A, it follows that the number of singular values of A is equal to rank(A). The geometric picture for SVD goes as follows: 1. First, via the unitary V> , change coordinates to the eigenbasis for A> A. 2. Via Σ, which has the same rank as A, scale by the σi . 3. Via the unitary U, rotate back. Note that in Step 1, the actions of the nonsingular vi don’t matter. Why? Because those vi correspond to eigenvalue 0, so they are in ker(A> A), which we established is the same as ker(A). More thoroughly, Theorem (FTLA, Part II) Let A = UΣV> be a singular value decomposition of A, and let rank(A) = r. Then (i) The first r columns of V, i.e. the right-singular vectors of A, form an orthonormal basis for range(A> ). (ii) The last n − r columns of V form an orthonormal basis for ker(A). (iii) The first r columns of U, i.e. the left-singular vectors of A, form an orthonormal basis for range(A). (iv) The last m − r columns of U form an orthonormal basis for ker(A> ). Proof Easy to fill in the details from the comment directly preceding the theorem combined with FTLA, Part I.  When passing to the eigenbasis of A> A via V> , we effectively ignore vectors in ker(A), as Σ will kill them. As for the relevant coordinates, we let the singular vectors transform them, Σ scale them, and then U bring them back into the image where they belong. 3 It’s worth noting that A> A and AA> have the same nonzero eigenvalues. 4 3 Application: Principal Component Analysis (PCA) Motivation When fitting a model (e.g. an OLS model) to data, we hope to represent the data in a simpler way. A basketball player is a complex object; a vector representing (points, assists, rebounds, eye color, birthday) is not. However, we hope to use only the most relevant features for (i) fast computation and (ii) to build a more stable model (i.e. reduce the variance of the model). How can we find the features that are important for predicting whether a player helps his team win (probably points, assists, rebounds), allowing us to ignore the features that are not (probably eye color, birthday)? Given a matrix of data X ∈ Rn×d containing n-many d-dimensional data points, PCA allows us to find a suitable subspace of Rd onto which we can project our data, leaving us with the most relevant features. How do we decide which features to drop? The idea is that we look at the data and keep only a small number (< d) of orthogonal directions (perhaps linear combinations of features) that capture the most variance of the data. Intuitively, the low-variance directions contain less information about the data, so we can throw them away, improving the model’s performance on new data without hurting predictive accuracy. Understanding Variance The first step of PCA is to center the data so that every feature has mean 0 amongst the data points. We do this because uncentered data would influence our choice of relevant directions (which will be unitlength arrows Pn starting at 0) in unwanted ways. Thus, we first center X by subtracting the mean vector Ex = n1 i=1 xi from each row. Next, we hope to find a unit vector v ∈ Rd that will capture the most “variance in the data.” What does this mean? We will ultimately project every data point xi onto v by taking x> i v, so we care about the variance of this quantity as i ranges through the data from 1 to n. A brief aside on variance in probabilistic terms: Given a random variable X with EX = 0, we have the following version of Chebyshev’s inequality: P(|X| ≥ a) ≤ Var(X) . a2 This inequality tells us that Var(X) is a measure of how likely X is to vary from its mean. Variance conveys information about the tail probabilities of X; it shows us how likely X is to take a highly unexpected value. In our case, highly unexpected values correspond to new information. That is why we want v that maximizes d the variance of the x> i v: it would mean that span(v) is the 1-dimensional subspace of R containing the most information about the data. Finding the First Principal Component We compute the sample variance of this projection amongst our n-many data points: n sample variance = 1 1 1X > 2 (x v) = ||Xv||2 = v> X> Xv. n i=1 i n n Thus, if we want to find the unit vector v maximizing the variance, we’ve walked into a constrained optimization problem: max v> X> Xv subject to ||v|| = 1 v∈Rd To solve this optimization problem, recall that for symmetric A ∈ Rd×d , for any v ∈ Rd with ||v|| = 1, λmin (A) ≤ v> Av ≤ λmax (A) 5 where for both bounds, equality holds iff v is a corresponding eigenvector. This immediately yields that the first loading vector v = v1 is a unit eigenvector corresponding to the maximal eigenvalue of X> X. Finding More Principal Components We often want more than one principal component. Given k −1 principal components, the problem of finding the k’th amounts to another constrained optimization problem: max v> X> Xv subject to ||v|| = 1 v∈Rd v> vi = 0 for i = 1, . . . , k − 1. The k’th loading vector vk is given by the following result: Theorem The solution to the above optimization problem is v = a unit eigenvector corresponding to the k’th largest eigenvalue of X> X. Proof By induction on k. See source at end of note.  All this tells us that we can compute the first k loading vectors by computing the SVD of X and taking the first k right-singular vectors. Projecting onto the PCA Coordinate System How do we project the data onto the subspace of Rd spanned by v1 , . . . , vk ? For each data point xi , we want to map k X (x> xi 7→ Projspan{v1 ,...,vk } (xi ) = i vj )vj . j=1 Computationally, it’s easier to handle all the xi at once through matrix multiplication. Let Vk = (v1 , . . . , vk ). e k ∈ Rn×k is given by Then the new data matrix X e k = XVk V> . X k e k are exactly what we wanted: the original data points projected onto the subspace spanned The rows of X by v1 , . . . , vk . Geometric View of PCA OLS could be viewed intuitively as finding the “line of best fit.” In it, we minimize the vertical distance between the data points and the fitted line. Similarly, we can view PCA as finding the “subspace of best fit” insofar as the k-dimensional subspace we project onto minimizes perpendicular distance between it and the original data points in Rd . To show this, we need to show that our first loading vector minimizes the reconstruction error n X ||xi − Pv (xi )||2 i=1 where Pv (xi ) = (x> i v)v represents the projection of xi onto the span of v. By the Pythagorean theorem, we have that ||xi − Pv (xi )||2 + ||Pv (xi )||2 = ||xi ||2 6 so that n X ||xi − Pv (xi )||2 = i=1 n X ||xi ||2 − ||Pv (xi )||2 i=1 = n X i=1 ||xi ||2 − n X 2 (x> i v) . i=1 The first term is constant in v, so minimizing reconstruction error amounts to minimizing which is precisely our objective in PCA. Pn > 2 i=1 (xi v) , Low-Rank Approximation Let ||·|| be any unitary-invariant norm on Rn×d . A family of such norms is the collection of induced `p -norms for matrices: ||Ax||p ||A||p = sup x6=0 ||x||p of which the operator (p = 1) and the spectral (p = 2) norms are examples. Take a matrix X ∈ Rn×d , where Rn×d is equipped with our unitary-invariant norm || · ||. If we seek the best rank-k approximation to X with respect to || · ||, then PCA comes in handy: Theorem (Eckart-Young-Mirsky) e k is the best rank-k approximation to X with respect to || · || in the sense that for any Our PCA solution X rank-r (r ≤ k) matrix Y ∈ Rn×d , e k|| ≤ ||X − Y||. ||X − X This theorem tells us that the process of projecting our data onto a subspace via PCA amounts to finding the best rank-k approximation of X. 7