Math 350 A Best Basis The Covariance Matrix - Study Guide

A Best Basis The problem we want to consider is this: Suppose you’re given some data that is expressed as p points in IRn . If we suppose that the data actually lies on a low dimensional subspace of IRn , is it possible to find that low dimensional basis? First, let’s look at an extreme example- Just three points, x(1) , x(2) , x(3) , each in IR4000 . If this is all our data, and we are able to store any basis we like, then a much more efficient basis would be the three vectors. That is, x(1) = 1 · x(1) + 0 · x(2) + 0 · x(3) so that x(1) would simply be stored as (1, 0, 0). Similarly, x(2) could be stored as (0, 1, 0), and the last vector as (0, 0, 1). In this example, we are trading the 4000 dimensional space for the lower, 3-dimensional space. In general, we will assume that there is a low dimensional subspace that encapsulates most of the data (rarely will data lie completely on a low dimensional subspace without some error). Our goal is to find a basis for that subspace. There are a couple of things that will make our job easier: • We will assume that the data has been mean-subtracted, so that the mean is zero. • The basis is orthonormal. • To find the “best” basis will require an error function. We will then minimize it. At the end of this section, you’ll see that the best k−dimensional basis for your data (regardless of k) is given by the first k eigenvectors of the covariance matrix, which are typically computed using the Singular Value Decomposition (SVD). The Covariance Matrix As we recall, the covariance matrix for data in IRn measures the covariance between the data in coordinate i and the data in coordinate j. If we have p vectors in IRn organized as usual in columns, x(1) , . . . , x(p) , then the k th coordinate of data point i will be denoted as (i) xk and if we assume that the mean in each coordinate is zero, p 1 X (i) x =0 p k=1 k then the covariance between coordinate i and coordinate j will be p 1 X (n) (n) x x Cij = p − 1 n=1 i j 1 (n) (n) Now, xi is the ith “row” of x(n) and xj matrix can be expressed as: is the j th “column” of x(n) , so the entire covariance n 1 X (n) (n)
T x x C= p − 1 p=1 If the matrix X stores the data in columns so that X is n × p, we can write: C= 1 XX T p−1 If the matrix X stores the data in rows so that X is p × n, we can express this as C= 1 XT X p−1 In either event, the covariance matrix is n × n. Finally, we recognize that the covariance matrix is symmetric, so the Spectral Theorem applies. In particular, there is an orthonormal matrix P and a diagonal matrix D so that C = P DP T where the columns of P form the eigenvectors associated with the diagonal elements of D (which are typically written largest to smallest). To connect this with the SVD of the data matrix X, if X is n × p (so that data is stored column-wise), and we write the reduced SVD as: X = U ΣV T Then   1 1 1 T T T T 2 XX = U ΣV V Σ U = U Σ UT C= p−1 p−1 p−1 We see a relationship between the singular values of X, σi , and the eigenvalues of the covariance matrix, λ̂i : 1 σi2 = λ̂i p−1 And, of course the most important conclusion: The eigenvectors of the covariance matrix are the left singular vectors of the data matrix (when the data is written column-wise and has been mean subtracted). Suppose that the matrix X is 3, 000, 000 × 200. Then the eigenvectors of the covariance matrix are in IR3,000,000 , and the matrix P (or matrix U ) would be a huge matrix to have to compute. There is a better way. Suppose that we know the right singular vectors in the matrix V (which reside in IR200 ). Is there a fast way to compute a few of the columns in U from that? Yes! Recall that this is how we originally defined the vectors in U : ui = 1 Xvi σi Or, more practically, compute Xvi and divide it by its norm. 2 Projections and the Mean We now look at something a little different, but we’ll tie it together later. Suppose you have your p data points in IRn that have not been mean subtracted, and you have a vector u onto which we want to project the data. First, if we project one point x onto our (unit) vector u, then the point becomes the scalar xT u. Similarly, projecting all the data gives us p real numbers: uT x(1) , uT x(2) , . . . , uT x(p) so the mean of the projected data is given by p p 1 X (n) x p n=1 1 X T (n) u x = uT p n=1 ! = uT x̄. Therefore, the mean of the projection is the projection of the mean. In particular, if the mean of a data set is zero, and the data is projected to a vector (or subspace), then the new mean is also zero. Projections and the Basis In this section, let’s think about projecting the data to various one dimensional subspaces and how that effects the variance of the data. We’ll keep our previous general data set, p points in IRn , and we’ll suppose that the data has been mean subtracted (although we could do the mean subtraction after the projection, as we showed in the previous section). If we project the data onto an arbitrary unit vector u, then the resulting variance (as scalars) will be: p S~u2 1 X T (n) 2 (u x ) = p − 1 n=1 p 1 X T (n) (n) T u x (x ) u = p − 1 n=1 p = uT 1 X (n) (n) T x (x ) p − 1 n=1 ! u = uT Cu This is actually a very key quantity, and will come up in the next section. We will look at this quantity more closely in a bit, but let’s look at what happens in one special case: Suppose that u is an eigenvector of the covariance C corresponding to the first eigenvalue of C, λ̂i using our previous notation. Then: v1T Cv1 = v1T λ̂1 v1 = λˆ1 Therefore, if we project all the data to the first eigenvector of C, the new variance will be the first eigenvalue of C. 3 Reconstruction Error and the Basis ~1, φ ~2, . . . , φ ~ n , we can Given a specific vector x ∈ IRn and an arbitrary orthonormal basis, φ write ~n ~ T x)φ ~ 1 + (φ ~ T x)φ ~ 2 + . . . + (φ ~ T x)φ x = (φ n 1 2 so that the magnitude of x can be written as: ~ T x)2 . ~ T x)2 + . . . + (φ ~ T x)2 + (φ kxk2 = (φ n 2 1 We can use the same algebraic manipulation that we used in the last section to rewrite this as: ~n. ~2 + . . . + φ ~ T xxT φ ~1 + φ ~ T xxT φ ~ T xxT φ kxk2 = φ 2 1 n ~1: We can break this up and define the error using one vector φ ~ T xxT φ ~ 1 + kxerr k2 kxk2 = φ 1 ~ 1 , we sum these together: Now do this for all p data points. For any single vector φ ~ 1 + kx(1) k2 ~ T x(1) (x(1) )T φ kx(1) k2 = φ err 1 (2) 2 T (2) (2) T ~ 1 + kx(2) k2 ~ x (x ) φ kx k = φ err 1 .... .. ~ 1 + kx(p) k2 ~ T x(p) (x(p) )T φ kx(p) k2 = φ 1 err p X ~T k =φ 1 (n) 2 kx n=1 p X ! (n) x (n) T (x ) ~1 + φ n=1 p X 2 kx(n) err k n=1 We can multiply everything by 1/p or 1/(p − 1) to make things work. That is, p p 1 X (n) 2 ~ T ~ 1 X (n) 2 kx k = φ1 C φ1 + kx k . p − 1 n=1 p − 1 n=1 err (1) In light of this equation, let us now define an error function using an arbitrary orthonormal ~1, . . . , φ ~ n . The error we get when using a one dimensional representation of our data basis, φ is given by p 1 X (n) 2 ~ ~ E(φ2 , . . . , φn ) = kx k . p − 1 n=1 err Notice that the left side of Equation 1 is constant. Therefore, minimizing our error function ~T C φ ~1. is equivalent to maximizing φ 1 Here now is our algorithm for find the “best” basis: 4 ~ 1 so that φ ~T C φ ~ 1 is maximized. 1. Find the unit vector φ 1 2. We “project out” this vector so that the ith data point now becomes: x(i) = x(i) − Projφ~1 (x(i) ) 3. Re-compute C. 4. Repeat from Step 1 until we have enough basis vectors. In practice, we will not need to do this- there is an easier way! The Best Basis and the Eigenvectors We have shown that finding the best basis reduces to maximizing the quantity: ~T C φ ~ φ max ~ ~T φ ~ =~0 φ6 φ where we divide by the magnitude (squared) to enforce the fact that we want a unit vector, and we want to stay away from the zero vector. We know that the eigenvectors of the covariance matrix form an orthonormal basis for n IR , so we can write any vector as a linear combination of them: ~ = c1 v1 + c2 v2 + · · · + cn vn = V c φ Using the eigenvector-eigenvalue factorization C = V DV T , we can now write the numerator as: ~T C φ ~ = (V c)T (V DV T )(V c) = cT (V T V )D(V T V )c = cT Dc = c2 λ1 + c2 λ2 + · · · + c2 λn φ 1 2 n Similarly, the denominator is: ~T φ ~ = c2 + c2 + · · · + c2 φ 1 2 n Let’s look at the coefficients of our expansion now. For λi , the coefficient in front is ρi = c2i c21 + c22 + · · · + c2n P where ρi ≥ 0 and ni=1 ρi = 1 (like a probability distribution). Let’s summarize where we ~T C φ ~ is equivalent to choosing ρ1 , ρ2 , · · · , ρn so that each are. We now see that maximizing φ ρi ≥ 0 and they sum to 1, to maximize the quantity: ρ1 λ1 + ρ2 λ2 + · · · + ρn λn 5 It is easy to see that, if the λi ≥ 0 and are ordered from largest to smallest, then: λn ≤ ρ1 λ1 + ρ2 λ2 + · · · + ρn λn ≤ λ1 . To maximize our given quantity, we set c1 = 0 and the rest of the coefficients to zero. This leads us to our main conclusion: ~ that maximizes the quantity The Best Basis The vector φ max ~ =~0 φ6 ~T C φ ~ φ ~T φ ~ φ is given by v1 , the eigenvector corresponding to the largest eigenvalue of C. Therefore, the best k−dimensional basis for a given set of data is found by taking the first k eigenvectors of the covariance matrix. We note in passing that the value of the maximum would then be the first eigenvalue of the covariance matrix. Alternatively, given the data in an n × p matrix X, the best k−dimensional basis is found by taking the first k columns of the U , the left singular vectors of the SVD of X. The Best Basis and The Rank Before we start the decomposition, it would be a good idea to decide on how the rank of X will be computed- Matlab uses a very very small tolerance on its rank cutoff, so it will often return the maximum rank. It is useful to look at the rank as that number of basis vectors required to preserve some percentage of the variance in the data. From our previous section on the covariance matrix, we had a relationship between the eigenvalues of the covariance matrix, λ̂i and the singular values of X: 1 σi2 = λ̂i p−1 so normalizing the set of eigenvalues is equivalent to doing it to the squared singular values: λ̂i λi = Pn j=1 λ̂j = 1 σ2 p−1 i Pn 1 2 j=1 p−1 σj σ2 = Pn i j=1 σj2 Now the λi are positive and sum to 1. How many dimensions we keep is usually problem dependent. There are some problems where keeping 60% of the variance is more than enough, and some problems will come up where we need 99% of the variance. Best Basis and Matlab In this section, we look at the Matlab commands that we use in finding and analyzing the best basis. 6 Suppose we have p points in IRn , and these are stored in an n × p matrix X. First, we mean-subtract the data (as vectors in IRn ). I will call the mean subtracted data X̂. Once we mean subtract, find the SVD. [n,p]=size(X); Xmean=mean(X,2); Xhat=X-repmat(Xmean,1,p); [U,S,V]=svd(Xhat,’econ’); temp=diag(S).^2; Lambda=temp/sum(temp); %Lambda contains the normalized evals Now we can plot the data in Lambda to see what the eigenvalues look like: plot(Lambda); If there is a break in the values where suddenly all the remaining eigenvalues are zero, then it is relatively easy to determine the rank. Suppose we want to find and plot the best two dimensional representation to the data in X. This low dimensional representation is provided by the coordinates, and we find them using matrix notation: Xcoords=U(:,1:2)’*Xhat; %Size of Xcoords is now 2 x p (p points in R^2) plot(Xcoords(1,:),Xcoords(2,:),’.’); To put the data back into IRn (as a 2 or k dimensional approximation), we would just multiply by the appropriate columns of U : Xrecon=U(:,1:2)*Xcoords; %Project the data back to R^n An Example in Detail Download the data author.mat from the class website. The data represents a short movie1 consisting of 109 frames, where each frame is 120 × 160, or 19200 pixels. Thus, loading the data will give you an array that is 19200 × 109 (this is what we were referring to as n × p in the notes). Also, we should note that the data in matrix, saved as Y1, consists of integer values between 0 and 255 (representing 255 shades of gray from black (0) to white (255)). This is important to mention: we’ll be translating these images into vectors of real numbers instead of integers, which may impact how the coloring takes place. Discussing this in detail now would take us too far away from the discussion, but is something to keep in mind. 1 The movie is actually a short sequence I took from my webcam when I was on sabbatical and needed some data quick! 7 Visualizing a vector as an array When visualizing these large scale vectors, remember that any vector in IR19200 could be visualized as a 120×160 array. In Matlab, this is performed by using the “reshape” command. We use imagesc, short for ”image scaled”, where Matlab will scale the values in the array to ready it for the color map. To see how the color mapping is performed, you can use the colorbar command. In this case, if ~u is a vector in IR19200 , then we can reshape it into an array and visualize it as an array using a gray scale color scheme by issuing the following commands, issued on one line imagesc(reshape(u,120,160)); colormap(gray); axis equal; colorbar To watch the “movie”, you can use the following code, which just goes through the data column by column and shows each as frame. I won’t use the color bar, and the pause command is there so we can actually see the result (otherwise, the frames go by way too fast!). for j=1:109 A=reshape(Y1(:,j),120,160); imagesc(A); colormap(gray); axis equal; axis off; pause(0.1) end A few sample frames and the movie’s mean frame are shown in Figure 1. You might find it interesting that that the mean is the background. Once we’ve found the mean, then we can proceed with the Singular Value Decomposition. Now consider the following code “snippet”: X=Y1-repmat(Ymean,1,109); [U,S,V]=svd(X,’econ’); tempS=diag(S).^2; tempS=tempS/sum(tempS) The last line stores the normalized eigenvalues (singular values squared) into the vector tempS. Now let’s consider the rank. Here are the first few of the values in that vector: 0.3436 0.1229 0.0847 0.0566 0.0407 There seems to be a natural break between the first two and the third. Keeping the data two dimensional would retain approximately 34.36 + 12.29 = 46.6% of the variance in the data. Keeping three dimensions takes us to about 55%. We need to go to 22 dimensions to retain 90%. In the script file, we’ll compare a two dimensional reconstruction, a 10 dimensional reconstruction and a 22 dimensional reconstruction of the data. This is done by using two, 10, then 22 basis vectors for the subspace in which the data lies. First, let’s examine the actual basis vectors. From our code snippet, the basis vectors of interest are the columns of U (the left singular vectors of X). Since each basis vector is in IR19200 , the basis vector(s) can also be visualized as an array. 8 Figure 1: The figure above has three sample frames (Frames 23, 60 and 80) and the movie’s mean. Sample Code %% Sample script for the "author" data % Some cleaning commands: close all; %Close all open figures clear %Clear the memory clc %Clear the command window % Download the data from the class website, so it can be loaded here: load author Y1=double(Y1); %Changes the data from "integers" to "floating point" %The data was saved as integers to save on memory. % We can visualize the "movie" figure(1) for j=1:109 A=reshape(Y1(:,j),120,160); imagesc(A); colormap(gray); axis equal; axis off; pause(0.1) end %% Mean subtraction Ymean=mean(Y1,2); 9 figure(2) imagesc(reshape(Ymean,120,160)); colormap(gray); axis equal; axis off; title(’The Movie Mean’); %% The SVD % We’ll mean subtract, but keep the original data intact. % the mean subtracted data in matrix X. X=Y1-repmat(Ymean,1,109); [U,S,V]=svd(X,’econ’); tempS=diag(S).^2; tempS=tempS/sum(tempS); %Normalized eigenvalues. We’ll call put figure(3) plot(tempS); title(’The Eigenvalues of X^TX (Singular Values squared)’) %% Projection to R^2 Xcoords=U(:,1:2)’*X; %The matrix Xcoords is now 2 x 109 figure(4) plot(Xcoords(1,:),Xcoords(2,:),’.’); title(’The movie data presented in the best two dimensional basis’); %% The Reconstruction back into 19200 dim space Xrecon=U(:,1:2)*Xcoords; %Add back the mean if needed. Exercises 1. Suppose we have p data points in IRn . Show that the variance of the data, projected to a standard basis vector ~ei , returns the standard variance in that dimension. 2. By hand, be able to reproduce the computation that shows that the variance of the scalar projection of the p data points in IRn to a unit vector u is given by uT Cu 3. Suppose we have two o.n. vectors u, v ∈ IRn . If we still have our p points in IRn , compute a quantity similar to the last exercise that gives the covariance between the data projected to u and the data projected to v. In the special case that u, v are eigenvectors of the covariance matrix, how does this quantity simplify? 10 4. Relate the variance to the error. That is, be able to reproduce the argument that ! p p p X X X ~1 + ~T kx(n) k2 x(n) (x(n) )T φ kx(n) k2 = φ err 1 n=1 n=1 n=1 5. In Matlab, run the sample code for author.mat, and answer the following: (a) Verify numerically that the variance of the projected data to the best basis vector (first one) is given by the first eigenvalue of the covariance matrix. (b) Verify numerically that we have a good basis in using the first two eigenvectors of the covariance matrix (the first two columns of U in the SVD) by creating 500 random orthonormal pairs of vectors and comparing the reconstruction errors. Here’s a start for your code: %Assume X is given as in the previous "author" code. for j=1:500 [Q,temp]=qr(randn(19200,2),0); %Q that has o.n. cols, discard temp Error=???? (A matrix whose columns are the x_err) ErrJ=norm(Error,’fro’); end Project: Eigenfaces Suppose we have a large collection of photos of people’s faces, and that each photo is u × v. Then we think of each “face” as a vector in IRn (where n = uv). If we look for a “best basis” in this space, we will get basis vectors that are also in IRn , and these can be viewed as an array of u × v pixels. Therefore, if the basis vectors are the eigenvectors of the covariance matrix, we might call such vectors “eigenfaces”. We already know how to compute such a basis. There is data on our class website where each face has 262 × 294 = 77, 028 pixels. Our project is to: 1. Compute the mean vector and represent it as a face. 2. Plot a sample of 5 faces and the mean face together in one plot. 3. Compute the first six eigenfaces (and represent them as faces). Put them together in one plot. 4. What might be a good approximation to the rank of the matrix (if our goal is to make the faces recognizable, about 74% of the variance). 5. Plot the data using the best two dimensional representation using a new figure. Plot the figures representing females as asterisks and the figures representing males as triangles. Do you see anything interesting? 11