Linear Algebra Fundamentals

Linear Algebra Fundamentals It can be argued that all of linear algebra can be understood using the four fundamental subspaces associated with a matrix. Because they form the foundation on which we later work, we want an explicit method for analyzing these subspaces- That method will be the Singular Value Decomposition (SVD). It is unfortunate that most first courses in linear algebra do not cover this material, so we do it here. Again, we cannot stress the importance of this decomposition enough- We will apply this technique throughout the rest of this text. 0.1 Representation, Basis and Dimension Let us quickly review some notation and basic ideas from linear algebra: Suppose that the matrix V is composed of the columns v 1 , . . . , v k , and that these columns form a basis basis for some subspace, H, in IRn (notice that this implies k ≤ n). Then every data point in H can be written as a linear combination of the basis vectors. In particular, if x ∈ H, then we can write: . x = c1 v 1 + . . . + ck v k = V c so that every data point in our subset of IRn    x1  x2     x =  .  ←→   ..  xn is identified with a point in IRk :  c1 ..  = c .  ck The vector c, which contains the coordinates of x, is the low dimensional representation of the point x. That is, the data point x resides in IRn , but c is in IRk , where k ≤ n. Furthermore, we would say that the subspace H (a subspace of IRn ) is isomorphic to IRk . We’ll recall the definition: Definition 0.1.1. Any one-to-one (and onto) linear map is called an isomorphism. In particular, any change of coordinates is an isomorphism. Spaces that 1 2 are isomorphic have essentially the same algebraic structure- adding vectors in one space is corresponds to adding vectors in the second space, and scalar multiplication in one space is the same as scalar multiplication in the second. Definition 0.1.2. Let H be a subspace of vector space X. Then H has dimension k if a basis for H requires k vectors. Given a linearly independent spanning set (the columns of V ) to compute the coordinates of a data point with respect to that basis requires a matrix inversion (or more generally, Gaussian elimination) to solve the equation: x=Vc In the case where we have n basis vectors of IRn , then V is an invertible matrix, and we write: c = V −1 x If we have fewer than n basis vectors, V will not be square, and thus not invertible in the usual sense. However, if x is contained in the span of the basis, then we will be able to solve for the coordinates of x. Example 0.1.1. Let the subspace H be formed by the span of the vectors v 1 , v 2 given below. Given the point x1 , x2 below, find which one belongs to H, and if it does, give its coordinates.         4 7 2 1 v 1 =  2  v 2 =  −1  x1 =  4  x2 =  3  −1 0 1 −1 SOLUTION: Rather than row-reduce twice, we’ll do matrix below.    1 0 1 2 7 4  2 −1 4 3 → 0 1 −1 1 0 −1 0 0 it once on the augmented 3 2 1  2 1  0 How should this be interpreted? The second vector, x2 is in H, as it can be expressed as 2v 1 +v 2 . Its low dimensional representation (its coordinate vector) is thus [2, 1]T . The first vector, x1 , cannot be expressed as a linear combination of v 1 and v 2 , so it does not belong to H. If the basis is orthonormal, we do not need to perform any row reduction. Let us recall a few more definitions: Definition 0.1.3. A real n × n matrix Q is said to be orthogonal if QT Q = I 0.1. REPRESENTATION, BASIS AND DIMENSION 3 This is the property that makes an orthonormal basis nice to work with- it’s inverse is its transpose. Thus, it is easy to compute the coordinates of a vector x with respect to this basis. That is, suppose that x = c1 u1 + . . . + ck uk Then the coordinate cj is just a dot product: x · u j = 0 + . . . + 0 + cj u j · u j + 0 + . . . 0 ⇒ cj = x · uj We can also interpret each individual coordinate as the projection of x onto the appropriate basis vector. Recall that the orthogonal projection of x onto a vector u is the following: u·x u Proju (x) = u·u If u is unit length, the denominator is 1 and we have: Proju (x) = (uT x)u = (uuT )x = u(uT x) Writing the coefficients in matrix form, with the columns of U being the orthonormal vectors forming the basis, we have: c = [x]U = U T x Additionally, the projection of x onto the subspace spanned by the (orthonormal) columns of a matrix U is: ProjU (x) = U c = U U T x (1) Example 0.1.2. We’ll change our previous example slightly so that u1 and u2 are orthonormal. Find the coordinates of x1 with respect to this basis.       1 2 −1 1 1 u1 = √  2  u2 = √  −1  x1 =  0  5 0 6 1 2 SOLUTION: c = UT x  ⇒ c= √ 1/√5 2/ 6 √ 2/√5 −1/ 6 √0 1/ 6   √   −1 3 √5  8 = −2 6 −2  The reader should verify that this is accurate. We summarize our discussion with the following theorem: Representation Theorem. Suppose H is a subspace of IRn with orthonormal basis vectors given by the k columns of a matrix U (so U is n × k). Then, given x ∈ H, 4 • The low dimensional representation of x with respect to U is the vector of coordinates, c ∈ IRk : c = UT x • The reconstruction of x as a vector in IRn is: x̂ = U U T x where, if the subspace formed by U contains x, then x = x̂- Notice in this case, the projection of x into the columnspace of U is the same as x. This last point may seem trivial since we started by saying that x ∈ U , however, soon we’ll be loosening that requirement. Example 0.1.3. Let x = [3, 2, 3]T and let the basis vectors be u1 = √12 [1, 0, 1]T and let u2 = [0, 1, 0]T . Compute the low dimensional representation of x, and its reconstruction (to verify that x is in the right subspace). SOLUTION: The low dimensional representation is given by:    √  √  √  3 1/ 2 0 1/ 2  3 2 T  2 = c=U x= 0 1 0 2 3 And the reconstruction (verify the arithmetic) is:      1/2 0 1/2 3 3 0 1 0  2  =  2  x̂ = U U T x =  1/2 0 1/2 3 3 For future reference, you might notice that U U T is not the identity, but U U is the 2 × 2 identity: √    √   √  1/ 2 0 1/ 2 0 1/ 2  = 1 0 0 1 UT U = √ 0 1 0 1 0 1/ 2 0 T Projections are important part of our work in modeling data- so much so that we’ll spend a bit of time formalizing the ideas in the next section. 0.2 Special Mappings: The Projectors In the previous section, we looked at projecting one vector onto a subspace by using Equation 1. In this section, we think about the projection as a function whose domain and range will be subspaces of IRn . The defining equation for such a function comes from the idea that if one projects a vector, then projecting it again will leave it unchanged. 0.2. SPECIAL MAPPINGS: THE PROJECTORS 5 Figure 1: Projections P1 and P2 in the first and second graphs (respectively). Asterisks denote the original data point, and circles represent their destination, the projection of the asterisk onto the vector [1, 1]T . The line segment follows the direction P x − x. Note that P1 does not project in an orthogonal fashion, while the second matrix P2 does. Definition 0.2.1. A Projector is a square matrix P so that: P2 = P In particular, Px is the projection of x. Example 0.2.1. The following are two projectors. Their matrix representations are given by:     1 1 1 1 0 P1 = P2 = 1 0 2 1 1 Some samples of the projections are given in Figure 1, where we see that both project to the subspace spanned by [1, 1]T . Let’s consider the action of these matrices on an arbitrary point:           1 0 x x 1 0 x x P1 x = = , P1 (P1 x) = = 1 0 y x 1 0 x x     x+y    1 1 1 x+y 1 x 2 P2 x = = x+y = y 1 2 1 1 2 2 You should verify that P22 x = P2 (P2 (x)) = x. You can deduce along which direction a point is projected by drawing a straight line from the point x to the point Px. In general, this direction will depend on the point. We denote this direction by the vector Px − x. From the previous examples, we see that Px − x is given by:    −x+y    x − y −1 0 2 P1 x − x = , and P2 x − x = = x−y 1 x−y 2 2 6 You’ll notice that in the case of P2 , P2 x − x = (P2 − I)x is orthogonal to P2 x. Definition 0.2.2. P is said to be an orthogonal projector if it is a projector, and the range of P is orthogonal to the range of (I − P). We can show orthogonality by taking an arbitrary point in the range, Px and an arbitrary point in (I − P), (I − P)y, and show the dot product is 0. There is a property of real projectors that make them nice to work with: They are also symmetric matrices: Theorem 0.2.1. The (real) projector P is an orthogonal projector iff P = PT . For a proof, see for example, [?]. Caution: An orthogonal projector need not be an orthogonal matrix. Notice that the projector P2 from Figure 1 was not an orthogonal matrix (that is, P2 P2T 6= I). We have two primary sources for projectors: Projecting to a vector: Let a be an arbitrary, real, non-zero vector. We show that aaT Pa = kak2 is a rank one orthogonal projector onto the span of a: • The matrix aaT has rank one, since every column is a multiple of a. • The given matrix is a projector: P2 = aaT aaT 1 aaT · = a(aT a)aT = =P 2 2 4 kak kak kak kak2 • The matrix is an orthogonal projector, since PT = P. Projecting to a Subspace: Let Q = [q 1 , q 2 , . . . , q k ] be a matrix with orthonormal columns. Then P = QQT is an orthogonal projector to the column space of Q. This generalizes the result of the previous exercise. Note that if Q was additionally a square matrix, QQT = I. Note that this is exactly the property that we discussed in the last example of the previous section. Exercises 1. Show that the plane H defined by:       1 1   H = α1  1  + α2  −1  such that α1 , α2 ∈ IR   1 0 0.3. THE FOUR FUNDAMENTAL SUBSPACES 7 is isormorphic to IR2 . 2. Let the subspace G be the plane defined below, and consider the vector x, where:         1 1 3   G = α1  3  + α2  −1  such that α1 , α2 ∈ IR x= 0    −2 0 2 (a) Find the projector P that takes an arbitrary vector and projects it (orthogonally) to the plane G. (b) Find the orthogonal projection of the given x onto the plane G. (c) Find the distance from the plane G to the vector x. 3. If the low dimensional representation of a vector x is [9, −1]T and the basis vectors are [1, 0, 1]T and [3, 1, 1]T , then what was the original vector x? (HINT: it is easy to compute it directly) 4. If the vector x = [10, 4, 2]T and the basis vectors are [1, 0, 1]T and [3, 1, 1]T , then what is the low dimensional representation for x? 5. Let a = [−1, 3]T . Find a square matrix P so that P x is the orthogonal projection of x onto the span of a. 0.3 The Four Fundamental Subspaces Given any m × n matrix A, we consider the mapping A : IRn → IRm by: x → Ax = y The four subspaces allow us to completely understand the domain and range of the mapping. We will first define them, then look at some examples. Definition 0.3.1. The Four Fundamental Subspaces • The row space of A is a subspace of IRn formed by taking all possible linear combinations of the rows of A. Formally,  Row(A) = x ∈ IRn | x = AT y y ∈ IRm • The null space of A is a subspace of IRn formed by Null(A) = {x ∈ IRn | Ax = 0} • The column space of A is a subspace of IRm formed by taking all possible linear combinations of the columns of A. Col(A) = {y ∈ IRm | y = Ax ∈ IRn } The column space is also the image of the mapping. Notice that Ax is simply a linear combination of the columns of A: Ax = x1 a1 + x2 a2 + · · · + xn an 8 • Finally, we define the null space of AT can be defined in the obvious way (see the Exercises). The fundamental subspaces subdivide the domain and range of the mapping in a particularly nice way: Theorem 0.3.1. Let A be an m × n matrix. Then • The nullspace of A is orthogonal to the row space of A • The nullspace of AT is orthogonal to the columnspace of A Proof: See the Exercises. Before going further, let us recall how to construct a basis for the column space, row space and nullspace of a matrix A. We’ll do it with a particular matrix: Example 0.3.1. Construct a basis for the column space, row space and nullspace of the matrix A below:   2 0 −2 2 5 7 3  A =  −2 3 −5 −8 −2 SOLUTION: The row reduced form of A is:   1 0 −1 1  0 1 1 1  0 0 0 0 The first two columns of the original matrix form a basis for the columnspace (which is a subspace of IR3 ):     2  2  Col(A) = span  −2  ,  −2    3 3 A basis for the row space is found by using the row reduced rows corresponding to the pivots (and is a subspace of IR4 ). You should also verify that you can find a basis for the null space of A, given below (also a subspace of IR4 ). If you’re having any difficulties here, be sure to look it up in a linear algebra text:         1 0  1 −1            0   1   −1   −1  ,   ,  Row(A) = span  Null(A) = span    1   1   0     −1        1 1 0 1 We will often refer to the dimensions of the four subspaces. We recall that there is a term for the dimension of the column space- That is, the rank. Definition 0.3.2. The rank of a matrix A is the number of independent columns of A. 0.4. EXERCISES 9 In our previous example, the rank of A is 2. Also from our example, we see that the rank is the dimension of the column space, and that this is the same as the dimension of the row space (all three numbers correspond to the number of pivots in the row reduced form of A). Finally, a handy theorem for counting is the following. The Rank Theorem. Let the m × n matrix A have rank r. Then r + dim (Null(A)) = n This theorem says that the number of pivot columns plus the other columns (which correspond to free variables) is equal to the total number of columns. Example 0.3.2. The Dimensions of the Subspaces. Given a matrix A that is m × n with rank k, then the dimensions of the four subspaces are shown below. • dim (Row(A)) = k • dim (Col(A)) = k
• dim Null(AT ) = m − k • dim (Null(A)) = n − k There are some interesting implications of these theorems to matrices of data- For example, suppose A is m × n. With no other information, we do not know whether we should consider this matrix as n points in IRm , or m points in IRn . In one sense, it doesn’t matter! The theorems we’ve discussed shows that the dimension of the columnspace is equal to the dimension of the rowspace. Later on, we’ll find out that if we can find a basis for the columnspace, it is easy to find a basis for the rowspace. We’ll need some more machinery first. The Best Approximation Theorem If W is a subspace of IRn and x ∈ n IR , then the point closest to x in W is the orthogonal projection of x into W . We prove this in the exercises below. 0.4 Exercises 1. Show that Null(A) ⊥ Row(A). Hint: One way to show this is to take an arbitrary x1 6= ~0 ∈ Null(A) and show that it is orthogonal to every row of A. 2. If A is m × n, how big can the rank of A possibly be? 3. Show that multiplication by an orthogonal matrix preserves lengths: kQxk2 = kxk2 (Hint: Use properties of inner products). Conclude that multiplication by Q represents a rigid rotation. 4. Prove the Pythagorean Theorem by induction: Given a set of n orthogonal vectors {xi } n n X X 2 k xi k2 = kxi k22 i=1 i=1 10 5. Let A be an m × n matrix where m > n, and let A have rank n. Let x, y ∈ IRm , such that y is the orthogonal projection of x onto the columnspace of A. We want a formula for the projector P : IRm → IRm so that Px = y. (a) Why is the projector not P = AAT ? (b) Since y − x is orthogonal to the range of A, show that AT (y − x) = 0 (2) (c) Show that there exists v so that Equation (2) can be written as: AT (Av − x) = 0 (3) (d) Argue that AT A is invertible, so that Equation (3) implies that v = AT A
−1 AT x (e) Finally, show that this implies that P = A AT A
−1 AT Note: If A has rank k < m, then we will need something different, since AT A will not be full rank. The missing piece is the singular value decomposition, to be discussed later. 6. The Orthogonal Decomposition Theorem: if x ∈ IRn and W is a (nonzero) subspace of IRn , then x can be written uniquely as x=w+z where w ∈ W and z ∈ W ⊥ . p To prove this, let {ui }i=1 be an orthonormal basis for W , define w = (x, u1 )u1 + . . . + (x, up )up , and define z = x − w. Then: (a) Show that z ∈ W ⊥ by showing that it is orthogonal to every ui . (b) To show that the decomposition is unique, suppose it is not. That is, there are two decompositions: x = w1 + z 1 , x = w2 + z 2 Show this implies that w1 − w2 = z 2 − z 1 , and that this vector is in both W and W ⊥ . What can we conclude from this? 7. Use the previous exercises to prove the The Best Approximation Theorem If W is a subspace of IRn and x ∈ IRn , then the point closest to x in W is the orthogonal projection of x into W . 0.5. THE DECOMPOSITION THEOREMS 0.5 11 The Decomposition Theorems 0.5.1 The Eigenvector/Eigenvalue Decomposition 1. Definition: Let A be an n × n matrix. Then an eigenvector-eigenvalue pair is a vector v 6= 0, and a scalar λ where Av = λv ⇒ (A − λI)v = 0 (4) 2. Remark: If Equation (4) has a nontrivial solution, then det (A − λI) = 0 which leads to solving for the roots of a polynomial of degree n. This polynomial is called the characteristic polynomial. 3. Remark: We solve for the eigenvalues first, then solve for the nullspace of (A − λI) (which is called the eigenspace, Eλ ), by solving (A − λI)v = 0 4. Remark: Note that it is possible that one eigenvalue is repeated. This may or may not correspond with the same number of eigenvectors. 5. Definition: If eigenvalue λ is repeated k times as a root to the characteristic equation, then the algebraic multiplicity of λ is k. 6. Definition: If the eigenspace Eλ has dimension k, then λ has geometric multiplicity k. 7. Example: Compute the eigenvalues and eigenvectors for: (i) the 2 × 2 identity matrix, (ii) The matrix (in Matlab notation): [1 2;0 1] 8. Theorem: If aλ is the algebraic multiplicity of λ and gλ is the geometric multiplicity, then aλ ≥ gλ We won’t prove this here. 9. Definition: If, for some eigenvalue λ of A, we have that aλ > gλ , A is said to be defective. 10. Theorem: If A is square, and P is square and invertible, then A and P −1 AP have the same eigenvalues. 11. Exercise: Prove the previous theorem. 12. Remark: One method of characterizing eigenvalues in terms of the determinant and trace of a matrix: det(A) = n Y i=1 λi trace(A) = ∞ X i=1 λi 12 13. Remark: We will be especially interested in symmetric matrices. The rest of this section is devoted to them. 14. Definition: A matrix A is orthogonally diagonalizeable if there is an orthogonal matrix Q and diagonal matrix D so that so that A = QDQT . 15. The Spectral Theorem: If A is an n × n symmetric matrix, then: (a) A has n real eigenvalues (counting multiplicity). (b) For all λ, aλ = gλ . (c) The eigenspaces for distinct eigenvalues are mutually orthogonal. (d) A is orthogonally diagonalizeable, with D = diag(λ1 , λ2 , . . . , λn ). Some remarks about the Spectral Theorem: • We assume that inside each eigenspace, we have an orthonormal basis of eigenvectors. This is not a restriction, since we can always construct such a basis using Gram-Schmidt. • If a matrix is real and symmetric, the Spectral Theorem says that its eigenvectors form an orthonormal basis for IRn . • The full proof is beyond the scope of this course, but we can prove some parts (given below). 16. The following is a proof of part (c). Supply justification for each step: Let v 1 , v 2 be eigenvectors from distinct eigenvalues, λ1 , λ2 . We show that v 1 · v 2 = 0: λ1 v 1 · v 2 = (Av 1 )T v 2 = v T1 AT v 2 = v T1 Av 2 = λ2 v 1 · v 2 Now, (λ1 − λ2 )v 1 · v 2 = 0. 17. The Spectral Decomposition: Since A is orthogonally diagonalizable, then   T  q1 λ1 0 . . . 0  0 λ2 . . . 0   q T2     A = [q 1 q 2 . . . q n ]  . .. . . ..   ..   ..   . . . .  0 so that: 0 ... q Tn λn    A = [λ1 q 1 λ2 q 2 . . . λn q n ]   q T1 q T2 .. . q Tn      so finally: A = λ1 q 1 q T1 + λ2 q 2 q T2 + . . . + λn q n q Tn That is, A is a sum of n rank one matrices, each of which is a projection matrix. 0.5. THE DECOMPOSITION THEOREMS 13 18. Matlab Exercise: Verify the spectral decomposition for a symmetric matrix. Type the following into Matlab (the lines that begin with a % denote comments that do not have to be typed in). %Construct a random, symmetric, 6 x 6 matrix: for i=1:6 for j=1:i A(i,j)=rand; A(j,i)=A(i,j); end end %Compute the eigenvalues of A: [Q,L]=eig(A); %NOTE: A = Q L Q’ %L is a diagonal matrix %Now form the spectral sum S=zeros(6,6); for i=1:6 S=S+L(i,i)*Q(:,i)*Q(:,i)’; end max(max(S-A)) Note that the maximum of S − A should be a very small number! (By the spectral decomposition theorem). 0.5.2 The Singular Value Decomposition There is a special matrix factorization that is extremely useful, both in applications and in proving theorems. This is mainly due to two facts, which we shall investigate in this section: (1) We can use this factorization on any matrix, (2) The factorization defines explicitly the rank of the matrix, and all four matrix subspaces. In what follows, assume that A is an m × n matrix (so A maps IRn to IRm ). 1. Remark: Although A itself is not symmetric, AT A is n × n and symn metric. Therefore, it is orthogonally diagonalizable. Let {λi }i=1 and V = [v 1 , v 2 , . . . , v n ] be the eigenvalues and orthonormal eigenvectors. 2. Exercise: Show that λi ≥ 0 for i = 1..n by showing that kAv i k22 = λi . 3. Definition: We define the singular values of A by: p σi = λi where λi is an eigenvalue of AT A. 14 4. Remark: In the rest of the section, we will assume any list (or diagonal matrix) of eigenvalues of AT A (or singular values of A) will be ordered from highest to lowest: λ1 ≥ λ2 ≥ . . . ≥ λn . 5. Exercise: Let vi , vj be non-zero vectors for λi , λj . We note that, if λi 6= λj , then the eigenvectors are orthogonal by the Spectral Theorem. If λi = λj , it is possible to choose the vectors so we have an orthonormal basis for Eλ . Prove that, if v i and v j are distinct eigenvectors of AT A, then their corresponding images, Av i and Av j , are orthogonal. SOLUTION: (Avi )T Avj = v i AT Av j = λj v Ti v j = 0 6. Exercise: Prove that, if x = α1 v 1 + . . . αn v n , then kAxk2 = α12 λ1 + . . . + αn2 λn where λi , vi are the evals/evecs of AT A. SOLUTION: Just work it out using the dot product and orthogonality (or the Pythagorean Theorem), and Exercises 2 and 5. kAxk2 = kα1 Av 1 + . . . αn Av n k2 = α12 kAv 1 k2 + . . . αn2 kv n k2 n 7. Exercise: Let W be the subspace generated by the basis {v j }j=k+1 , where v j are the eigenvectors associated with the zero eigenvalues of AT A (therefore, we are assuming that the first k eigenvalues are NOT zero). Show that W = Null(A). SOLUTION: We show that W is the nullspace of A: Let x 6= 0 ∈ W . Then kAxk = xT AT Ax = 0 since AT Ax = 0. Therefore, Ax = 0, and x is in the nullspace of A. Since x was chosen arbitrarily, W is contained in nullspace of A. Now, if we take a vector in the null space of A, is it possible that it is NOT contained in W ? If x ∈ null space, then kAxk = 0. Furthermore, we have the expansion x = α1 v 1 + . . . + αr v r and therefore kAxk2 = α12 λ1 + . . . + αr2 λr which is zero only if every α is zero. Thus, W is the nullspace of A. 0.5. THE DECOMPOSITION THEOREMS 15 8. Exercise: Prove that if the rank of AT A is r, then so is the rank of A. SOLUTION: If the rank of AT A is r, the dimension of the nullspace of A is n − r by our previous exercise. By the rank theorem, rank(A) + dim(Null(A)) = n Thus, the rank of A is also r. 9. Remark: Define ui = 1 1 Av i = Av i kAv i k2 σi and let U be the matrix whose ith column is ui . 10. Remark: This definition only makes sense sense for the first r vectors v (otherwise, Av i = 0). Thus, we’ll have to extend the basis to span all of IRm . 11. Exercise: Sketch how you might do this. 12. Exercise: Show that ui is an eigenvector of AAT whose eigenvalue is also λi . SOLUTION: Show that AAT ui = λi ui : AAT ui = AAT  so that: AAT ui = 1 Av i σi  = 1 A(AT Av i ) σi λi λi Av i = σi ui σi σi 13. Exercise: Show that AT ui = σi v i SOLUTION: We show this for non-zero σi : T A ui = A T  1 Av i σi  = λi v i = σi v i σi 14. Remark: So far, we have shown how to construct two matrices, U and V given a matrix A. That is, the matrix V is constructed by the eigenvectors of AT A, and the matrix U can be constructed using the v’s or by finding the eigenvectors of AAT . 15. Exercise: Let A be m × n. Define the m × n matrix Σ = diag(σ1 , . . . , σn ) where σi is the ith singular value of the matrix A. Show that AV = U Σ 16 Figure 2: The geometric meaning of the right and left singular vectors of the SVD decomposition. Note that Av i = σi ui . The mapping x → Ax will map the unit circle on the left to the ellipse on the right. 16. The Singular Value Decomposition (SVD) Let A be any m × n matrix of rank r. Then A = U ΣV T where U, Σ, V are the matrices defined in the previous exercises. That is, U is an orthogonal m × m matrix, Σ is a diagonal m × n matrix, and V is an orthogonal n × n matrix. The u’s are called the left singular vectors and the v’s are called the right singular vectors. 17. Remark: Keep in mind the following relationship between the right and left singular vectors: Av i = σi ui T = σi v i A ui 18. Computing The Four Subspaces to a matrix A. Let A = U ΣV T be the SVD of A which has rank r. Be sure that the singular values are ordered from highest to lowest. Then: r (a) A basis for the columnspace of A, Col(A) is {ui }i=1 n (b) A basis for nullspace of A, Null(A) is {v i }i=r+1 r (c) A basis for the rowspace of A, Row(A) is {v i }i=1 m (d) A basis for the nullspace of AT , Null(AT ) is {ui }i=r+1 19. We can also visualize the right and left singular values as in Figure 2. We think of the v i as a special orthogonal basis in Rn (Domain) that maps to the ellipse whose axes are defined by σi ui . 20. The SVD is one of the premier tools of linear algebra, because it allows us to completely reveal everything we need to know about a matrix mapping: The rank, the basis of the nullspace, a basis for the column space, the basis 0.5. THE DECOMPOSITION THEOREMS Col(AT ) 17 v1 u1 v2 .. . vk u2 .. . ⇐⇒ Col(A) uk N(A) vk+1 .. . =⇒ 0 vn 0 uk+1 .. . ⇐= N(AT ) um IRn IRm Figure 3: The SVD of A ([U,S,V]=svd(A)) completely and explicitly describes the 4 fundamental subspaces associated with the matrix, as shown. We have a one to one correspondence between the rowspace and columnspace of A, the remaining v’s map to zero, and the remaining u’s map to zero (under AT ). for the nullspace of AT , and of the row space. This is depicted in Figure 3. 21. Lastly, the SVD provides a decomposition of any linear mapping into two “rotations” and a scaling. This will become important later when we try to deduce a mapping matrix from data (See the section on signal separation). 22. Exercise: Compute the SVD by hand of the following matrices:     0 2 1 1  0 0  0 0 0 0 SOLUTION: For the first matrix,    2 0 1 AAT = ⇒ U= 0 0 0 0 1   D= 2 0 0 0  18  AT A = 1 1 1 1  1 0 1 0  1 V =√ 2 ⇒  −1 1 1 1   D= 2 0 0 0  Therefore,   = 1 0 0 1  √ 2 0 0 0  1 √ 2  1 1 −1 1  23. Remark: If m or n is very large, it might not make sense to keep the full matrix U and V . 24. The Reduced SVD Let A be m × n with rank r. Then we can write: A = Ũ Σ̃V˜T where Ũ is an m × r matrix with orthogonal columns, Σ̃ is an r × r square matrix, and Ṽ is an n × r matrix. 25. Theorem: (Actually, this is just another way to express the SVD). Let A = U ΣV T be the SVD of A, which has rank r. Then: A= r X σi ui v Ti i=1 Therefore, we can approximate A by the sum of rank one matrices. 26. Matlab and the SVD Matlab has the SVD built in. The function specifications are: [U,S,V]=svd(A) and [U,S,V]=svd(A,0) where the first function call returns the full SVD, and the second call returns a reduced SVD- see Matlab’s help file for the details on the second call. 27. Matlab Exercise: Image Processing and the SVD. First, in Matlab, load the clown picture: load clown This loads a matrix X and a colormap, map, into the workspace. To see the clown, type: image(X); colormap(map) We now perform a Singular Value Decomposition on the clown. Type in: [U,S,V]=svd(X); How many vectors are needed to retain a good picture of the clown? Try performing a k−dimensional reconstruction of the image by typing: H=U(:,1:k)*S(1:k,1:k)*V(:,1:k)’; image(H) for k = 5, 10, 20 and 30. What do you see? 0.5. THE DECOMPOSITION THEOREMS 19 Generalized Inverses Let a matrix A be m × n with rank r. In the general case, A does not have an inverse. Is there a way of restricting the domain and range of the mapping y = Ax so that the map is invertible? We know that the columnspace and rowspace of A have the same dimensions. Therefore, there exists a 1-1 and onto map between these spaces, and this is our restriction. To “solve” y = Ax, we replace y by its orthogonal projection to the columnspace of A, ŷ. This gives the least squares solution, which makes the problem solvable. To get a unique solution, we replace x by its projection to the rowspace of A, x̂. The problem ŷ = Ax̂ now has a solution, and that solution is unique. We can rewrite this problem now in terms of the reduced SVD of A: x̂ = V V T x, ŷ = U U T y Now we can write: U U T y = U ΣV T V V T x
so that V Σ−1 U T y = V V T x (Exercise: Verify that these computations are correct!) Given an m × n matrix A, define its pseudoinverse, A† by: A† = V Σ−1 U T We have shown that the least squares solution to y = Ax is given by: x̂ = A† y where x̂ is in the rowspace of A, and its image, Ax̂ is the projection of y into the columnspace of A. Geometrically, we can understand these computations in terms of the four fundamental subspaces. T 6R(A ) Ax = y - 6R(A) qy - N (AT ) - N (A) IRn IRm 20 In this case, there is no value of x ∈ IRn which will map onto y, since y is outside the columnspace of A. To get a solution, we project y onto the columnspace of A as shown below: T 6R(A ) 6R(A) qy ŷ = U U T y - N (AT ) - N (A) Now it is possible to find an x that will map onto ŷ, but if the nullspace of A is nontrivial, then all of the points on the dotted line will also map to ŷ T 6R(A ) x̂ = V V T xr 6R(A) rx ŷ qy - N (AT ) - N (A) Finally, we must choose a unique value of x for the mapping- We choose the x inside the rowspace of A. This is a very useful idea, and it is one we will explore in more detail later. For now, notice that to get this solution, we analyzed our four fundamental subspaces in terms of the basis vectors given by the SVD. Exercises 1. Consider  2 1 3 1 −1 2      x1  x2  = 5 1 x3 (a) Before solving this problem, what are the dimensions of the four fundamental subspaces? (b) Use Matlab to compute the SVD of the matrix A, and solve the problem by computing the pseudoinverse of A directly. 0.5. THE DECOMPOSITION THEOREMS 21 (c) Check your answer explicitly and verify that x̂ and ŷ are in the rowspace and columnspace. (Hint: If a vector x is already in the rowspace, what is V V T x?) 2. Consider       2 −1 7 −3 4  1 −1 3  x1 0 1 −2    x2  2 −5 12    x3  −2 0 −4 x4 1 −3 7      =    5 1 0 −2 6       (a) Find the dimensions of the four fundamental subspaces by using the SVD of A (in Matlab). (b) Solve the problem. (c) Check your answer explicitly and verify that x̂ and ŷ are in the rowspace and columnspace. 3. Write the following in Matlab to reproduce Figure 2: theta=linspace(0,2*pi,30); z=exp(i*theta); X=[real(z);imag(z)]; %The domain points m=1/sqrt(2); A=(m*[1,1;1,-1])*[1,0;0,3]; Y=A*X; %The image of the circle t=linspace(0,1); vec1=[0;0]*(1-t)+[0;1]*t; vec2=[0;0]*(1-t)+[1;0]*t; %The basis vectors v Avec1=A*vec1; Avec2=A*vec2; %Image of the basis vectors figure(1) %The domain plot(X(1,:),X(2,:),’k’,vec1(1,:),vec1(2,:),’k’, vec2(1,:),vec2(2,:),’k’); axis equal figure(2) %The image plot(Y(1,:),Y(2,:),’k’,Avec1(1,:),Avec1(2,:),’k’, Avec2(1,:),Avec2(2,:),’k’); axis equal 4. In the previous example, what was the matrix A? The vectors v? The vectors u? The singular values σ1 , σ2 ? Once you’ve written these down from the program, perform the SVD of A in Matlab. Are the vectors the same that you wrote down? NOTE: These show that the singular vectors are not unique- they vary by ±v, or ±u. 22 0.6 Interactions Between Subspaces and the SVD Suppose that a matrix A is p×n and B is q ×n. Then we have four fundamental subspaces for each of A, B. In particular, the row spaces and null spaces of A and B are all in IRn . Note that this interpretation is looking at the matrix A as containing p sample data points from IRn , and the matrix B as containing q sample data points from IRn . It is natural to ask about the interaction of these subspaces of IRn , which are enumerated below: • The rowspace of A separate from B. This is also the intersection of the rowspace of A with the nullspace of B. • The rowspace of B separate from A. This is also the intersection of the rowspace of B with the nullspace of A. • The intersection of the rowspaces of A and B. • The intersection of the nullspaces of A and B. How can we obtain a basis for the intersection of nullspaces? This is in fact fairly easy as long as p, q are relatively small. We construct a new matrix Z that is p + q × n:   A Z= B and find the nullspace of this via the SVD of Z. Exercise: Prove that, if v ∈ IRn , and Zv = 0, then v ∈ Null(A) ∩ Null(B). Exercise: Prove directly that, if v ∈ IRn , and Zv 6= 0, then v ∈ Row(A) ∪ Row(B). Note that the first exercise proves the second in that: c (Null(A) ∩ Null(B)) = Null(A)c ∪ Null(B)c = Row(A) ∪ Row(B) where c is the complement of the set. We can find bases for the Row(A) and Row(B) directly- Let us define the two reduced SVDs: A = UA ΣA VAT , B = UB ΣB VBT Then the columns of VA , VB form a basis for the rowspace of A and the rowspace of B respectively. Before continuing further in the discussion, lets consider the idea of distance and angles between two subspaces. 0.6. INTERACTIONS BETWEEN SUBSPACES AND THE SVD 0.6.1 23 Angles Between Subspaces As in Golub and Van Loan, let F, G be two subspaces of IRm . Without loss of generality, assume: p = dim(F ) ≥ dim(G) = q ≥ 1   The principal angles θ1 , . . . , θq ∈ 0, π2 between F and G are defined recursively by: cos(θk ) = max max uT v = uTk vk u∈F v∈G subject to the additional constraints that u, v be unit length and they are orthogonal to the previously found vectors u1 , . . . , uk−1 and v1 , . . . vk−1 . The vectors u1 , . . . uq and v1 , . . . , vq are called the principal vectors between subspaces F and G. Some remarks about this definition: • Since u, v are normalized, the maximum value of uT v is 1, corresponding to θ = 0. Thus, principal vectors associated with this value of θ are the same. These vectors will also give the an orthonormal basis for the intersection of F, G. • Since θ is restricted, the smallest value of uT v is zero, corresponding to θ = π/2. The corresponding vectors of F and G will be orthogonal. • We can define the distance between subspaces F and G if p = q by using the largest principal angle, θp : dist(F, G) = q 1 − cos2 (θp ) = sin(θp ) For example, if F and G are the same subspace, then θk = 0 for k = 1, . . . , p, and the distance between them is 0. On the other extreme, if F and G are orthogonal, then θk = π/2 for k = 1, . . . , p and the distance between them is 1. 0.6.2 Computing the principal angles and vectors. Here we give an intuitive idea behind the algorithm; for a more details see Golub and Van Loan (Chapter 12). Given two sets of orthonormal basis vectors for subspaces F and G (we’ll use VA and VB found earlier), we can write: uT v = y T VAT VB z so that u = VA y and v = VB z. Thus, if U, V are the matrices whose columns are the principal vectors, then U T V = Y T (VAT VB )Z = diag(cos(θ1 ), . . . , cos(θq )) = D 24 Notice that this is the SVD of VAT VB : VAT VB = Y DZ T Note again that we have not proven anything- this was just an observation. For a proof, we would need additional facts about the SVD of a matrix that would take us too far afield. Example: Let F, G be subspaces of IR3 , where F is the x − y plane and G is the x − z plane. Clearly, there is a one-dimensional intersection. We show this using the computations in the previous discussion.         0  0   1  1 F = span  0  ,  1  , G = span  0  ,  0      0 0 0 1 We can use these basis vectors as VA , VB respectively. Now,       1 0 1 0 1 0 0 T   0 0 = VA VB = 0 0 0 1 0 0 1 Since this matrix is already in diagonal form,     1 0 cos(θ1 ) 0 =I IT 0 0 0 cos(θ2 ) so that Y = Z = I and U = VA I = VA , V = VB I = VB Therefore, u1 = v1 corresponding to the subspace intersection, and the distance between the subspaces is 1. We also see that there is no nontrivial intersection between the nullspaces of F, G. We would see this if we took the SVD of [VAT VBT ]T as suggested at the beginning of this section, since there would be no zero singular values. Similarly, we can find the intersection of F with G⊥ . In this case,     0   1 0 0  0 1 = VAT VB = = Y · 1 · 1 = Y · 1 · ZT 0 1 0 1 0 Since this “matrix” is already in diagonal form,     0 0 U = VA Y =  1  , V = VB · 1 =  1  0 0